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

Added fast Jacobian estimation and robustness improvements in OptHOM

parent 136a0425
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 289 additions and 186 deletions
Numeric/JacobianBasis.cpp
+ 204
− 154
View file @ 2356d75a
...@@ -30,6 +30,8 @@ JacobianBasis::JacobianBasis(int tag) ...@@ -30,6 +30,8 @@ JacobianBasis::JacobianBasis(int tag)
const int jacobianOrder = JacobianBasis::jacobianOrder(parentType, order); const int jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1); const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
fullMatrix<double> lagPoints; // Sampling points
switch (parentType) { switch (parentType) {
case TYPE_PNT : case TYPE_PNT :
lagPoints = gmshGeneratePointsLine(0); lagPoints = gmshGeneratePointsLine(0);
...@@ -38,19 +40,19 @@ JacobianBasis::JacobianBasis(int tag) ...@@ -38,19 +40,19 @@ JacobianBasis::JacobianBasis(int tag)
lagPoints = gmshGeneratePointsLine(jacobianOrder); lagPoints = gmshGeneratePointsLine(jacobianOrder);
break; break;
case TYPE_TRI : case TYPE_TRI :
lagPoints = gmshGeneratePointsTriangle(jacobianOrder); lagPoints = gmshGeneratePointsTriangle(jacobianOrder,false);
break; break;
case TYPE_QUA : case TYPE_QUA :
lagPoints = gmshGeneratePointsQuadrangle(jacobianOrder); lagPoints = gmshGeneratePointsQuadrangle(jacobianOrder,false);
break; break;
case TYPE_TET : case TYPE_TET :
lagPoints = gmshGeneratePointsTetrahedron(jacobianOrder); lagPoints = gmshGeneratePointsTetrahedron(jacobianOrder,false);
break; break;
case TYPE_PRI : case TYPE_PRI :
lagPoints = gmshGeneratePointsPrism(jacobianOrder); lagPoints = gmshGeneratePointsPrism(jacobianOrder,false);
break; break;
case TYPE_HEX : case TYPE_HEX :
lagPoints = gmshGeneratePointsHexahedron(jacobianOrder); lagPoints = gmshGeneratePointsHexahedron(jacobianOrder,false);
break; break;
case TYPE_PYR : case TYPE_PYR :
lagPoints = generateJacPointsPyramid(jacobianOrder); lagPoints = generateJacPointsPyramid(jacobianOrder);
...@@ -118,7 +120,35 @@ JacobianBasis::JacobianBasis(int tag) ...@@ -118,7 +120,35 @@ JacobianBasis::JacobianBasis(int tag)
primGradShapeBarycenterY(j) = barDPsi[j][1]; primGradShapeBarycenterY(j) = barDPsi[j][1];
primGradShapeBarycenterZ(j) = barDPsi[j][2]; primGradShapeBarycenterZ(j) = barDPsi[j][2];
} }
delete[] barDPsi; delete[] barDPsi;
// Compute "fast" Jacobian evaluation matrices (at 1st order nodes + barycenter)
// numJacNodesFast = numPrimMapNodes+1;
// fullMatrix<double> lagPointsFast(numJacNodesFast,3); // Sampling points
// lagPointsFast.copy(primMapBasis->points,0,numPrimMapNodes,0,3,0,0); // 1st order nodes
// lagPointsFast(numPrimMapNodes,0) = barycenter[0]; // Last point = barycenter
// lagPointsFast(numPrimMapNodes,1) = barycenter[1];
// lagPointsFast(numPrimMapNodes,2) = barycenter[2];
fullMatrix<double> lagPointsFast(numJacNodesFast,3); // Sampling points
lagPointsFast(0,0) = barycenter[0]; // Last point = barycenter
lagPointsFast(0,1) = barycenter[1];
lagPointsFast(0,2) = barycenter[2];
fullMatrix<double> allDPsiFast;
mapBasis->df(lagPointsFast, allDPsiFast);
gradShapeMatXFast.resize(numJacNodesFast, numMapNodes);
gradShapeMatYFast.resize(numJacNodesFast, numMapNodes);
gradShapeMatZFast.resize(numJacNodesFast, numMapNodes);
for (int i=0; i<numJacNodesFast; i++) {
for (int j=0; j<numMapNodes; j++) {
gradShapeMatXFast(i, j) = allDPsiFast(j, 3*i);
gradShapeMatYFast(i, j) = allDPsiFast(j, 3*i+1);
gradShapeMatZFast(i, j) = allDPsiFast(j, 3*i+2);
}
}
} }
// Computes (unit) normals to straight line element at barycenter (with norm of gradient as return value) // Computes (unit) normals to straight line element at barycenter (with norm of gradient as return value)
...@@ -197,31 +227,97 @@ double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ) const ...@@ -197,31 +227,97 @@ double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ) const
} }
// Calculate (signed) Jacobian at mapping's nodes for one element, with normal vectors to namespace {
// straight element for regularization
void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const // Calculate (signed) Jacobian for one element, given the gradients of the shape functions
// and vectors for regularization of line and surface Jacobians in 3D.
inline void computeSignedJacobian(int dim, int numJacNodes, const fullMatrix<double> &gradShapeMatX,
const fullMatrix<double> &gradShapeMatY, const fullMatrix<double> &gradShapeMatZ,
const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
fullVector<double> &jacobian)
{ {
switch (bezier->getDim()) { switch (dim) {
case 0 : {
for (int i = 0; i < numJacNodes; i++) jacobian(i) = 1.;
break;
}
case 1 : {
fullMatrix<double> dxyzdX(numJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX);
for (int i = 0; i < numJacNodes; i++) {
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2);
const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2);
jacobian(i) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ);
}
break;
}
case 2 : {
fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX);
gradShapeMatY.mult(nodesXYZ, dxyzdY);
for (int i = 0; i < numJacNodes; i++) {
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
const double &dxdZ = normals(0,0), &dydZ = normals(0,1), &dzdZ = normals(0,2);
jacobian(i) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ);
}
break;
}
case 3 : {
fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3), dxyzdZ(numJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX);
gradShapeMatY.mult(nodesXYZ, dxyzdY);
gradShapeMatZ.mult(nodesXYZ, dxyzdZ);
for (int i = 0; i < numJacNodes; i++) {
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2);
jacobian(i) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ);
}
break;
}
}
}
} // namespace
// Calculate (signed) Jacobian for one element, with normal vectors to straight element for regularization.
// Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
const int dim = bezier->getDim();
switch (dim) {
case 1 : { case 1 : {
fullMatrix<double> normals(2,3); fullMatrix<double> normals(2,3);
getPrimNormals1D(nodesXYZ,normals); getPrimNormals1D(nodesXYZ,normals);
getSignedJacobian(nodesXYZ,normals,jacobian); computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normals,jacobian);
break; break;
} }
case 2 : { case 2 : {
fullMatrix<double> normal(1,3); fullMatrix<double> normal(1,3);
getPrimNormal2D(nodesXYZ,normal); getPrimNormal2D(nodesXYZ,normal);
getSignedJacobian(nodesXYZ,normal,jacobian); computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normal,jacobian);
break; break;
} }
case 0 : case 0 :
case 3 : { case 3 : {
fullMatrix<double> dum; fullMatrix<double> dum;
getSignedJacobian(nodesXYZ,dum,jacobian); computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,dum,jacobian);
break; break;
} }
...@@ -229,18 +325,22 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVe ...@@ -229,18 +325,22 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVe
} }
// Calculate scaled (signed) Jacobian at mapping's nodes for one element, with normal vectors to // Calculate scaled (signed) Jacobian for one element, with normal vectors to straight element for regularization
// straight element for regularization and scaling // and scaling. Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const void JacobianBasis::getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{ {
switch (bezier->getDim()) { const int dim = bezier->getDim();
switch (dim) {
case 1 : { case 1 : {
fullMatrix<double> normals(2,3); fullMatrix<double> normals(2,3);
const double scale = 1./getPrimNormals1D(nodesXYZ,normals); const double scale = 1./getPrimNormals1D(nodesXYZ,normals);
normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale; // Faster to scale 1 normal than afterwards normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale; // Faster to scale 1 normal than afterwards
getSignedJacobian(nodesXYZ,normals,jacobian); computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normals,jacobian);
break; break;
} }
...@@ -248,7 +348,7 @@ void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVe ...@@ -248,7 +348,7 @@ void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVe
fullMatrix<double> normal(1,3); fullMatrix<double> normal(1,3);
const double scale = 1./getPrimNormal2D(nodesXYZ,normal); const double scale = 1./getPrimNormal2D(nodesXYZ,normal);
normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
getSignedJacobian(nodesXYZ,normal,jacobian); computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normal,jacobian);
break; break;
} }
...@@ -256,7 +356,7 @@ void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVe ...@@ -256,7 +356,7 @@ void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVe
case 3 : { case 3 : {
fullMatrix<double> dum; fullMatrix<double> dum;
const double scale = 1./getPrimJac3D(nodesXYZ); const double scale = 1./getPrimJac3D(nodesXYZ);
getSignedJacobian(nodesXYZ,dum,jacobian); computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,dum,jacobian);
jacobian.scale(scale); jacobian.scale(scale);
break; break;
} }
...@@ -265,55 +365,84 @@ void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVe ...@@ -265,55 +365,84 @@ void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVe
} }
// Calculate (signed) Jacobian at mapping's nodes for one element, given vectors for regularization // Calculate (signed) Jacobian for several elements.
// of line and surface Jacobians in 3D // Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, // TODO: Optimize and test 1D & 2D
const fullMatrix<double> &normals, fullVector<double> &jacobian) const void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
{ {
switch (bezier->getDim()) { switch (bezier->getDim()) {
case 0 : { case 0 : {
for (int i = 0; i < numJacNodes; i++) jacobian(i) = 1.; const int numEl = nodesX.size2();
for (int iEl = 0; iEl < numEl; iEl++)
for (int i = 0; i < nJacNodes; i++) jacobian(i,iEl) = 1.;
break; break;
} }
case 1 : { case 1 : {
fullMatrix<double> dxyzdX(numJacNodes,3); const int numEl = nodesX.size2();
gradShapeMatX.mult(nodesXYZ, dxyzdX); fullMatrix<double> dxdX(nJacNodes,numEl), dydX(nJacNodes,numEl), dzdX(nJacNodes,numEl);
for (int i = 0; i < numJacNodes; i++) { gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2); for (int iEl = 0; iEl < numEl; iEl++) {
fullMatrix<double> nodesXYZ(numPrimJacNodes,3);
for (int i = 0; i < numPrimJacNodes; i++) {
nodesXYZ(i,0) = nodesX(i,iEl);
nodesXYZ(i,1) = nodesY(i,iEl);
nodesXYZ(i,2) = nodesZ(i,iEl);
}
fullMatrix<double> normals(2,3);
getPrimNormals1D(nodesXYZ,normals);
const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2); const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2);
const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2); const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2);
jacobian(i) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ); for (int i = 0; i < nJacNodes; i++)
jacobian(i,iEl) = calcDet3D(dxdX(i,iEl),dydX(i,iEl),dzdX(i,iEl),
dxdY,dydY,dzdY,
dxdZ,dydZ,dzdZ);
} }
break; break;
} }
case 2 : { case 2 : {
fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3); const int numEl = nodesX.size2();
gradShapeMatX.mult(nodesXYZ, dxyzdX); fullMatrix<double> dxdX(nJacNodes,numEl), dydX(nJacNodes,numEl), dzdX(nJacNodes,numEl);
gradShapeMatY.mult(nodesXYZ, dxyzdY); fullMatrix<double> dxdY(nJacNodes,numEl), dydY(nJacNodes,numEl), dzdY(nJacNodes,numEl);
for (int i = 0; i < numJacNodes; i++) { gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2); gSMatY.mult(nodesX, dxdY); gSMatY.mult(nodesY, dydY); gSMatY.mult(nodesZ, dzdY);
const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2); for (int iEl = 0; iEl < numEl; iEl++) {
const double &dxdZ = normals(0,0), &dydZ = normals(0,1), &dzdZ = normals(0,2); fullMatrix<double> nodesXYZ(numPrimJacNodes,3);
jacobian(i) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ); for (int i = 0; i < numPrimJacNodes; i++) {
nodesXYZ(i,0) = nodesX(i,iEl);
nodesXYZ(i,1) = nodesY(i,iEl);
nodesXYZ(i,2) = nodesZ(i,iEl);
}
fullMatrix<double> normal(1,3);
getPrimNormal2D(nodesXYZ,normal);
const double &dxdZ = normal(0,0), &dydZ = normal(0,1), &dzdZ = normal(0,2);
for (int i = 0; i < nJacNodes; i++)
jacobian(i,iEl) = calcDet3D(dxdX(i,iEl),dydX(i,iEl),dzdX(i,iEl),
dxdY(i,iEl),dydY(i,iEl),dzdY(i,iEl),
dxdZ,dydZ,dzdZ);
} }
break; break;
} }
case 3 : { case 3 : {
fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3), dxyzdZ(numJacNodes,3); const int numEl = nodesX.size2();
gradShapeMatX.mult(nodesXYZ, dxyzdX); fullMatrix<double> dxdX(nJacNodes,numEl), dydX(nJacNodes,numEl), dzdX(nJacNodes,numEl);
gradShapeMatY.mult(nodesXYZ, dxyzdY); fullMatrix<double> dxdY(nJacNodes,numEl), dydY(nJacNodes,numEl), dzdY(nJacNodes,numEl);
gradShapeMatZ.mult(nodesXYZ, dxyzdZ); fullMatrix<double> dxdZ(nJacNodes,numEl), dydZ(nJacNodes,numEl), dzdZ(nJacNodes,numEl);
for (int i = 0; i < numJacNodes; i++) { gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2); gSMatY.mult(nodesX, dxdY); gSMatY.mult(nodesY, dydY); gSMatY.mult(nodesZ, dzdY);
const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2); gSMatZ.mult(nodesX, dxdZ); gSMatZ.mult(nodesY, dydZ); gSMatZ.mult(nodesZ, dzdZ);
const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2); for (int iEl = 0; iEl < numEl; iEl++)
jacobian(i) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ); for (int i = 0; i < nJacNodes; i++)
} jacobian(i,iEl) = calcDet3D(dxdX(i,iEl),dydX(i,iEl),dzdX(i,iEl),
dxdY(i,iEl),dydY(i,iEl),dzdY(i,iEl),
dxdZ(i,iEl),dydZ(i,iEl),dzdZ(i,iEl));
break; break;
} }
...@@ -321,17 +450,20 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, ...@@ -321,17 +450,20 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ,
} }
// Calculate (signed) Jacobian at mapping's nodes for one element, given vectors for regularization // Calculate (signed) Jacobian and its gradients for one element, with normal vectors to straight element
// of line and surface Jacobians in 3D // for regularization. Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ, void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &normals, const fullMatrix<double> &gSMatY,
fullMatrix<double> &JDJ) const const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals,
fullMatrix<double> &JDJ) const
{ {
switch (bezier->getDim()) { switch (bezier->getDim()) {
case 0 : { case 0 : {
for (int i = 0; i < numJacNodes; i++) { for (int i = 0; i < nJacNodes; i++) {
for (int j = 0; j < numMapNodes; j++) { for (int j = 0; j < numMapNodes; j++) {
JDJ (i,j) = 0.; JDJ (i,j) = 0.;
JDJ (i,j+1*numMapNodes) = 0.; JDJ (i,j+1*numMapNodes) = 0.;
...@@ -343,14 +475,14 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ, ...@@ -343,14 +475,14 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
} }
case 1 : { case 1 : {
fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3); fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX); gSMatX.mult(nodesXYZ, dxyzdX);
for (int i = 0; i < numJacNodes; i++) { for (int i = 0; i < nJacNodes; i++) {
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2); const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2); const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2);
const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2); const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2);
for (int j = 0; j < numMapNodes; j++) { for (int j = 0; j < numMapNodes; j++) {
const double &dPhidX = gradShapeMatX(i,j); const double &dPhidX = gSMatX(i,j);
JDJ (i,j) = dPhidX * dydY * dzdZ + dPhidX * dzdY * dydZ; JDJ (i,j) = dPhidX * dydY * dzdZ + dPhidX * dzdY * dydZ;
JDJ (i,j+1*numMapNodes) = dPhidX * dzdY * dxdZ - dPhidX * dxdY * dzdZ; JDJ (i,j+1*numMapNodes) = dPhidX * dzdY * dxdZ - dPhidX * dxdY * dzdZ;
JDJ (i,j+2*numMapNodes) = dPhidX * dxdY * dydZ - dPhidX * dydY * dxdZ; JDJ (i,j+2*numMapNodes) = dPhidX * dxdY * dydZ - dPhidX * dydY * dxdZ;
...@@ -361,16 +493,16 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ, ...@@ -361,16 +493,16 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
} }
case 2 : { case 2 : {
fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3); fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX); gSMatX.mult(nodesXYZ, dxyzdX);
gradShapeMatY.mult(nodesXYZ, dxyzdY); gSMatY.mult(nodesXYZ, dxyzdY);
for (int i = 0; i < numJacNodes; i++) { for (int i = 0; i < nJacNodes; i++) {
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2); const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2); const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
const double &dxdZ = normals(0,0), &dydZ = normals(0,1), &dzdZ = normals(0,2); const double &dxdZ = normals(0,0), &dydZ = normals(0,1), &dzdZ = normals(0,2);
for (int j = 0; j < numMapNodes; j++) { for (int j = 0; j < numMapNodes; j++) {
const double &dPhidX = gradShapeMatX(i,j); const double &dPhidX = gSMatX(i,j);
const double &dPhidY = gradShapeMatY(i,j); const double &dPhidY = gSMatY(i,j);
JDJ (i,j) = JDJ (i,j) =
dPhidX * dydY * dzdZ + dzdX * dPhidY * dydZ + dPhidX * dydY * dzdZ + dzdX * dPhidY * dydZ +
dPhidX * dzdY * dydZ - dydX * dPhidY * dzdZ; dPhidX * dzdY * dydZ - dydX * dPhidY * dzdZ;
...@@ -389,18 +521,18 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ, ...@@ -389,18 +521,18 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
} }
case 3 : { case 3 : {
fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3), dxyzdZ(numJacNodes,3); fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3), dxyzdZ(nJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX); gSMatX.mult(nodesXYZ, dxyzdX);
gradShapeMatY.mult(nodesXYZ, dxyzdY); gSMatY.mult(nodesXYZ, dxyzdY);
gradShapeMatZ.mult(nodesXYZ, dxyzdZ); gSMatZ.mult(nodesXYZ, dxyzdZ);
for (int i = 0; i < numJacNodes; i++) { for (int i = 0; i < nJacNodes; i++) {
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2); const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2); const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2); const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2);
for (int j = 0; j < numMapNodes; j++) { for (int j = 0; j < numMapNodes; j++) {
const double &dPhidX = gradShapeMatX(i,j); const double &dPhidX = gSMatX(i,j);
const double &dPhidY = gradShapeMatY(i,j); const double &dPhidY = gSMatY(i,j);
const double &dPhidZ = gradShapeMatZ(i,j); const double &dPhidZ = gSMatZ(i,j);
JDJ (i,j) = JDJ (i,j) =
dPhidX * dydY * dzdZ + dzdX * dPhidY * dydZ + dPhidX * dydY * dzdZ + dzdX * dPhidY * dydZ +
dydX * dzdY * dPhidZ - dzdX * dydY * dPhidZ - dydX * dzdY * dPhidZ - dzdX * dydY * dPhidZ -
...@@ -425,7 +557,7 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ, ...@@ -425,7 +557,7 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ, void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &nodesXYZStraight, const fullMatrix<double> &nodesXYZStraight,
fullVector<double> &lambdaJ , fullMatrix<double> &gradLambdaJ) const fullVector<double> &lambdaJ, fullMatrix<double> &gradLambdaJ) const
{ {
// jacobian of the straight elements (only triangles for now) // jacobian of the straight elements (only triangles for now)
...@@ -476,88 +608,6 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ, ...@@ -476,88 +608,6 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
} }
// Calculate (signed) Jacobian at mapping's nodes for several elements
// TODO: Optimize and test 1D & 2D
void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
{
switch (bezier->getDim()) {
case 0 : {
const int numEl = nodesX.size2();
for (int iEl = 0; iEl < numEl; iEl++)
for (int i = 0; i < numJacNodes; i++) jacobian(i,iEl) = 1.;
break;
}
case 1 : {
const int numEl = nodesX.size2();
fullMatrix<double> dxdX(numJacNodes,numEl), dydX(numJacNodes,numEl), dzdX(numJacNodes,numEl);
gradShapeMatX.mult(nodesX, dxdX); gradShapeMatX.mult(nodesY, dydX); gradShapeMatX.mult(nodesZ, dzdX);
for (int iEl = 0; iEl < numEl; iEl++) {
fullMatrix<double> nodesXYZ(numPrimJacNodes,3);
for (int i = 0; i < numPrimJacNodes; i++) {
nodesXYZ(i,0) = nodesX(i,iEl);
nodesXYZ(i,1) = nodesY(i,iEl);
nodesXYZ(i,2) = nodesZ(i,iEl);
}
fullMatrix<double> normals(2,3);
getPrimNormals1D(nodesXYZ,normals);
const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2);
const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2);
for (int i = 0; i < numJacNodes; i++)
jacobian(i,iEl) = calcDet3D(dxdX(i,iEl),dydX(i,iEl),dzdX(i,iEl),
dxdY,dydY,dzdY,
dxdZ,dydZ,dzdZ);
}
break;
}
case 2 : {
const int numEl = nodesX.size2();
fullMatrix<double> dxdX(numJacNodes,numEl), dydX(numJacNodes,numEl), dzdX(numJacNodes,numEl);
fullMatrix<double> dxdY(numJacNodes,numEl), dydY(numJacNodes,numEl), dzdY(numJacNodes,numEl);
gradShapeMatX.mult(nodesX, dxdX); gradShapeMatX.mult(nodesY, dydX); gradShapeMatX.mult(nodesZ, dzdX);
gradShapeMatY.mult(nodesX, dxdY); gradShapeMatY.mult(nodesY, dydY); gradShapeMatY.mult(nodesZ, dzdY);
for (int iEl = 0; iEl < numEl; iEl++) {
fullMatrix<double> nodesXYZ(numPrimJacNodes,3);
for (int i = 0; i < numPrimJacNodes; i++) {
nodesXYZ(i,0) = nodesX(i,iEl);
nodesXYZ(i,1) = nodesY(i,iEl);
nodesXYZ(i,2) = nodesZ(i,iEl);
}
fullMatrix<double> normal(1,3);
getPrimNormal2D(nodesXYZ,normal);
const double &dxdZ = normal(0,0), &dydZ = normal(0,1), &dzdZ = normal(0,2);
for (int i = 0; i < numJacNodes; i++)
jacobian(i,iEl) = calcDet3D(dxdX(i,iEl),dydX(i,iEl),dzdX(i,iEl),
dxdY(i,iEl),dydY(i,iEl),dzdY(i,iEl),
dxdZ,dydZ,dzdZ);
}
break;
}
case 3 : {
const int numEl = nodesX.size2();
fullMatrix<double> dxdX(numJacNodes,numEl), dydX(numJacNodes,numEl), dzdX(numJacNodes,numEl);
fullMatrix<double> dxdY(numJacNodes,numEl), dydY(numJacNodes,numEl), dzdY(numJacNodes,numEl);
fullMatrix<double> dxdZ(numJacNodes,numEl), dydZ(numJacNodes,numEl), dzdZ(numJacNodes,numEl);
gradShapeMatX.mult(nodesX, dxdX); gradShapeMatX.mult(nodesY, dydX); gradShapeMatX.mult(nodesZ, dzdX);
gradShapeMatY.mult(nodesX, dxdY); gradShapeMatY.mult(nodesY, dydY); gradShapeMatY.mult(nodesZ, dzdY);
gradShapeMatZ.mult(nodesX, dxdZ); gradShapeMatZ.mult(nodesY, dydZ); gradShapeMatZ.mult(nodesZ, dzdZ);
for (int iEl = 0; iEl < numEl; iEl++)
for (int i = 0; i < numJacNodes; i++)
jacobian(i,iEl) = calcDet3D(dxdX(i,iEl),dydX(i,iEl),dzdX(i,iEl),
dxdY(i,iEl),dydY(i,iEl),dzdY(i,iEl),
dxdZ(i,iEl),dydZ(i,iEl),dzdZ(i,iEl));
break;
}
}
}
fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order) fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
{ {
int nbMonomials = (order+3)*((order+3)+1)*(2*(order+3)+1)/6 - 5; int nbMonomials = (order+3)*((order+3)+1)*(2*(order+3)+1)/6 - 5;
......
Numeric/JacobianBasis.h
+ 52
− 13
View file @ 2356d75a
...@@ -16,19 +16,36 @@ class JacobianBasis { ...@@ -16,19 +16,36 @@ class JacobianBasis {
private: private:
const bezierBasis *bezier; const bezierBasis *bezier;
fullMatrix<double> lagPoints; // sampling point
fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ; fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
fullMatrix<double> gradShapeMatXFast, gradShapeMatYFast, gradShapeMatZFast;
fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY, primGradShapeBarycenterZ; fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY, primGradShapeBarycenterZ;
fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac. fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
int numJacNodes, numPrimJacNodes; int numJacNodes, numPrimJacNodes;
int numMapNodes, numPrimMapNodes; int numMapNodes, numPrimMapNodes;
static const int numJacNodesFast = 1;
void getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
void getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
void getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
void getSignedJacAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
fullMatrix<double> &JDJ) const;
public : public :
JacobianBasis(int tag); JacobianBasis(int tag);
// get methods // Get methods
inline int getNumJacNodes() const { return numJacNodes; } inline int getNumJacNodes() const { return numJacNodes; }
inline int getNumJacNodesFast() const { return numJacNodesFast; }
inline int getNumMapNodes() const { return numMapNodes; } inline int getNumMapNodes() const { return numMapNodes; }
inline int getNumPrimJacNodes() const { return numPrimJacNodes; } inline int getNumPrimJacNodes() const { return numPrimJacNodes; }
inline int getNumPrimMapNodes() const { return numPrimMapNodes; } inline int getNumPrimMapNodes() const { return numPrimMapNodes; }
...@@ -36,22 +53,44 @@ class JacobianBasis { ...@@ -36,22 +53,44 @@ class JacobianBasis {
inline int getNumSubNodes() const { return bezier->subDivisor.size1(); } inline int getNumSubNodes() const { return bezier->subDivisor.size1(); }
inline int getNumLagCoeff() const { return bezier->getNumLagCoeff(); } inline int getNumLagCoeff() const { return bezier->getNumLagCoeff(); }
// // Jacobian evaluation methods
double getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const; double getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
double getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const; double getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
double getPrimJac3D(const fullMatrix<double> &nodesXYZ) const; double getPrimJac3D(const fullMatrix<double> &nodesXYZ) const;
void getSignedJacobian(const fullMatrix<double> &nodesXYZ, inline void getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals, fullVector<double> &jacobian) const; const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
void getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ, getSignedJacAndGradientsGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,normals,JDJ);
const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const; }
inline void getSignedJacAndGradientsFast(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
getSignedJacAndGradientsGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
gradShapeMatZFast,nodesXYZ,normals,JDJ);
}
void getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ, void getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &nodesXYZStraight, const fullMatrix<double> &nodesXYZStraight,
fullVector<double> &lambdaJ , fullMatrix<double> &gradLambdaJ) const; fullVector<double> &lambdaJ , fullMatrix<double> &gradLambdaJ) const;
void getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const; inline void getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
void getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY, getSignedJacobianGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,jacobian);
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const; }
void getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const; inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
}
inline void getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
getSignedJacobianGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,
gradShapeMatZ,nodesX,nodesY,nodesZ,jacobian);
}
inline void getSignedJacobianFast(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
gradShapeMatZFast,nodesX,nodesY,nodesZ,jacobian);
}
inline void getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
getScaledJacobianGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,jacobian);
}
inline void getScaledJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
getScaledJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
}
// //
inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const { inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
bezier->matrixLag2Bez.mult(jac,bez); bezier->matrixLag2Bez.mult(jac,bez);
...@@ -66,7 +105,7 @@ class JacobianBasis { ...@@ -66,7 +105,7 @@ class JacobianBasis {
bezier->subDivisor.mult(bez,result); bezier->subDivisor.mult(bez,result);
} }
// // Jacobian basis order and pyramidal basis
static int jacobianOrder(int parentType, int order); static int jacobianOrder(int parentType, int order);
static fullMatrix<double> generateJacMonomialsPyramid(int order); static fullMatrix<double> generateJacMonomialsPyramid(int order);
static fullMatrix<double> generateJacPointsPyramid(int order); static fullMatrix<double> generateJacPointsPyramid(int order);
......
contrib/HighOrderMeshOptimizer/OptHOM.cpp
+ 2
− 2
View file @ 2356d75a
...@@ -68,8 +68,8 @@ static inline double compute_f1(double v, double barrier) ...@@ -68,8 +68,8 @@ static inline double compute_f1(double v, double barrier)
} }
OptHOM::OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix, OptHOM::OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix,
bool fixBndNodes) : bool fixBndNodes, bool fastJacEval) :
mesh(els, toFix, fixBndNodes) mesh(els, toFix, fixBndNodes, fastJacEval)
{ {
_optimizeMetricMin = false; _optimizeMetricMin = false;
} }
......
contrib/HighOrderMeshOptimizer/OptHOM.h
+ 1
− 1
View file @ 2356d75a
...@@ -44,7 +44,7 @@ class OptHOM ...@@ -44,7 +44,7 @@ class OptHOM
public: public:
Mesh mesh; Mesh mesh;
OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix, OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix,
bool fixBndNodes); bool fixBndNodes, bool fastJacEval = false);
// returns 1 if the mesh has been optimized with success i.e. all jacobians // returns 1 if the mesh has been optimized with success i.e. all jacobians
// are in the range; returns 0 if the mesh is valid (all jacobians positive, // are in the range; returns 0 if the mesh is valid (all jacobians positive,
// JMIN > 0) but JMIN < barrier_min || JMAX > barrier_max; returns -1 if the // JMIN > 0) but JMIN < barrier_min || JMAX > barrier_max; returns -1 if the
......
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+ 20
− 12
View file @ 2356d75a
...@@ -35,7 +35,8 @@ ...@@ -35,7 +35,8 @@
#include "ParamCoord.h" #include "ParamCoord.h"
#include "OptHomMesh.h" #include "OptHomMesh.h"
Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBndNodes) Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBndNodes, bool fastJacEval) :
_fastJacEval(fastJacEval)
{ {
_dim = (*els.begin())->getDim(); _dim = (*els.begin())->getDim();
...@@ -65,7 +66,7 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn ...@@ -65,7 +66,7 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn
it != els.end(); ++it, ++iEl) { it != els.end(); ++it, ++iEl) {
_el[iEl] = *it; _el[iEl] = *it;
const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace(); const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace();
_nBezEl[iEl] = jac->getNumJacNodes(); _nBezEl[iEl] = _fastJacEval ? jac->getNumJacNodesFast() : jac->getNumJacNodes();
_nNodEl[iEl] = jac->getNumMapNodes(); _nNodEl[iEl] = jac->getNumMapNodes();
for (int iVEl = 0; iVEl < jac->getNumMapNodes(); iVEl++) { for (int iVEl = 0; iVEl < jac->getNumMapNodes(); iVEl++) {
MVertex *vert = _el[iEl]->getVertex(iVEl); MVertex *vert = _el[iEl]->getVertex(iVEl);
...@@ -73,10 +74,12 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn ...@@ -73,10 +74,12 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn
_el2V[iEl].push_back(iV); _el2V[iEl].push_back(iV);
const int nPCV = _pc->nCoord(vert); const int nPCV = _pc->nCoord(vert);
bool isFV = false; bool isFV = false;
if (fixBndNodes) if (nPCV > 0) {
isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end()); if (fixBndNodes)
else isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
isFV = (vert->onWhat()->dim() >= 1) && (toFix.find(vert) == toFix.end()); else
isFV = (vert->onWhat()->dim() >= 1) && (toFix.find(vert) == toFix.end());
}
if (isFV) { if (isFV) {
int iFV = addFreeVert(vert,iV,nPCV,toFix); int iFV = addFreeVert(vert,iV,nPCV,toFix);
_el2FV[iEl].push_back(iFV); _el2FV[iEl].push_back(iFV);
...@@ -285,10 +288,9 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ, ...@@ -285,10 +288,9 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ,
std::vector<double> &gSJ) std::vector<double> &gSJ)
{ {
const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace(); const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace();
const int &numJacNodes = jacBasis->getNumJacNodes(); const int &numJacNodes = _nBezEl[iEl];
const int &numMapNodes = jacBasis->getNumMapNodes(); const int &numMapNodes = _nNodEl[iEl];
fullMatrix<double> JDJ (numJacNodes,3*numMapNodes+1); fullMatrix<double> JDJ(numJacNodes,3*numMapNodes+1), BDB(numJacNodes,3*numMapNodes+1);
fullMatrix<double> BDB (numJacNodes,3*numMapNodes+1);
// Coordinates of nodes // Coordinates of nodes
fullMatrix<double> nodesXYZ(numMapNodes,3), normals(_dim,3); fullMatrix<double> nodesXYZ(numMapNodes,3), normals(_dim,3);
...@@ -301,11 +303,17 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ, ...@@ -301,11 +303,17 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ,
// Calculate Jacobian and gradients, scale if 3D (already scaled by // Calculate Jacobian and gradients, scale if 3D (already scaled by
// regularization normals in 2D) // regularization normals in 2D)
jacBasis->getSignedJacAndGradients(nodesXYZ,_scaledNormEl[iEl],JDJ); if (_fastJacEval)
jacBasis->getSignedJacAndGradientsFast(nodesXYZ,_scaledNormEl[iEl],JDJ);
else
jacBasis->getSignedJacAndGradients(nodesXYZ,_scaledNormEl[iEl],JDJ);
if (_dim == 3) JDJ.scale(_invStraightJac[iEl]); if (_dim == 3) JDJ.scale(_invStraightJac[iEl]);
// Transform Jacobian and gradients from Lagrangian to Bezier basis // Transform Jacobian and gradients from Lagrangian to Bezier basis
jacBasis->lag2Bez(JDJ,BDB); if (_fastJacEval)
BDB = JDJ;
else
jacBasis->lag2Bez(JDJ,BDB);
// Scaled jacobian // Scaled jacobian
for (int l = 0; l < numJacNodes; l++) sJ [l] = BDB (l,3*numMapNodes); for (int l = 0; l < numJacNodes; l++) sJ [l] = BDB (l,3*numMapNodes);
......
contrib/HighOrderMeshOptimizer/OptHomMesh.h
+ 3
− 1
View file @ 2356d75a
...@@ -41,7 +41,7 @@ ...@@ -41,7 +41,7 @@
class Mesh class Mesh
{ {
public: public:
Mesh(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes); Mesh(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes, bool fastJacEval);
inline const int &nPC() { return _nPC; } inline const int &nPC() { return _nPC; }
inline int nVert() { return _vert.size(); } inline int nVert() { return _vert.size(); }
...@@ -75,6 +75,8 @@ private: ...@@ -75,6 +75,8 @@ private:
int _dim; int _dim;
// Total nb. of parametric coordinates // Total nb. of parametric coordinates
int _nPC; int _nPC;
// Use fast Jacobian estimation?
bool _fastJacEval;
// List of elements // List of elements
std::vector<MElement*> _el; std::vector<MElement*> _el;
// Normals to 2D elements for Jacobian regularization and scaling // Normals to 2D elements for Jacobian regularization and scaling
......
contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+ 6
− 2
View file @ 2356d75a
...@@ -308,8 +308,12 @@ static void optimizeConnectedBlobs ...@@ -308,8 +308,12 @@ static void optimizeConnectedBlobs
//std::ostringstream ossI1; //std::ostringstream ossI1;
//ossI1 << "initial_ITER_" << i << ".msh"; //ossI1 << "initial_ITER_" << i << ".msh";
//temp.mesh.writeMSH(ossI1.str().c_str()); //temp.mesh.writeMSH(ossI1.str().c_str());
int success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX, int success = -1;
false, samples, p.itMax, p.optPassMax); if (temp.mesh.nPC() == 0)
Msg::Info("Blob %i has no degree of freedom, skipping", i+1);
else
success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
if (success >= 0 && p.BARRIER_MIN_METRIC > 0) { if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
Msg::Info("Jacobian optimization succeed, starting svd optimization"); Msg::Info("Jacobian optimization succeed, starting svd optimization");
success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX, success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX,
......
contrib/HighOrderMeshOptimizer/ParamCoord.h
+ 1
− 1
View file @ 2356d75a
...@@ -92,7 +92,7 @@ public: ...@@ -92,7 +92,7 @@ public:
class ParamCoordParent : public ParamCoord class ParamCoordParent : public ParamCoord
{ {
public: public:
int nCoord(MVertex* vert) { return vert->onWhat()->dim(); } int nCoord(MVertex* vert) { return vert->onWhat()->haveParametrization() ? vert->onWhat()->dim() : 0; }
virtual void exportParamCoord(MVertex *v, const SPoint3 &uvw); virtual void exportParamCoord(MVertex *v, const SPoint3 &uvw);
virtual SPoint3 getUvw(MVertex* vert); virtual SPoint3 getUvw(MVertex* vert);
virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw); virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment