Skip to content
Snippets Groups Projects
Commit 3c62eecb authored by Thomas Toulorge's avatar Thomas Toulorge
Browse files

Fixed computation and optimization of NCJ quality measure

parent e70595d5
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Mesh/qualityMeasures.cpp
+ 267
− 254
View file @ 3c62eecb
...@@ -21,146 +21,83 @@ ...@@ -21,146 +21,83 @@
namespace { namespace {
// Compute unit vector // Compute unit vector and gradients w.r.t. x0, y0, z0
inline void unitVec(const double &x0, const double &y0, const double &z0, // Remark: gradients w.r.t. x1, y1, z1 not computed as they are just the opposite
const double &x1, const double &y1, const double &z1, inline void unitVecAndGrad(const SPoint3 &p0, const SPoint3 &p1,
double &x01n, double &y01n, double &z01n) SVector3 &vec, std::vector<SVector3> &grad)
{ {
const double x01 = x1-x0, y01 = y1-y0, z01 = z1-z0; vec = SVector3(p0, p1);
const double invL01 = 1./sqrt(x01*x01 + y01*y01 + z01*z01); const double n = vec.normalize(), invN = 1/n;
x01n = x01*invL01; y01n = y01*invL01; z01n = z01*invL01; grad[0] = invN*vec[0]*vec; grad[0][0] -= invN; // dv01/dx0
grad[1] = invN*vec[1]*vec; grad[1][1] -= invN; // dv01/dy0
grad[2] = invN*vec[2]*vec; grad[2][2] -= invN; // dv01/dz0
} }
// Compute unit vector and gradients w.r.t. x0, y0, z0, x1, y1, z1 // Given vectors (0, 1) and (0, 2), their gradients and opposite of gradients,
inline void unitVecAndGrad(const double &x0, const double &y0, const double &z0, // and the unit normal vector, compute NCJ from area of triangle defined by both
const double &x1, const double &y1, const double &z1, // vectors and gradients w.r.t. x0, y0, z0, x1, y1, z1, x2, y2, z2
fullMatrix<double> &res) inline void NCJAndGrad2D(const SVector3 &v01,
const std::vector<SVector3> &dv01dp0,
const std::vector<SVector3> &dv01dp1,
const SVector3 &v02, const std::vector<SVector3> &dv02dp0,
const std::vector<SVector3> &dv02dp2,
const SVector3 &normal,
double &NCJ, std::vector<double> &dNCJ)
{ {
const double x01 = x1-x0, y01 = y1-y0, z01 = z1-z0; const SVector3 dvndx0 = crossprod(v01, dv02dp0[0]) + crossprod(dv01dp0[0], v02); // v01 x dv02/dx0 + dv01/dx0 x v02
const double x01Sq = x01*x01, y01Sq = y01*y01, z01Sq = z01*z01; const SVector3 dvndy0 = crossprod(v01, dv02dp0[1]) + crossprod(dv01dp0[1], v02); // v01 x dv02/dy0 + dv01/dy0 x v02
const double invL01Sq = 1./(x01Sq+y01Sq+z01Sq), invL01 = sqrt(invL01Sq); const SVector3 dvndz0 = crossprod(v01, dv02dp0[2]) + crossprod(dv01dp0[2], v02); // v01 x dv02/dz0 + dv01/dz0 x v02
const double invL01Cu = invL01Sq*invL01; const SVector3 dvndx1 = crossprod(dv01dp1[0], v02); // dv01/dx1 x v02 (= -dv01/dx0 x v02)
const SVector3 dvndy1 = crossprod(dv01dp1[1], v02); // dv01/dy1 x v02 (= -dv01/dy0 x v02)
double &dx01ndx0 = res(0, 0), &dx01ndy0 = res(0, 1), &dx01ndz0 = res(0, 2), const SVector3 dvndz1 = crossprod(dv01dp1[2], v02); // dv01/dz1 x v02 (= -dv01/dz0 x v02)
&dx01ndx1 = res(0, 3), &dx01ndy1 = res(0, 4), &dx01ndz1 = res(0, 5), &x01n = res(0, 6); const SVector3 dvndx2 = crossprod(v01, dv02dp2[0]); // v01 x dv02/dx2 (= v01 x -dv02/dx0)
double &dy01ndx0 = res(1, 0), &dy01ndy0 = res(1, 1), &dy01ndz0 = res(1, 2), const SVector3 dvndy2 = crossprod(v01, dv02dp2[1]); // v01 x dv02/dy2 (= v01 x -dv02/dy0)
&dy01ndx1 = res(1, 3), &dy01ndy1 = res(1, 4), &dy01ndz1 = res(1, 5), &y01n = res(1, 6); const SVector3 dvndz2 = crossprod(v01, dv02dp2[2]); // v01 x dv02/dz2 (= v01 x -dv02/dz0)
double &dz01ndx0 = res(2, 0), &dz01ndy0 = res(2, 1), &dz01ndz0 = res(2, 2),
&dz01ndx1 = res(2, 3), &dz01ndy1 = res(2, 4), &dz01ndz1 = res(2, 5), &z01n = res(2, 6); SVector3 vn = crossprod(v01, v02);
NCJ = dot(vn, normal); // NCJ
dx01ndx0 = x01Sq*invL01Cu-invL01; dNCJ[0] = dot(dvndx0, normal); // dNCJ/dx0
dx01ndy0 = x01*y01*invL01Cu; dNCJ[1] = dot(dvndy0, normal); // dNCJ/dy0
dx01ndz0 = x01*z01*invL01Cu; dNCJ[2] = dot(dvndz0, normal); // dNCJ/dz0
dx01ndx1 = -dx01ndx0; dNCJ[3] = dot(dvndx1, normal); // dNCJ/dx1
dx01ndy1 = -dx01ndy0; dNCJ[4] = dot(dvndy1, normal); // dNCJ/dy1
dx01ndz1 = -dx01ndz0; dNCJ[5] = dot(dvndz1, normal); // dNCJ/dz1
x01n = x01*invL01; dNCJ[6] = dot(dvndx2, normal); // dNCJ/dx2
dNCJ[7] = dot(dvndy2, normal); // dNCJ/dy2
dy01ndx0 = dx01ndy0; dNCJ[8] = dot(dvndz2, normal); // dNCJ/dz2
dy01ndy0 = y01Sq*invL01Cu-invL01;
dy01ndz0 = y01*z01*invL01Cu;
dy01ndx1 = -dy01ndx0;
dy01ndy1 = -dy01ndy0;
dy01ndz1 = -dy01ndz0;
y01n = y01*invL01;
dz01ndx0 = dx01ndz0;
dz01ndy0 = dy01ndz0;
dz01ndz0 = z01Sq*invL01Cu-invL01;
dz01ndx1 = -dx01ndz0;
dz01ndy1 = -dz01ndy0;
dz01ndz1 = -dz01ndz0;
z01n = z01*invL01;
} }
// Given vectors a and b, compute area of triangle defined by both vectors (times a factor) //// Revert vector and gradients
inline double triArea(const double fact, //inline void revertVG(const fullMatrix<double> &vg, fullMatrix<double> &res)
const double &xa, const double &ya, const double &za, //{
const double &xb, const double &yb, const double &zb) // res(0, 0) = -vg(0, 3); res(0, 1) = -vg(0, 4); res(0, 2) = -vg(0, 5); res(0, 6) = -vg(0, 6);
{ // res(1, 0) = -vg(1, 3); res(1, 1) = -vg(1, 4); res(1, 2) = -vg(1, 5); res(1, 6) = -vg(1, 6);
const double zn = xa*yb - xb*ya, yn = za*xb - xa*zb, xn = ya*zb - za*yb; // res(2, 0) = -vg(2, 3); res(2, 1) = -vg(2, 4); res(2, 2) = -vg(2, 5); res(2, 6) = -vg(2, 6);
return fact * sqrt(xn*xn + yn*yn + zn*zn); //}
}
// Given vectors (0, 1) and (0, 2) and their gradients, compute area of triangle // Scatter the NCJ gradients at vertex iV w.r.t vertices i0, i1 and i2
// defined by both vectors and grad. w.r.t. x0, y0, z0, x1, y1, z1, x2, y2, z2 (times a factor) // in the vector of gradients for 2D element of nV vertices
inline void triAreaAndGrad(const double fact, const fullMatrix<double> &vga, template<int iV, int nV, int i0, int i1, int i2>
const fullMatrix<double> &vgb, fullVector<double> &res) inline void scatterNCJGrad(const std::vector<double> &dNCJi, std::vector<double> &dNCJ)
{ {
const double &dxadx0 = vga(0, 0), &dxady0 = vga(0, 1), &dxadz0 = vga(0, 2), dNCJ[(iV*nV+i0)*3] = dNCJi[0]; dNCJ[(iV*nV+i0)*3+1] = dNCJi[1]; dNCJ[(iV*nV+i0)*3+2] = dNCJi[2];
&dxadx1 = vga(0, 3), &dxady1 = vga(0, 4), &dxadz1 = vga(0, 5), &xa = vga(0, 6); dNCJ[(iV*nV+i1)*3] = dNCJi[3]; dNCJ[(iV*nV+i1)*3+1] = dNCJi[4]; dNCJ[(iV*nV+i1)*3+2] = dNCJi[5];
const double &dyadx0 = vga(1, 0), &dyady0 = vga(1, 1), &dyadz0 = vga(1, 2), dNCJ[(iV*nV+i2)*3] = dNCJi[6]; dNCJ[(iV*nV+i2)*3+1] = dNCJi[7]; dNCJ[(iV*nV+i2)*3+2] = dNCJi[8];
&dyadx1 = vga(1, 3), &dyady1 = vga(1, 4), &dyadz1 = vga(1, 5), &ya = vga(1, 6);
const double &dzadx0 = vga(2, 0), &dzady0 = vga(2, 1), &dzadz0 = vga(2, 2),
&dzadx1 = vga(2, 3), &dzady1 = vga(2, 4), &dzadz1 = vga(2, 5), &za = vga(2, 6);
const double &dxbdx0 = vgb(0, 0), &dxbdy0 = vgb(0, 1), &dxbdz0 = vgb(0, 2),
&dxbdx2 = vgb(0, 3), &dxbdy2 = vgb(0, 4), &dxbdz2 = vgb(0, 5), &xb = vgb(0, 6);
const double &dybdx0 = vgb(1, 0), &dybdy0 = vgb(1, 1), &dybdz0 = vgb(1, 2),
&dybdx2 = vgb(1, 3), &dybdy2 = vgb(1, 4), &dybdz2 = vgb(1, 5), &yb = vgb(1, 6);
const double &dzbdx0 = vgb(2, 0), &dzbdy0 = vgb(2, 1), &dzbdz0 = vgb(2, 2),
&dzbdx2 = vgb(2, 3), &dzbdy2 = vgb(2, 4), &dzbdz2 = vgb(2, 5), &zb = vgb(2, 6);
const double xn = ya*zb - za*yb;
const double dxndx0 = ya*dzbdx0 + dyadx0*zb - yb*dzadx0 - dybdx0*za;
const double dxndy0 = ya*dzbdy0 + dyady0*zb - yb*dzady0 - dybdy0*za;
const double dxndz0 = ya*dzbdz0 + dyadz0*zb - yb*dzadz0 - dybdz0*za;
const double dxndx1 = dyadx1*zb - yb*dzadx1;
const double dxndy1 = dyady1*zb - yb*dzady1;
const double dxndz1 = dyadz1*zb - yb*dzadz1;
const double dxndx2 = ya*dzbdx2 - dybdx2*za;
const double dxndy2 = ya*dzbdy2 - dybdy2*za;
const double dxndz2 = ya*dzbdz2 - dybdz2*za;
const double yn = za*xb - xa*zb;
const double dyndx0 = -xa*dzbdx0 - dxadx0*zb + xb*dzadx0 + dxbdx0*za;
const double dyndy0 = -xa*dzbdy0 - dxady0*zb + xb*dzady0 + dxbdy0*za;
const double dyndz0 = -xa*dzbdz0 - dxadz0*zb + xb*dzadz0 + dxbdz0*za;
const double dyndx1 = xb*dzadx1 - dxadx1*zb;
const double dyndy1 = xb*dzady1 - dxady1*zb;
const double dyndz1 = xb*dzadz1 - dxadz1*zb;
const double dyndx2 = dxbdx2*za - xa*dzbdx2;
const double dyndy2 = dxbdy2*za - xa*dzbdy2;
const double dyndz2 = dxbdz2*za - xa*dzbdz2;
const double zn = xa*yb - xb*ya;
const double dzndx0 = xa*dybdx0 + dxadx0*yb - xb*dyadx0 - dxbdx0*ya;
const double dzndy0 = xa*dybdy0 + dxady0*yb - xb*dyady0 - dxbdy0*ya;
const double dzndz0 = xa*dybdz0 + dxadz0*yb - xb*dyadz0 - dxbdz0*ya;
const double dzndx1 = dxadx1*yb - xb*dyadx1;
const double dzndy1 = dxady1*yb - xb*dyady1;
const double dzndz1 = dxadz1*yb - xb*dyadz1;
const double dzndx2 = xa*dybdx2 - dxbdx2*ya;
const double dzndy2 = xa*dybdy2 - dxbdy2*ya;
const double dzndz2 = xa*dybdz2 - dxbdz2*ya;
double &dNCJdx0 = res(0), &dNCJdy0 = res(1), &dNCJdz0 = res(2),
&dNCJdx1 = res(3), &dNCJdy1 = res(4), &dNCJdz1 = res(5),
&dNCJdx2 = res(6), &dNCJdy2 = res(7), &dNCJdz2 = res(8), &NCJ = res(9);
const double n = sqrt(xn*xn + yn*yn + zn*zn), invN = 1./n, fact2 = fact*invN;
dNCJdx0 = (dxndx0*xn + dyndx0*yn + dzndx0*zn) * fact2;
dNCJdy0 = (dxndy0*xn + dyndy0*yn + dzndy0*zn) * fact2;
dNCJdz0 = (dxndz0*xn + dyndz0*yn + dzndz0*zn) * fact2;
dNCJdx1 = (dxndx1*xn + dyndx1*yn + dzndx1*zn) * fact2;
dNCJdy1 = (dxndy1*xn + dyndy1*yn + dzndy1*zn) * fact2;
dNCJdz1 = (dxndz1*xn + dyndz1*yn + dzndz1*zn) * fact2;
dNCJdx2 = (dxndx2*xn + dyndx2*yn + dzndx2*zn) * fact2;
dNCJdy2 = (dxndy2*xn + dyndy2*yn + dzndy2*zn) * fact2;
dNCJdz2 = (dxndz2*xn + dyndz2*yn + dzndz2*zn) * fact2;
NCJ = fact * n;
} }
// Revert vector and gradients // Scatter the NCJ gradients at vertex iV w.r.t vertices i0, i1, i2 and i3
inline void revertVG(const fullMatrix<double> &vg, fullMatrix<double> &res) // in the vector of gradients for 3D element of nV vertices
template<int iV, int nV, int i0, int i1, int i2, int i3>
inline void scatterNCJGrad(const std::vector<double> &dNCJi, std::vector<double> &dNCJ)
{ {
res(0, 0) = -vg(0, 3); res(0, 1) = -vg(0, 4); res(0, 2) = -vg(0, 5); res(0, 6) = -vg(0, 6); dNCJ[iV*nV+i0*3] = dNCJi[0]; dNCJ[iV*nV+i0*3+1] = dNCJi[1]; dNCJ[iV*nV+i0*3+2] = dNCJi[2];
res(1, 0) = -vg(1, 3); res(1, 1) = -vg(1, 4); res(1, 2) = -vg(1, 5); res(1, 6) = -vg(1, 6); dNCJ[iV*nV+i1*3] = dNCJi[3]; dNCJ[iV*nV+i1*3+1] = dNCJi[4]; dNCJ[iV*nV+i1*3+2] = dNCJi[5];
res(2, 0) = -vg(2, 3); res(2, 1) = -vg(2, 4); res(2, 2) = -vg(2, 5); res(2, 6) = -vg(2, 6); dNCJ[iV*nV+i2*3] = dNCJi[6]; dNCJ[iV*nV+i2*3+1] = dNCJi[7]; dNCJ[iV*nV+i2*3+2] = dNCJi[8];
dNCJ[iV*nV+i2*3] = dNCJi[9]; dNCJ[iV*nV+i2*3+1] = dNCJi[10]; dNCJ[iV*nV+i2*3+2] = dNCJi[11];
} }
...@@ -297,82 +234,133 @@ double qmTriangle::angles(MTriangle *e) ...@@ -297,82 +234,133 @@ double qmTriangle::angles(MTriangle *e)
void qmTriangle::NCJRange(const MTriangle *el, double &valMin, double &valMax) void qmTriangle::NCJRange(const MTriangle *el, double &valMin, double &valMax)
{ {
fullVector<double> ncj(3); const JacobianBasis *jac = el->getJacobianFuncSpace();
const MVertex *v0 = el->getVertex(0), *v1 = el->getVertex(1), *v2 = el->getVertex(2); fullMatrix<double> primNodesXYZ(3, 3);
NCJ(v0->x(), v0->y(), v0->z(), v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z(), ncj); // SVector3 geoNorm(0.,0.,0.);
valMin = *std::min_element(ncj.getDataPtr(), ncj.getDataPtr()+3); // std::map<MElement*,GEntity*>::const_iterator itEl2ent = element2entity.find(_el[iEl]);
valMax = *std::max_element(ncj.getDataPtr(), ncj.getDataPtr()+3); // GEntity *ge = (itEl2ent == element2entity.end()) ? 0 : itEl2ent->second;
// const bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization();
for (int i=0; i<jac->getNumPrimMapNodes(); i++) {
const MVertex *v = el->getVertex(i);
primNodesXYZ(i,0) = v->x();
primNodesXYZ(i,1) = v->y();
primNodesXYZ(i,2) = v->z();
// if (hasGeoNorm && (_vert[iV]->onWhat() == ge)) {
// double u, v;
// _vert[iV]->getParameter(0,u);
// _vert[iV]->getParameter(1,v);
// geoNorm += ((GFace*)ge)->normal(SPoint2(u,v));
// }
}
// if (hasGeoNorm && (geoNorm.normSq() == 0.)) {
// SPoint2 param = ((GFace*)ge)->parFromPoint(_el[iEl]->barycenter(true),false);
// geoNorm = ((GFace*)ge)->normal(param);
// }
fullMatrix<double> nM(1, 3);
jac->getPrimNormal2D(primNodesXYZ, nM);
SVector3 normal(nM(0, 0), nM(0, 1), nM(0, 2));
std::vector<double> ncj(3);
NCJ(el->getVertex(0)->point(), el->getVertex(1)->point(),
el->getVertex(2)->point(), normal, ncj);
valMin = *std::min_element(ncj.begin(), ncj.end());
valMax = *std::max_element(ncj.begin(), ncj.end());
} }
void qmTriangle::NCJ(const double &x0, const double &y0, const double &z0, void qmTriangle::NCJ(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
const double &x1, const double &y1, const double &z1, const SVector3 &normal, std::vector<double> &NCJ)
const double &x2, const double &y2, const double &z2,
fullVector<double> &ncj)
{ {
// Compute unit vectors for each edge // Compute unit vectors for each edge
double x01n, y01n, z01n, x12n, y12n, z12n, x20n, y20n, z20n; SVector3 v01n(p0, p1), v12n(p1, p2), v20n(p2, p0);
unitVec(x0, y0, z0, x1, y1, z1, x01n, y01n, z01n); v01n.normalize();
unitVec(x1, y1, z1, x2, y2, z2, x12n, y12n, z12n); v12n.normalize();
unitVec(x2, y2, z2, x0, y0, z0, x20n, y20n, z20n); v20n.normalize();
// Compute NCJ at vertex from unit vectors a and b as 0.5*||a^b||/A_equilateral // Compute NCJ at vertex from unit vectors a and b as 0.5*||a^b||/A_equilateral
// Factor = 2./sqrt(3.) = 0.5/A_equilateral // Factor = 2./sqrt(3.) = 0.5/A_equilateral
static const double fact = 2./sqrt(3.); static const double fact = 2./sqrt(3.);
ncj(0) = triArea(fact, x01n, y01n, z01n, -x20n, -y20n, -z20n); NCJ[0] = fact * dot(crossprod(v01n, -v20n), normal);
ncj(1) = triArea(fact, x12n, y12n, z12n, -x01n, -y01n, -z01n); NCJ[1] = fact * dot(crossprod(v12n, -v01n), normal);
ncj(2) = triArea(fact, x20n, y20n, z20n, -x12n, -y12n, -z12n); NCJ[2] = fact * dot(crossprod(v20n, -v12n), normal);
} }
void qmTriangle::NCJAndGradients(const double &x0, const double &y0, const double &z0, // Compute NCJ and its gradients at corners
const double &x1, const double &y1, const double &z1, // Gradients packed in vector: (dNCJ0/dx0, dNCJ0/dy0, dNCJ0/dz0,
const double &x2, const double &y2, const double &z2, // dNCJ0/dx1, ... dNCJ0/dz3, dNCJ1/dx0, ..., dNCJ3/dz3)
fullMatrix<double> &res) void qmTriangle::NCJAndGradients(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
const SVector3 &normal,
std::vector<double> &NCJ, std::vector<double> &dNCJ)
{ {
// Factor = 2./sqrt(3.) = 0.5/A_equilateral // Factor = 2./sqrt(3.) = 0.5/A_equilateral
static const double fact = 2./sqrt(3.); static const double fact = 2./sqrt(3.);
// Compute unit vector and its gradients for each edge // Compute unit vector, its gradients and opposite grandients for edge (0, 1)
fullMatrix<double> vg01(3,7); SVector3 v01n, v10n;
unitVecAndGrad(x0, y0, z0, x1, y1, z1, vg01); std::vector<SVector3> dv01ndp0(3), dv01ndp1(3);
fullMatrix<double> vg12(3,7); unitVecAndGrad(p0, p1, v01n, dv01ndp0);
unitVecAndGrad(x1, y1, z1, x2, y2, z2, vg12); v10n = -v01n;
fullMatrix<double> vg20(3,7); for (int i=0; i<3; i++) dv01ndp1[i] = -dv01ndp0[i];
unitVecAndGrad(x2, y2, z2, x0, y0, z0, vg20); const std::vector<SVector3> &dv10ndp1 = dv01ndp0, &dv10ndp0 = dv01ndp1;
fullMatrix<double> vg10(3,7);
revertVG(vg01, vg10); // Compute unit vector, its gradients and opposite grandients for edge (1, 2)
fullMatrix<double> vg21(3,7); SVector3 v12n, v21n;
revertVG(vg12, vg21); std::vector<SVector3> dv12ndp1(3), dv12ndp2(3);
fullMatrix<double> vg02(3,7); unitVecAndGrad(p1, p2, v12n, dv12ndp1);
revertVG(vg20, vg02); v21n = -v12n;
for (int i=0; i<3; i++) dv12ndp2[i] = -dv12ndp1[i];
const std::vector<SVector3> &dv21ndp2 = dv12ndp1, &dv21ndp1 = dv12ndp2;
// Compute unit vector, its gradients and opposite grandients for edge (2, 0)
SVector3 v20n, v02n;
std::vector<SVector3> dv20ndp2(3), dv20ndp0(3);
unitVecAndGrad(p2, p0, v20n, dv20ndp2);
v02n = -v20n;
for (int i=0; i<3; i++) dv20ndp0[i] = -dv20ndp2[i];
const std::vector<SVector3> &dv02ndp0 = dv20ndp2, &dv02ndp2 = dv20ndp0;
// Compute NCJ at vertex 0 as 0.5*||u01^u02||/A_triEqui // Compute NCJ at vertex 0 as 0.5*||u01^u02||/A_triEqui
// and gradients w.r.t. x0, y0, z0, x1, y1, z1, x2, y2, z2 // and gradients w.r.t. x0, y0, z0, x1, y1, z1, x2, y2, z2
fullVector<double> res0(10); std::vector<double> dNCJ0(9);
triAreaAndGrad(fact, vg01, vg02, res0); NCJAndGrad2D(v01n, dv01ndp0, dv01ndp1, v02n, dv02ndp0, dv02ndp2, normal, NCJ[0], dNCJ0);
res(0, 0) = res0(0); res(0, 1) = res0(1); res(0, 2) = res0(2); // dNCJ[0] = dNCJ0[0]; dNCJ[1] = dNCJ0[1]; dNCJ[2] = dNCJ0[2];
res(0, 3) = res0(3); res(0, 4) = res0(4); res(0, 5) = res0(5); // dNCJ[3] = dNCJ0[3]; dNCJ[4] = dNCJ0[4]; dNCJ[5] = dNCJ0[5];
res(0, 6) = res0(6); res(0, 7) = res0(7); res(0, 8) = res0(8); // dNCJ[6] = dNCJ0[6]; dNCJ[7] = dNCJ0[7]; dNCJ[8] = dNCJ0[8];
res(0, 9) = res0(9); scatterNCJGrad<0, 3, 0, 1, 2>(dNCJ0, dNCJ);
// Compute NCJ at vertex 1 as 0.5*||u12^u10||/A_triEqui // Compute NCJ at vertex 1 as 0.5*||u12^u10||/A_triEqui
// and gradients w.r.t. x1, y1, z1, x2, y2, z2, x0, y0, z0 // and gradients w.r.t. x1, y1, z1, x2, y2, z2, x0, y0, z0
fullVector<double> res1(10); std::vector<double> dNCJ1(9);
triAreaAndGrad(fact, vg12, vg10, res1); NCJAndGrad2D(v12n, dv12ndp1, dv12ndp2, v10n, dv10ndp1, dv10ndp0, normal, NCJ[1], dNCJ1);
res(1, 0) = res1(6); res(1, 1) = res1(7); res(1, 2) = res1(8); // dNCJ[9] = dNCJ1[6]; dNCJ[10] = dNCJ1[7]; dNCJ[11] = dNCJ1[8];
res(1, 3) = res1(0); res(1, 4) = res1(1); res(1, 5) = res1(2); // dNCJ[10] = dNCJ1[0]; dNCJ[11] = dNCJ1[1]; dNCJ[12] = dNCJ1[2];
res(1, 6) = res1(3); res(1, 7) = res1(4); res(1, 8) = res1(5); // dNCJ[13] = dNCJ1[3]; dNCJ[14] = dNCJ1[4]; dNCJ[15] = dNCJ1[5];
res(1, 9) = res1(9); scatterNCJGrad<1, 3, 1, 2, 0>(dNCJ1, dNCJ);
// Compute NCJ at vertex 2 as 0.5*||u20^u21||/A_triEqui // Compute NCJ at vertex 2 as 0.5*||u20^u21||/A_triEqui
// Compute NCJ at vertex 2 and gradients w.r.t. x2, y2, z2, x0, y0, z0, x1, y1, z1 // Compute NCJ at vertex 2 and gradients w.r.t. x2, y2, z2, x0, y0, z0, x1, y1, z1
fullVector<double> res2(10); std::vector<double> dNCJ2(9);
triAreaAndGrad(fact, vg20, vg21, res2); NCJAndGrad2D(v20n, dv20ndp2, dv20ndp0, v21n, dv21ndp2, dv21ndp1, normal, NCJ[2], dNCJ2);
res(2, 0) = res2(3); res(2, 1) = res2(4); res(2, 2) = res2(5); // dNCJ[16] = dNCJ2[3]; dNCJ[17] = dNCJ2[4]; dNCJ[18] = dNCJ2[5];
res(2, 3) = res2(6); res(2, 4) = res2(7); res(2, 5) = res2(8); // dNCJ[19] = dNCJ2[6]; dNCJ[20] = dNCJ2[7]; dNCJ[21] = dNCJ2[8];
res(2, 6) = res2(0); res(2, 7) = res2(1); res(2, 8) = res2(2); // dNCJ[22] = dNCJ2[0]; dNCJ[23] = dNCJ2[1]; dNCJ[24] = dNCJ2[2];
res(2, 9) = res2(9); scatterNCJGrad<2, 3, 2, 0, 1>(dNCJ2, dNCJ);
for (int i=0; i<3; i++) NCJ[i] *= fact;
for (int i=0; i<27; i++) dNCJ[i] *= fact;
// for (int iV=0; iV<3; iV++) {
// std::cout << "DBGTT: Vertex " << iV << ":\n";
// std::cout << "DBGTT: -> NCJ = " << NCJ[iV] << "\n";
// for (unsigned ig=0; ig<3; ig++) {
// int ind = iV*9+ig*3;
// std::cout << "DBGTT: -> dNCJ/dp" << ig << " = (" << dNCJ[ind] << ", " << dNCJ[ind+1] << ", " << dNCJ[ind+2] << ")\n";
//// int ind2 = ig*3;
//// std::vector<double> dNCJLoc = (iV == 0) ? dNCJ0 : (iV == 1) ? dNCJ1 : dNCJ2;
//// std::cout << "DBGTT: -> dNCJ/dp" << ig << " (local) = (" << dNCJLoc[ind2] << ", " << dNCJLoc[ind2+1] << ", " << dNCJLoc[ind2+2] << ")\n";
// }
// }
} }
...@@ -478,105 +466,130 @@ double qmQuadrangle::angles(MQuadrangle *e) ...@@ -478,105 +466,130 @@ double qmQuadrangle::angles(MQuadrangle *e)
void qmQuadrangle::NCJRange(const MQuadrangle *el, double &valMin, double &valMax) void qmQuadrangle::NCJRange(const MQuadrangle *el, double &valMin, double &valMax)
{ {
fullVector<double> ncj(4); const JacobianBasis *jac = el->getJacobianFuncSpace();
const MVertex *v0 = el->getVertex(0), *v1 = el->getVertex(1), fullMatrix<double> primNodesXYZ(4, 3);
*v2 = el->getVertex(2), *v3 = el->getVertex(3); // SVector3 geoNorm(0.,0.,0.);
NCJ(v0->x(), v0->y(), v0->z(), v1->x(), v1->y(), v1->z(), // std::map<MElement*,GEntity*>::const_iterator itEl2ent = element2entity.find(_el[iEl]);
v2->x(), v2->y(), v2->z(), v3->x(), v3->y(), v3->z(), ncj); // GEntity *ge = (itEl2ent == element2entity.end()) ? 0 : itEl2ent->second;
valMin = *std::min_element(ncj.getDataPtr(), ncj.getDataPtr()+4); // const bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization();
valMax = *std::max_element(ncj.getDataPtr(), ncj.getDataPtr()+4); for (int i=0; i<jac->getNumPrimMapNodes(); i++) {
const MVertex *v = el->getVertex(i);
primNodesXYZ(i,0) = v->x();
primNodesXYZ(i,1) = v->y();
primNodesXYZ(i,2) = v->z();
// if (hasGeoNorm && (_vert[iV]->onWhat() == ge)) {
// double u, v;
// _vert[iV]->getParameter(0,u);
// _vert[iV]->getParameter(1,v);
// geoNorm += ((GFace*)ge)->normal(SPoint2(u,v));
// }
}
// if (hasGeoNorm && (geoNorm.normSq() == 0.)) {
// SPoint2 param = ((GFace*)ge)->parFromPoint(_el[iEl]->barycenter(true),false);
// geoNorm = ((GFace*)ge)->normal(param);
// }
fullMatrix<double> nM(1, 3);
jac->getPrimNormal2D(primNodesXYZ, nM);
SVector3 normal(nM(0, 0), nM(0, 1), nM(0, 2));
std::vector<double> ncj(4);
NCJ(el->getVertex(0)->point(), el->getVertex(1)->point(),
el->getVertex(2)->point(), el->getVertex(3)->point(),
normal, ncj);
valMin = *std::min_element(ncj.begin(), ncj.end());
valMax = *std::max_element(ncj.begin(), ncj.end());
} }
void qmQuadrangle::NCJ(const double &x0, const double &y0, const double &z0, void qmQuadrangle::NCJ(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
const double &x1, const double &y1, const double &z1, const SPoint3 &p3, const SVector3 &normal, std::vector<double> &ncj)
const double &x2, const double &y2, const double &z2,
const double &x3, const double &y3, const double &z3,
fullVector<double> &ncj)
{ {
// Compute unit vectors for each edge // Compute unit vectors for each edge
double x01n, y01n, z01n, x12n, y12n, z12n, x23n, y23n, z23n, x30n, y30n, z30n; SVector3 v01n(p0, p1), v12n(p1, p2), v23n(p2, p3), v30n(p3, p0);
unitVec(x0, y0, z0, x1, y1, z1, x01n, y01n, z01n); v01n.normalize();
unitVec(x1, y1, z1, x2, y2, z2, x12n, y12n, z12n); v12n.normalize();
unitVec(x2, y2, z2, x3, y3, z3, x23n, y23n, z23n); v23n.normalize();
unitVec(x3, y3, z3, x0, y0, z0, x30n, y30n, z30n); v30n.normalize();
// Compute NCJ at vertex from unit vectors a and b as 0.5*||a^b||/A_triRect // Compute NCJ at vertex from unit vectors a and b as 0.5*||a^b||/A_equilateral
// Factor = 1. = 0.5/A_triRect ncj[0] = dot(crossprod(v01n, -v30n), normal);
static const double fact = 1.; ncj[1] = dot(crossprod(v12n, -v01n), normal);
ncj(0) = triArea(fact, x01n, y01n, z01n, -x30n, -y30n, -z30n); ncj[2] = dot(crossprod(v23n, -v12n), normal);
ncj(1) = triArea(fact, x12n, y12n, z12n, -x01n, -y01n, -z01n); ncj[3] = dot(crossprod(v30n, -v23n), normal);
ncj(2) = triArea(fact, x23n, y23n, z23n, -x12n, -y12n, -z12n);
ncj(3) = triArea(fact, x30n, y30n, z30n, -x23n, -y23n, -z23n);
} }
void qmQuadrangle::NCJAndGradients(const double &x0, const double &y0, const double &z0, void qmQuadrangle::NCJAndGradients(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
const double &x1, const double &y1, const double &z1, const SPoint3 &p3, const SVector3 &normal,
const double &x2, const double &y2, const double &z2, std::vector<double> &NCJ, std::vector<double> &dNCJ)
const double &x3, const double &y3, const double &z3,
fullMatrix<double> &res)
{ {
// Factor = 1. = 0.5/A_triRect // Compute unit vector, its gradients and opposite gradients for edge (0,1)
static const double fact = 1.; SVector3 v01n, v10n;
std::vector<SVector3> dv01ndp0(3), dv01ndp1(3);
// Compute unit vector and its gradients for each edge unitVecAndGrad(p0, p1, v01n, dv01ndp0);
fullMatrix<double> vg01(3,7); v10n = -v01n;
unitVecAndGrad(x0, y0, z0, x1, y1, z1, vg01); for (int i=0; i<3; i++) dv01ndp1[i] = -dv01ndp0[i];
fullMatrix<double> vg12(3,7); const std::vector<SVector3> &dv10ndp1 = dv01ndp0, &dv10ndp0 = dv01ndp1;
unitVecAndGrad(x1, y1, z1, x2, y2, z2, vg12);
fullMatrix<double> vg23(3,7); // Compute unit vector, its gradients and opposite gradients for edge (1,2)
unitVecAndGrad(x2, y2, z2, x3, y3, z3, vg23); SVector3 v12n, v21n;
fullMatrix<double> vg30(3,7); std::vector<SVector3> dv12ndp1(3), dv12ndp2(3);
unitVecAndGrad(x3, y3, z3, x0, y0, z0, vg30); unitVecAndGrad(p1, p2, v12n, dv12ndp1);
fullMatrix<double> vg10(3,7); v21n = -v12n;
revertVG(vg01, vg10); for (int i=0; i<3; i++) dv12ndp2[i] = -dv12ndp1[i];
fullMatrix<double> vg21(3,7); const std::vector<SVector3> &dv21ndp2 = dv12ndp1, &dv21ndp1 = dv12ndp2;
revertVG(vg12, vg21);
fullMatrix<double> vg32(3,7); // Compute unit vector, its gradients and opposite gradients for edge (2,3)
revertVG(vg23, vg32); SVector3 v23n, v32n;
fullMatrix<double> vg03(3,7); std::vector<SVector3> dv23ndp2(3), dv23ndp3(3);
revertVG(vg30, vg03); unitVecAndGrad(p2, p3, v23n, dv23ndp2);
v32n = -v23n;
for (int i=0; i<3; i++) dv23ndp3[i] = -dv23ndp2[i];
const std::vector<SVector3> &dv32ndp3 = dv23ndp2, &dv32ndp2 = dv23ndp3;
// Compute unit vector, its gradients and opposite gradients for edge (3,0)
SVector3 v30n, v03n;
std::vector<SVector3> dv30ndp3(3), dv30ndp0(3);
unitVecAndGrad(p3, p0, v30n, dv30ndp3);
v03n = -v30n;
for (int i=0; i<3; i++) dv30ndp0[i] = -dv30ndp3[i];
const std::vector<SVector3> &dv03ndp0 = dv30ndp3, &dv03ndp3 = dv30ndp0;
// Compute NCJ at vertex 0 as 0.5*||u01^u03||/A_triRect // Compute NCJ at vertex 0 as 0.5*||u01^u03||/A_triRect
// and gradients w.r.t. x0, y0, z0, x1, y1, z1, x3, y3, z3 // and gradients w.r.t. x0, y0, z0, x1, y1, z1, x3, y3, z3
fullVector<double> res0(10); std::vector<double> dNCJ0(9);
triAreaAndGrad(fact, vg01, vg03, res0); NCJAndGrad2D( v01n, dv01ndp0, dv01ndp1, v03n, dv03ndp0, dv03ndp3, normal, NCJ[0], dNCJ0);
res(0, 0) = res0(0); res(0, 1) = res0(1); res(0, 2) = res0(2); scatterNCJGrad<0, 4, 0, 1, 3>(dNCJ0, dNCJ);
res(0, 3) = res0(3); res(0, 4) = res0(4); res(0, 5) = res0(5);
res(0, 6) = 0.; res(0, 7) = 0.; res(0, 8) = 0.;
res(0, 9) = res0(6); res(0, 10) = res0(7); res(0, 11) = res0(8);
res(0, 12) = res0(9);
// Compute NCJ at vertex 1 as 0.5*||u12^u10||/A_triRect // Compute NCJ at vertex 1 as 0.5*||u12^u10||/A_triRect
// and gradients w.r.t. x1, y1, z1, x2, y2, z2, x0, y0, z0 // and gradients w.r.t. x1, y1, z1, x2, y2, z2, x0, y0, z0
fullVector<double> res1(10); std::vector<double> dNCJ1(9);
triAreaAndGrad(fact, vg12, vg10, res1); NCJAndGrad2D( v12n, dv12ndp1, dv12ndp2, v10n, dv10ndp1, dv10ndp0, normal, NCJ[1], dNCJ1);
res(1, 0) = res1(6); res(1, 1) = res1(7); res(1, 2) = res1(8); scatterNCJGrad<1, 4, 1, 2, 0>(dNCJ1, dNCJ);
res(1, 3) = res1(0); res(1, 4) = res1(1); res(1, 5) = res1(2);
res(1, 6) = res1(3); res(1, 7) = res1(4); res(1, 8) = res1(5);
res(1, 9) = 0.; res(1, 10) = 0.; res(1, 11) = 0.;
res(1, 12) = res1(9);
// Compute NCJ at vertex 2 as 0.5*||u23^u21||/A_triRect // Compute NCJ at vertex 2 as 0.5*||u23^u21||/A_triRect
// and gradients w.r.t. x2, y2, z2, x3, y3, z3, x1, y1, z1 // and gradients w.r.t. x2, y2, z2, x3, y3, z3, x1, y1, z1
fullVector<double> res2(10); std::vector<double> dNCJ2(9);
triAreaAndGrad(fact, vg23, vg21, res2); NCJAndGrad2D( v23n, dv23ndp2, dv23ndp3, v21n, dv21ndp2, dv21ndp1, normal, NCJ[2], dNCJ2);
res(2, 0) = 0.; res(2, 1) = 0.; res(2, 2) = 0.; scatterNCJGrad<2, 4, 2, 3, 1>(dNCJ2, dNCJ);
res(2, 3) = res2(6); res(2, 4) = res2(7); res(2, 5) = res2(8);
res(2, 6) = res2(0); res(2, 7) = res2(1); res(2, 8) = res2(2);
res(2, 9) = res2(3); res(2, 10) = res2(4); res(2, 11) = res2(5);
res(2, 12) = res2(9);
// Compute NCJ at vertex 3 as 0.5*||u30^u32||/A_triRect // Compute NCJ at vertex 3 as 0.5*||u30^u32||/A_triRect
// and gradients w.r.t. x3, y3, z3, x0, y0, z0, x2, y2, z2 // and gradients w.r.t. x3, y3, z3, x0, y0, z0, x2, y2, z2
fullVector<double> res3(10); std::vector<double> dNCJ3(9);
triAreaAndGrad(fact, vg23, vg21, res3); NCJAndGrad2D(v30n, dv30ndp3, dv30ndp0, v32n, dv32ndp3, dv32ndp2, normal, NCJ[3], dNCJ3);
res(3, 0) = res3(3); res(3, 1) = res3(4); res(3, 2) = res3(5); scatterNCJGrad<3, 4, 3, 0, 2>(dNCJ3, dNCJ);
res(3, 3) = 0.; res(3, 4) = 0.; res(3, 5) = 0.;
res(3, 6) = res3(6); res(3, 7) = res3(7); res(3, 8) = res3(8); // for (int iV=0; iV<4; iV++) {
res(3, 9) = res3(0); res(3, 10) = res3(1); res(3, 11) = res3(2); // std::cout << "DBGTT: Vertex " << iV << ":\n";
res(3, 12) = res3(9); // std::cout << "DBGTT: -> NCJ = " << NCJ[iV] << "\n";
// for (unsigned ig=0; ig<4; ig++) {
// int ind = iV*12+ig*3;
// std::cout << "DBGTT: -> dNCJ/dp" << ig << " = (" << dNCJ[ind] << ", " << dNCJ[ind+1] << ", " << dNCJ[ind+2] << ")\n";
//// int ind2 = ig*3;
//// std::vector<double> dNCJLoc = (iV == 0) ? dNCJ0 : (iV == 1) ? dNCJ1 : dNCJ2;
//// std::cout << "DBGTT: -> dNCJ/dp" << ig << " (local) = (" << dNCJLoc[ind2] << ", " << dNCJLoc[ind2+1] << ", " << dNCJLoc[ind2+2] << ")\n";
// }
// }
} }
......
Mesh/qualityMeasures.h
+ 25
− 29
View file @ 3c62eecb
...@@ -6,8 +6,12 @@ ...@@ -6,8 +6,12 @@
#ifndef _QUALITY_MEASURES_H_ #ifndef _QUALITY_MEASURES_H_
#define _QUALITY_MEASURES_H_ #define _QUALITY_MEASURES_H_
#include "fullMatrix.h" //#include "fullMatrix.h"
#include <vector>
//#include "SPoint3.h"
class SPoint3;
class SVector3;
class BDS_Point; class BDS_Point;
class BDS_Face; class BDS_Face;
class MVertex; class MVertex;
...@@ -35,14 +39,11 @@ public: ...@@ -35,14 +39,11 @@ public:
static double minNCJ(const MTriangle *e); static double minNCJ(const MTriangle *e);
static void NCJRange(const MTriangle *e, double &valMin, double &valMax); static void NCJRange(const MTriangle *e, double &valMin, double &valMax);
static inline int numNCJVal() { return 3; } static inline int numNCJVal() { return 3; }
static void NCJ(const double &x0, const double &y0, const double &z0, static void NCJ(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
const double &x1, const double &y1, const double &z1, const SVector3 &normal, std::vector<double> &NCJ);
const double &x2, const double &y2, const double &z2, static void NCJAndGradients(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
fullVector<double> &ncj); const SVector3 &normal,
static void NCJAndGradients(const double &x0, const double &y0, const double &z0, std::vector<double> &NCJ, std::vector<double> &dNCJ);
const double &x1, const double &y1, const double &z1,
const double &x2, const double &y2, const double &z2,
fullMatrix<double> &ncj);
}; };
...@@ -55,16 +56,11 @@ public: ...@@ -55,16 +56,11 @@ public:
static double minNCJ(const MQuadrangle *e); static double minNCJ(const MQuadrangle *e);
static void NCJRange(const MQuadrangle *e, double &valMin, double &valMax); static void NCJRange(const MQuadrangle *e, double &valMin, double &valMax);
static inline int numNCJVal() { return 4; } static inline int numNCJVal() { return 4; }
static void NCJ(const double &x0, const double &y0, const double &z0, static void NCJ(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
const double &x1, const double &y1, const double &z1, const SPoint3 &p3, const SVector3 &normal, std::vector<double> &ncj);
const double &x2, const double &y2, const double &z2, static void NCJAndGradients(const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
const double &x3, const double &y3, const double &z3, const SPoint3 &p3, const SVector3 &normal,
fullVector<double> &ncj); std::vector<double> &NCJ, std::vector<double> &dNCJ);
static void NCJAndGradients(const double &x0, const double &y0, const double &z0,
const double &x1, const double &y1, const double &z1,
const double &x2, const double &y2, const double &z2,
const double &x3, const double &y3, const double &z3,
fullMatrix<double> &ncj);
}; };
...@@ -99,16 +95,16 @@ public: ...@@ -99,16 +95,16 @@ public:
static double minNCJ(const MTetrahedron *e); static double minNCJ(const MTetrahedron *e);
// static void NCJRange(const MTetrahedron *e, double &valMin, double &valMax); // static void NCJRange(const MTetrahedron *e, double &valMin, double &valMax);
static inline int numNCJVal() { return 4; } static inline int numNCJVal() { return 4; }
static void NCJ(const double &x0, const double &y0, const double &z0, // static void NCJ(const double &x0, const double &y0, const double &z0,
const double &x1, const double &y1, const double &z1, // const double &x1, const double &y1, const double &z1,
const double &x2, const double &y2, const double &z2, // const double &x2, const double &y2, const double &z2,
const double &x3, const double &y3, const double &z3, // const double &x3, const double &y3, const double &z3,
fullVector<double> &ncj); // fullVector<double> &ncj);
static void NCJAndGradients(const double &x0, const double &y0, const double &z0, // static void NCJAndGradients(const double &x0, const double &y0, const double &z0,
const double &x1, const double &y1, const double &z1, // const double &x1, const double &y1, const double &z1,
const double &x2, const double &y2, const double &z2, // const double &x2, const double &y2, const double &z2,
const double &x3, const double &y3, const double &z3, // const double &x3, const double &y3, const double &z3,
fullMatrix<double> &ncj); // fullMatrix<double> &ncj);
}; };
......
contrib/MeshOptimizer/MeshOptPatch.cpp
+ 42
− 43
View file @ 3c62eecb
...@@ -269,7 +269,7 @@ void Patch::initScaledJac() ...@@ -269,7 +269,7 @@ void Patch::initScaledJac()
if ((_dim == 2) && _scaledNormEl.empty()) { if ((_dim == 2) && _scaledNormEl.empty()) {
_scaledNormEl.resize(nEl()); _scaledNormEl.resize(nEl());
// for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl); // for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(iEl); for (int iEl = 0; iEl < nEl(); iEl++) calcNormalEl2D(iEl, true);
} }
else if (_invStraightJac.empty()) { else if (_invStraightJac.empty()) {
_invStraightJac.resize(nEl(),1.); _invStraightJac.resize(nEl(),1.);
...@@ -293,7 +293,7 @@ void Patch::initMetricMin() ...@@ -293,7 +293,7 @@ void Patch::initMetricMin()
// TODO: Re-introduce normal to geometry? // TODO: Re-introduce normal to geometry?
//void Mesh::calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl) //void Mesh::calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl)
void Patch::calcScaledNormalEl2D(int iEl) void Patch::calcNormalEl2D(int iEl, bool scale)
{ {
const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace(); const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace();
...@@ -304,9 +304,9 @@ void Patch::calcScaledNormalEl2D(int iEl) ...@@ -304,9 +304,9 @@ void Patch::calcScaledNormalEl2D(int iEl)
// const bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization(); // const bool hasGeoNorm = ge && (ge->dim() == 2) && ge->haveParametrization();
for (int i=0; i<jac->getNumPrimMapNodes(); i++) { for (int i=0; i<jac->getNumPrimMapNodes(); i++) {
const int &iV = _el2V[iEl][i]; const int &iV = _el2V[iEl][i];
primNodesXYZ(i,0) = _xyz[iV].x(); primNodesXYZ(i, 0) = _xyz[iV].x();
primNodesXYZ(i,1) = _xyz[iV].y(); primNodesXYZ(i, 1) = _xyz[iV].y();
primNodesXYZ(i,2) = _xyz[iV].z(); primNodesXYZ(i, 2) = _xyz[iV].z();
// if (hasGeoNorm && (_vert[iV]->onWhat() == ge)) { // if (hasGeoNorm && (_vert[iV]->onWhat() == ge)) {
// double u, v; // double u, v;
// _vert[iV]->getParameter(0,u); // _vert[iV]->getParameter(0,u);
...@@ -319,16 +319,20 @@ void Patch::calcScaledNormalEl2D(int iEl) ...@@ -319,16 +319,20 @@ void Patch::calcScaledNormalEl2D(int iEl)
// geoNorm = ((GFace*)ge)->normal(param); // geoNorm = ((GFace*)ge)->normal(param);
// } // }
fullMatrix<double> &elNorm = _scaledNormEl[iEl]; fullMatrix<double> uElNorm(1, 3);
elNorm.resize(1,3); fullMatrix<double> &elNorm = scale ? _scaledNormEl[iEl] : uElNorm;
const double norm = jac->getPrimNormal2D(primNodesXYZ,elNorm); elNorm.resize(1, 3);
double factor = 1./norm; const double norm = jac->getPrimNormal2D(primNodesXYZ, elNorm);
// if (hasGeoNorm) { if (scale) {
// const double scal = geoNorm(0)*elNorm(0,0)+geoNorm(1)*elNorm(0,1)+geoNorm(2)*elNorm(0,2); double factor = 1./norm;
// if (scal < 0.) factor = -factor; // if (hasGeoNorm) {
// } // const double scal = geoNorm(0)*elNorm(0,0)+geoNorm(1)*elNorm(0,1)+geoNorm(2)*elNorm(0,2);
elNorm.scale(factor); // Re-scaling normal here is faster than an extra scaling operation on the Jacobian // if (scal < 0.) factor = -factor;
// }
elNorm.scale(factor); // Re-scaling normal here is faster than an extra scaling operation on the Jacobian
}
else
_unitNormEl[iEl] = SVector3(uElNorm(0, 0), uElNorm(0, 1), uElNorm(0, 2));
} }
...@@ -510,6 +514,14 @@ void Patch::initNCJ() ...@@ -510,6 +514,14 @@ void Patch::initNCJ()
case TYPE_QUA: _nNCJEl[iEl] = qmQuadrangle::numNCJVal(); break; case TYPE_QUA: _nNCJEl[iEl] = qmQuadrangle::numNCJVal(); break;
} }
} }
// Set normals to 2D elements (with magnitude of inverse Jacobian) or initial
// Jacobians of 3D elements
if ((_dim == 2) && _unitNormEl.empty()) {
_unitNormEl.resize(nEl());
// for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
for (int iEl = 0; iEl < nEl(); iEl++) calcNormalEl2D(iEl, false);
}
} }
...@@ -519,27 +531,19 @@ void Patch::NCJ(int iEl, std::vector<double> &NCJ) ...@@ -519,27 +531,19 @@ void Patch::NCJ(int iEl, std::vector<double> &NCJ)
const int &numNodes = _nNodEl[iEl]; const int &numNodes = _nNodEl[iEl];
std::vector<int> &iVEl = _el2V[iEl]; std::vector<int> &iVEl = _el2V[iEl];
fullVector<double> val(numNCJVal);
switch(_el[iEl]->getType()) { // TODO: Complete with other types? switch(_el[iEl]->getType()) { // TODO: Complete with other types?
case TYPE_TRI: { case TYPE_TRI: {
const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1], &iV2 = _el2V[iEl][2]; const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1], &iV2 = _el2V[iEl][2];
qmTriangle::NCJ(_xyz[iV0].x(), _xyz[iV0].y(), _xyz[iV0].z(), qmTriangle::NCJ(_xyz[iV0], _xyz[iV1], _xyz[iV2], _unitNormEl[iEl], NCJ);
_xyz[iV1].x(), _xyz[iV1].y(), _xyz[iV1].z(),
_xyz[iV2].x(), _xyz[iV2].y(), _xyz[iV2].z(), val);
break; break;
} }
case TYPE_QUA: { case TYPE_QUA: {
const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1], const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1],
&iV2 = _el2V[iEl][2], &iV3 = _el2V[iEl][3]; &iV2 = _el2V[iEl][2], &iV3 = _el2V[iEl][3];
qmQuadrangle::NCJ(_xyz[iV0].x(), _xyz[iV0].y(), _xyz[iV0].z(), qmQuadrangle::NCJ(_xyz[iV0], _xyz[iV1], _xyz[iV2], _xyz[iV3], _unitNormEl[iEl], NCJ);
_xyz[iV1].x(), _xyz[iV1].y(), _xyz[iV1].z(),
_xyz[iV2].x(), _xyz[iV2].y(), _xyz[iV2].z(),
_xyz[iV3].x(), _xyz[iV3].y(), _xyz[iV3].z(), val);
break; break;
} }
} }
for (int l = 0; l < numNCJVal; l++) NCJ[l] = val(l);
} }
...@@ -549,29 +553,23 @@ void Patch::NCJAndGradients(int iEl, std::vector<double> &NCJ, std::vector<doubl ...@@ -549,29 +553,23 @@ void Patch::NCJAndGradients(int iEl, std::vector<double> &NCJ, std::vector<doubl
const int &numNodes = _nNodEl[iEl]; const int &numNodes = _nNodEl[iEl];
std::vector<int> &iVEl = _el2V[iEl]; std::vector<int> &iVEl = _el2V[iEl];
fullMatrix<double> gradVal(numNCJVal, 3*numNodes+1); std::vector<double> dNCJ(numNCJVal*3*numNodes);
switch(_el[iEl]->getType()) { // TODO: Complete with other types? switch(_el[iEl]->getType()) { // TODO: Complete with other types?
case TYPE_TRI: { case TYPE_TRI: {
const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1], &iV2 = _el2V[iEl][2]; const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1], &iV2 = _el2V[iEl][2];
qmTriangle::NCJAndGradients(_xyz[iV0].x(), _xyz[iV0].y(), _xyz[iV0].z(), qmTriangle::NCJAndGradients(_xyz[iV0], _xyz[iV1], _xyz[iV2],
_xyz[iV1].x(), _xyz[iV1].y(), _xyz[iV1].z(), _unitNormEl[iEl], NCJ, dNCJ);
_xyz[iV2].x(), _xyz[iV2].y(), _xyz[iV2].z(), gradVal);
break; break;
} }
case TYPE_QUA: { case TYPE_QUA: {
const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1], const int &iV0 = _el2V[iEl][0], &iV1 = _el2V[iEl][1],
&iV2 = _el2V[iEl][2], &iV3 = _el2V[iEl][3]; &iV2 = _el2V[iEl][2], &iV3 = _el2V[iEl][3];
qmQuadrangle::NCJAndGradients(_xyz[iV0].x(), _xyz[iV0].y(), _xyz[iV0].z(), qmQuadrangle::NCJAndGradients(_xyz[iV0], _xyz[iV1], _xyz[iV2], _xyz[iV3],
_xyz[iV1].x(), _xyz[iV1].y(), _xyz[iV1].z(), _unitNormEl[iEl], NCJ, dNCJ);
_xyz[iV2].x(), _xyz[iV2].y(), _xyz[iV2].z(),
_xyz[iV3].x(), _xyz[iV3].y(), _xyz[iV3].z(), gradVal);
break; break;
} }
} }
// NCJ
for (int l = 0; l < numNCJVal; l++) NCJ[l] = gradVal(l,3*numNodes);
// Gradients of the NCJ // Gradients of the NCJ
int iPC = 0; int iPC = 0;
std::vector<SPoint3> gXyzV(numNCJVal); std::vector<SPoint3> gXyzV(numNCJVal);
...@@ -579,14 +577,15 @@ void Patch::NCJAndGradients(int iEl, std::vector<double> &NCJ, std::vector<doubl ...@@ -579,14 +577,15 @@ void Patch::NCJAndGradients(int iEl, std::vector<double> &NCJ, std::vector<doubl
for (int i = 0; i < numNodes; i++) { for (int i = 0; i < numNodes; i++) {
int &iFVi = _el2FV[iEl][i]; int &iFVi = _el2FV[iEl][i];
if (iFVi >= 0) { if (iFVi >= 0) {
for (int l = 0; l < numNCJVal; l++)
gXyzV [l] = SPoint3(gradVal(l,i+0*numNodes), gradVal(l,i+1*numNodes),
gradVal(l,i+2*numNodes));
_coordFV[iFVi]->gXyz2gUvw(_uvw[iFVi],gXyzV,gUvwV);
for (int l = 0; l < numNCJVal; l++) { for (int l = 0; l < numNCJVal; l++) {
gNCJ[indGNCJ(iEl,l,iPC)] = gUvwV[l][0]; int ind = (l*numNodes+i)*3;
if (_nPCFV[iFVi] >= 2) gNCJ[indGNCJ(iEl,l,iPC+1)] = gUvwV[l][1]; gXyzV[l] = SPoint3(dNCJ[ind], dNCJ[ind+1], dNCJ[ind+2]);
if (_nPCFV[iFVi] == 3) gNCJ[indGNCJ(iEl,l,iPC+2)] = gUvwV[l][2]; }
_coordFV[iFVi]->gXyz2gUvw(_uvw[iFVi], gXyzV, gUvwV);
for (int l = 0; l < numNCJVal; l++) {
gNCJ[indGNCJ(iEl, l, iPC)] = gUvwV[l][0];
if (_nPCFV[iFVi] >= 2) gNCJ[indGNCJ(iEl, l, iPC+1)] = gUvwV[l][1];
if (_nPCFV[iFVi] == 3) gNCJ[indGNCJ(iEl, l, iPC+2)] = gUvwV[l][2];
} }
iPC += _nPCFV[iFVi]; iPC += _nPCFV[iFVi];
} }
......
contrib/MeshOptimizer/MeshOptPatch.h
+ 2
− 1
View file @ 3c62eecb
...@@ -144,10 +144,11 @@ private: ...@@ -144,10 +144,11 @@ private:
std::vector<fullMatrix<double> > _scaledNormEl; // Normals to 2D elements for Jacobian regularization and scaling std::vector<fullMatrix<double> > _scaledNormEl; // Normals to 2D elements for Jacobian regularization and scaling
std::vector<double> _invStraightJac; // Initial Jacobians for 3D elements std::vector<double> _invStraightJac; // Initial Jacobians for 3D elements
// void calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl); // void calcScaledNormalEl2D(const std::map<MElement*,GEntity*> &element2entity, int iEl);
void calcScaledNormalEl2D(int iEl); void calcNormalEl2D(int iEl, bool scale);
// Mesh quality // Mesh quality
std::vector<int> _nNCJEl; // Number of NCJ values for an el. std::vector<int> _nNCJEl; // Number of NCJ values for an el.
std::vector<SVector3> _unitNormEl; // Normals to 2D elements for NCJ regularization
}; };
......
contrib/MeshQualityOptimizer/MeshQualityObjContribNCJ.h
+ 1
− 1
View file @ 3c62eecb
...@@ -86,7 +86,7 @@ void ObjContribNCJ<FuncType>::updateMinMax() ...@@ -86,7 +86,7 @@ void ObjContribNCJ<FuncType>::updateMinMax()
for (int iEl = 0; iEl < _mesh->nEl(); iEl++) { for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
std::vector<double> NCJ(_mesh->nNCJEl(iEl)); // Scaled Jacobians std::vector<double> NCJ(_mesh->nNCJEl(iEl)); // Scaled Jacobians
_mesh->NCJ(iEl, NCJ); _mesh->NCJ(iEl, NCJ);
for (int l = 0; l < _mesh->nBezEl(iEl); l++) { // Check each Bezier coeff. for (int l = 0; l < _mesh->nNCJEl(iEl); l++) { // Check each Bezier coeff.
_min = std::min(_min, NCJ[l]); _min = std::min(_min, NCJ[l]);
_max = std::max(_max, NCJ[l]); _max = std::max(_max, NCJ[l]);
} }
......
contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
+ 1
− 0
View file @ 3c62eecb
...@@ -84,6 +84,7 @@ void MeshQualityOptimizer(GModel *gm, MeshQualOptParameters &p) ...@@ -84,6 +84,7 @@ void MeshQualityOptimizer(GModel *gm, MeshQualOptParameters &p)
ObjContribScaledNodeDispSq<ObjContribFuncSimple> nodeDistFunc(p.weightFixed, p.weightFree); ObjContribScaledNodeDispSq<ObjContribFuncSimple> nodeDistFunc(p.weightFixed, p.weightFree);
ObjContribNCJ<ObjContribFuncBarrierMovMin> minNCJBarFunc(1.); ObjContribNCJ<ObjContribFuncBarrierMovMin> minNCJBarFunc(1.);
minNCJBarFunc.setTarget(p.minTargetNCJ, 1.); minNCJBarFunc.setTarget(p.minTargetNCJ, 1.);
// minNCJBarFunc.setTarget(p.minTargetNCJ, 0.866025404);
ObjContribNCJ<ObjContribFuncBarrierFixMinMovMax> minMaxNCJBarFunc(1.); ObjContribNCJ<ObjContribFuncBarrierFixMinMovMax> minMaxNCJBarFunc(1.);
minMaxNCJBarFunc.setTarget(p.maxTargetNCJ, 1.); minMaxNCJBarFunc.setTarget(p.maxTargetNCJ, 1.);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment