OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Corrected interpolation for rational splines (NURBS)

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1070 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
rho 2011-06-23 13:22:47 +00:00 committed by Knut Morten Okstad
parent 45bad7f88b
commit 5bb6bfec3b
4 changed files with 96 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
src/ASM/SplineField2D.C
View File

@ -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);
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]*weight;
temp_ind += ndim;
}
}
}
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[temp_ind] += b0[iid+vder - uder]*(*val_p)*weight;
temp_ind += ndim;
}
for (int dd = dim;dd < kdim;dd++,co_p += kdim) {
int temp_ind = 1;
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_ind += ndim;
}
}
}
coefind += 1;
}

45
src/ASM/SplineField3D.C
View File

@ -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;
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;
}

29
src/ASM/SplineFields2D.C
View File

@ -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);
}
@ -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;
}
}

50
src/ASM/SplineFields3D.C
View File

@ -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);
}
@ -292,6 +300,7 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& 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);
@ -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;
const double *val_p = &values[coefind*ncomp];
const double *co_p = &co[coefind*kdim] + dim;
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];
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 {