From 5bb6bfec3b38e52091a7a4f29871a63442edc9ed Mon Sep 17 00:00:00 2001 From: rho Date: Thu, 23 Jun 2011 13:22:47 +0000 Subject: [PATCH] Corrected interpolation for rational splines (NURBS) git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1070 e10b68d5-8a6e-419e-a041-bce267b0401d --- src/ASM/SplineField2D.C | 33 +++++++++++++----------- src/ASM/SplineField3D.C | 49 +++++++++++++++++++----------------- src/ASM/SplineFields2D.C | 31 +++++++++++++++-------- src/ASM/SplineFields3D.C | 54 +++++++++++++++++++++++----------------- 4 files changed, 96 insertions(+), 71 deletions(-) diff --git a/src/ASM/SplineField2D.C b/src/ASM/SplineField2D.C index e9ba41b9..a3caf990 100644 --- a/src/ASM/SplineField2D.C +++ b/src/ASM/SplineField2D.C @@ -83,6 +83,7 @@ double SplineField2D::valueFE(const FiniteElement& fe) const register const double* val_ptr = &values[val_start_ix]; if (rational) { + double tempw; double w = 0.0; const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * kdim + dim; register const double* w_ptr = &(surf->rcoefs_begin()[w_start_ix]); @@ -90,17 +91,19 @@ double SplineField2D::valueFE(const FiniteElement& fe) const for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) { register const double bval_v = *bval_v_ptr; tempVal = 0.0; + tempw = 0.0; for (register double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) { register const double bval_u = *bval_u_ptr; - tempVal += bval_u * (*val_ptr++); + tempVal += bval_u * (*val_ptr++)*(*w_ptr); - w += bval_u * (*w_ptr); + tempw += bval_u * (*w_ptr); w_ptr += kdim; } val += tempVal * bval_v; val_ptr += unum - uorder; + w += tempw * bval_v; w_ptr += kdim * (unum - uorder); } @@ -207,31 +210,33 @@ bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const // Compute the tensor product value int coefind = uleft-uorder+1 + unum*(vleft-vorder+1); int derivs_plus1=derivs+1; + double weight = 1.0; for (int jj = 0; jj < vorder; ++jj) { int jjd=jj*(derivs_plus1); std::fill(temp.begin(), temp.end(), 0.0); for (int ii = 0; ii < uorder; ++ii) { int iid=ii*(derivs_plus1); - const double *val_p = &values[coefind]; - const double *co_p = &co[coefind*kdim] + dim; - int temp_ind = 0; - for (int vder = 0; vder < derivs_plus1; ++vder) - for (int uder = 0; uder < vder+1; ++uder) { - temp[temp_ind] += b0[iid+vder - uder]*(*val_p); - temp_ind += ndim; - } - - for (int dd = dim;dd < kdim;dd++,co_p += kdim) { + + if (rational) { int temp_ind = 1; + weight = co[coefind*kdim+dim]; for (int vder = 0; vder < derivs_plus1; ++vder) { for (int uder = 0; uder < vder+1; ++uder) { - temp[temp_ind] += b0[iid+vder - uder]*(*co_p); + temp[temp_ind] += b0[iid+vder - uder]*weight; temp_ind += ndim; } } } - + + const double *val_p = &values[coefind]; + int temp_ind = 0; + for (int vder = 0; vder < derivs_plus1; ++vder) + for (int uder = 0; uder < vder+1; ++uder) { + temp[temp_ind] += b0[iid+vder - uder]*(*val_p)*weight; + temp_ind += ndim; + } + coefind += 1; } diff --git a/src/ASM/SplineField3D.C b/src/ASM/SplineField3D.C index 14da8f41..2130fdaa 100644 --- a/src/ASM/SplineField3D.C +++ b/src/ASM/SplineField3D.C @@ -93,32 +93,35 @@ double SplineField3D::valueFE(const FiniteElement& fe) const double w = 1.0; if (rational) { + double tempw, tempw2; + w = 0.0; const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * kdim + dim; const double *w_ptr = &(vol->rcoefs_begin()[w_start_ix]); for (double* bval_w_ptr = Bw.begin(); bval_w_ptr != Bw.end(); ++bval_w_ptr) { const double bval_w = *bval_w_ptr; - temp = 0.0; + temp = tempw = 0.0; for (double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) { const double bval_v = *bval_v_ptr; - temp2 = 0.0; + temp2 = tempw2 = 0.0; for (double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) { const double bval_u = *bval_u_ptr; - temp2 += bval_u * (*val_ptr++); - w += bval_u*(*w_ptr); + temp2 += bval_u * (*val_ptr++) * (*w_ptr); + tempw2 += bval_u*(*w_ptr); w_ptr += kdim; } - temp += temp2 * bval_v; + temp += temp2 * bval_v; + tempw += tempw2 * bval_v; val_ptr += unum - uorder; w_ptr += kdim * (unum - uorder); } tempResult += temp * bval_w; - w *= bval_w; + w += tempw * bval_w; val_ptr += unum * (vnum - vorder); w_ptr += kdim * unum * (vnum - vorder); @@ -251,6 +254,7 @@ bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const int coefind = uleft-uorder+1 + unum*(vleft-vorder+1 + vnum*(wleft-worder+1)); int derivs_plus1=derivs+1; + double weight = 1.0; for (int k = 0; k < worder; ++k) { int kd=k*derivs_plus1; std::fill(temp.begin(), temp.end(), 0.0); @@ -261,32 +265,31 @@ bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const for (int i = 0; i < uorder; ++i) { int id=i*derivs_plus1; - const double *val_p = &values[coefind]; - const double *co_p = &co[coefind*kdim] + dim; - - int temp_ind = 0; - for (int wder = 0; wder <= derivs; ++wder) { - for (int vder = 0; vder <= wder; ++vder) { - for (int uder = 0; uder <= vder; ++uder) { - temp2[temp_ind] - += b0[id + wder - vder]*(*val_p); - temp_ind+=kdim; - } - } - } - - for (int d = dim;d < kdim;d++, co_p += kdim) { + + if (rational) { int temp_ind = 1; + weight = co[coefind*kdim+dim]; for (int wder = 0; wder <= derivs; ++wder) { for (int vder = 0; vder <= wder; ++vder) { for (int uder = 0; uder <= vder; ++uder) { - temp2[temp_ind] - += b0[id + wder - vder]*(*co_p); + temp2[temp_ind] += b0[id + wder - vder]*weight; temp_ind+=ndim; } } } } + + int temp_ind = 0; + const double *val_p = &values[coefind]; + for (int wder = 0; wder <= derivs; ++wder) { + for (int vder = 0; vder <= wder; ++vder) { + for (int uder = 0; uder <= vder; ++uder) { + temp2[temp_ind] += b0[id + wder - vder]*(*val_p)*weight; + temp_ind+=kdim; + } + } + } + coefind += 1; } diff --git a/src/ASM/SplineFields2D.C b/src/ASM/SplineFields2D.C index 259d6958..42438355 100644 --- a/src/ASM/SplineFields2D.C +++ b/src/ASM/SplineFields2D.C @@ -97,25 +97,31 @@ bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const std::fill(result.begin(), result.end(), double(0)); if (rational) { + double wtemp; double w = 0.0; + const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * kdim + dim; register const double* w_ptr = &(surf->rcoefs_begin()[w_start_ix]); for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) { register const double bval_v = *bval_v_ptr; std::fill(tempVal.begin(), tempVal.end(), 0); + wtemp = 0.0; for (register double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) { register const double bval_u = *bval_u_ptr; for (vtemp = tempVal.begin(); vtemp != tempVal.end(); ++vtemp) { - *vtemp += bval_u * (*val_ptr++); + *vtemp += bval_u * (*val_ptr++)*(*w_ptr); } - w += bval_u * (*w_ptr); + wtemp += bval_u * (*w_ptr); w_ptr += kdim; } + vtemp = tempVal.begin(); for (register double* v = result.begin(); v != result.end(); ++v) { *v += (*vtemp++) * bval_v; } + w += wtemp * bval_v; + val_ptr += ncomp * (unum - uorder); w_ptr += kdim * (unum - uorder); } @@ -216,7 +222,7 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const const std::vector::iterator co = rational ? surf->rcoefs_begin() : surf->coefs_begin(); int kdim = dim + (rational ? 1 : 0); int ndim = ncomp + (rational ? 1 : 0); - + // Make temporary storage for the basis values and a temporary // computation cache. Go::ScratchVect b0(uorder*(derivs+1)); @@ -243,28 +249,31 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const // Compute the tensor product value int coefind = uleft-uorder+1 + unum*(vleft-vorder+1); int derivs_plus1=derivs+1; + double weight = 1.0; for (int jj = 0; jj < vorder; ++jj) { int jjd=jj*(derivs_plus1); std::fill(temp.begin(), temp.end(), 0.0); for (int ii = 0; ii < uorder; ++ii) { int iid=ii*(derivs_plus1); - const double *val_p = &values[coefind*ncomp]; - const double *co_p = &co[coefind*kdim] + dim; - for (int dd = 0; dd < ncomp; ++dd,++val_p) { - int temp_ind=dd; + + if (rational) { + weight = co[coefind*kdim + dim]; + int temp_ind = ncomp; for (int vder = 0; vder < derivs_plus1; ++vder) { for (int uder = 0; uder < vder+1; ++uder) { - temp[temp_ind] += b0[iid+vder - uder]*(*val_p); + temp[temp_ind] += b0[iid+vder - uder]*weight; temp_ind += ndim; } } } - for (int dd = dim;dd < kdim;dd++,co_p += kdim) { - int temp_ind = ncomp; + + const double *val_p = &values[coefind*ncomp]; + for (int dd = 0; dd < ncomp; ++dd,++val_p) { + int temp_ind=dd; for (int vder = 0; vder < derivs_plus1; ++vder) { for (int uder = 0; uder < vder+1; ++uder) { - temp[temp_ind] += b0[iid+vder - uder]*(*co_p); + temp[temp_ind] += b0[iid+vder - uder]*(*val_p)*(weight); temp_ind += ndim; } } diff --git a/src/ASM/SplineFields3D.C b/src/ASM/SplineFields3D.C index 19fe03c6..32dc2e46 100644 --- a/src/ASM/SplineFields3D.C +++ b/src/ASM/SplineFields3D.C @@ -111,36 +111,44 @@ bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const double w = 1.0; if (rational) { + double tempw, tempw2; w = 0.0; + const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * kdim + dim; const double *w_ptr = &(vol->rcoefs_begin()[w_start_ix]); for (double* bval_w_ptr = Bw.begin(); bval_w_ptr != Bw.end(); ++bval_w_ptr) { const double bval_w = *bval_w_ptr; std::fill(tempPt.begin(), tempPt.end(), 0); + tempw = 0.0; for (double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) { const double bval_v = *bval_v_ptr; std::fill(tempPt2.begin(), tempPt2.end(), 0); + tempw2 = 0.0; for (double* bval_u_ptr = Bu.begin(); bval_u_ptr != Bu.end(); ++bval_u_ptr) { const double bval_u = *bval_u_ptr; for (vtemp = tempPt2.begin(); vtemp != tempPt2.end(); ++vtemp) { - *vtemp += bval_u * (*val_ptr++); + *vtemp += bval_u * (*val_ptr++)*(*w_ptr); } - w += bval_u*(*w_ptr); + tempw2 += bval_u*(*w_ptr); w_ptr += kdim; } + vtemp = tempPt2.begin(); for (double* v = tempPt.begin(); v != tempPt.end(); ++v) { *v += (*vtemp++) * bval_v; } + + tempw += tempw2*bval_v; val_ptr += ncomp * (unum - uorder); w_ptr += kdim * (unum - uorder); } + vtemp = tempPt.begin(); for (double* v = tempResult.begin(); v != tempResult.end(); ++v) { *v += (*vtemp++) * bval_w; } - w *= bval_w; + w += tempw * bval_w; val_ptr += ncomp * unum * (vnum - vorder); w_ptr += kdim * unum * (vnum - vorder); } @@ -291,7 +299,8 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const // Compute the tensor product value int coefind = uleft-uorder+1 + unum*(vleft-vorder+1 + vnum*(wleft-worder+1)); int derivs_plus1=derivs+1; - + + double weight = 1.0; for (int k = 0; k < worder; ++k) { int kd=k*derivs_plus1; std::fill(temp.begin(), temp.end(), 0.0); @@ -302,36 +311,35 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const for (int i = 0; i < uorder; ++i) { int id=i*derivs_plus1; + + if (rational) { + int temp_ind = ncomp; + weight = co[coefind*kdim+dim]; + + for (int wder = 0; wder <= derivs; ++wder) { + for (int vder = 0; vder <= wder; ++vder) { + for (int uder = 0; uder <= vder; ++uder) { + temp2[temp_ind] += b0[id + wder - vder]*weight; + temp_ind+=ndim; + } + } + } + } + const double *val_p = &values[coefind*ncomp]; - const double *co_p = &co[coefind*kdim] + dim; - for (int d = 0; d < ncomp; ++d,++val_p) { int temp_ind = d; for (int wder = 0; wder <= derivs; ++wder) { for (int vder = 0; vder <= wder; ++vder) { for (int uder = 0; uder <= vder; ++uder) { - temp2[temp_ind] - += b0[id + wder - vder]*(*val_p); + temp2[temp_ind] += b0[id + wder - vder]*(*val_p)*weight; temp_ind+=kdim; } } } } - - for (int d = dim;d < kdim;d++, co_p += kdim) { - int temp_ind = ncomp; - - for (int wder = 0; wder <= derivs; ++wder) { - for (int vder = 0; vder <= wder; ++vder) { - for (int uder = 0; uder <= vder; ++uder) { - temp2[temp_ind] - += b0[id + wder - vder]*(*co_p); - temp_ind+=ndim; - } - } - } - } + coefind += 1; } @@ -373,7 +381,7 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const Go::volume_ratder(&restemp[0], dim, derivs, &restemp2[0]); for (int i = 0; i < totpts; ++i) { for (int d = 0; d < ncomp; ++d) { - pts[i][d] = restemp2[i*ncomp + d]; + pts[i][d] = restemp2[i*ndim + d]; } } } else {