Bugfix: 3D-LR projection (CGL2) bezier evaluation was wrong

This commit is contained in:
Kjetil Andre Johannessen 2017-09-29 10:44:24 +02:00 committed by Arne Morten Kvarving
parent 9329ac73f6
commit d3bab37cd1

View File

@ -34,7 +34,6 @@
#include "Vec3Oper.h"
#include <array>
ASMu3D::ASMu3D (unsigned char n_f)
: ASMunstruct(3,3,n_f), lrspline(nullptr), tensorspline(nullptr),
myGeoBasis(1), bezierExtract(myBezierExtract)
@ -1966,9 +1965,6 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
const int p1 = lrspline->order(0);
const int p2 = lrspline->order(1);
const int p3 = lrspline->order(2);
Go::BsplineBasis basis1 = getBezierBasis(p1);
Go::BsplineBasis basis2 = getBezierBasis(p2);
Go::BsplineBasis basis3 = getBezierBasis(p3);
Matrix B(p1*p2*p3, 4); // Bezier evaluation points and derivatives
// Evaluate the secondary solution field at each point
@ -1983,6 +1979,9 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
<<") not found."<< std::endl;
return false;
}
Go::BsplineBasis basis1 = getBezierBasis(p1, lrspline->getElement(iel)->umin(), lrspline->getElement(iel)->umax());
Go::BsplineBasis basis2 = getBezierBasis(p2, lrspline->getElement(iel)->vmin(), lrspline->getElement(iel)->vmax());
Go::BsplineBasis basis3 = getBezierBasis(p3, lrspline->getElement(iel)->wmin(), lrspline->getElement(iel)->wmax());
// Evaluate the basis functions at current parametric point
FiniteElement fe(lrspline->getElement(iel)->nBasisFunctions());
@ -1991,10 +1990,6 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
fe.w = gpar[2][i];
fe.iel = iel+1;
double du = lrspline->getElement(iel)->umax() - lrspline->getElement(iel)->umin();
double dv = lrspline->getElement(iel)->vmax() - lrspline->getElement(iel)->vmin();
double dw = lrspline->getElement(iel)->wmax() - lrspline->getElement(iel)->wmin();
double u[2*p1];
double v[2*p2];
double w[2*p3];
@ -2007,9 +2002,9 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
for(int jj=0; jj<p2; jj++) {
for(int ii=0; ii<p1; ii++, ib++) {
B(ib,1) = u[2*ii ]*v[2*jj ]*w[2*kk ];
B(ib,2) = u[2*ii+1]*v[2*jj ]*w[2*kk ]*2.0/du;
B(ib,3) = u[2*ii ]*v[2*jj+1]*w[2*kk ]*2.0/dv;
B(ib,4) = u[2*ii ]*v[2*jj ]*w[2*kk+1]*2.0/dw;
B(ib,2) = u[2*ii+1]*v[2*jj ]*w[2*kk ];
B(ib,3) = u[2*ii ]*v[2*jj+1]*w[2*kk ];
B(ib,4) = u[2*ii ]*v[2*jj ]*w[2*kk+1];
}
}
}