diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 9f956eb87c6488dcb2685e0137fd032726c5fa4d..922b312d3003eddb4ae99524a32fe46c71f5e77c 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -30,6 +30,8 @@ JacobianBasis::JacobianBasis(int tag)
const int jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
+ fullMatrix<double> lagPoints; // Sampling points
+
switch (parentType) {
case TYPE_PNT :
lagPoints = gmshGeneratePointsLine(0);
@@ -38,19 +40,19 @@ JacobianBasis::JacobianBasis(int tag)
lagPoints = gmshGeneratePointsLine(jacobianOrder);
break;
case TYPE_TRI :
- lagPoints = gmshGeneratePointsTriangle(jacobianOrder);
+ lagPoints = gmshGeneratePointsTriangle(jacobianOrder,false);
break;
case TYPE_QUA :
- lagPoints = gmshGeneratePointsQuadrangle(jacobianOrder);
+ lagPoints = gmshGeneratePointsQuadrangle(jacobianOrder,false);
break;
case TYPE_TET :
- lagPoints = gmshGeneratePointsTetrahedron(jacobianOrder);
+ lagPoints = gmshGeneratePointsTetrahedron(jacobianOrder,false);
break;
case TYPE_PRI :
- lagPoints = gmshGeneratePointsPrism(jacobianOrder);
+ lagPoints = gmshGeneratePointsPrism(jacobianOrder,false);
break;
case TYPE_HEX :
- lagPoints = gmshGeneratePointsHexahedron(jacobianOrder);
+ lagPoints = gmshGeneratePointsHexahedron(jacobianOrder,false);
break;
case TYPE_PYR :
lagPoints = generateJacPointsPyramid(jacobianOrder);
@@ -118,7 +120,35 @@ JacobianBasis::JacobianBasis(int tag)
primGradShapeBarycenterY(j) = barDPsi[j][1];
primGradShapeBarycenterZ(j) = barDPsi[j][2];
}
+
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)
@@ -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
-// straight element for regularization
-void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
+namespace {
+
+// 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 : {
fullMatrix<double> normals(2,3);
getPrimNormals1D(nodesXYZ,normals);
- getSignedJacobian(nodesXYZ,normals,jacobian);
+ computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normals,jacobian);
break;
}
case 2 : {
fullMatrix<double> normal(1,3);
getPrimNormal2D(nodesXYZ,normal);
- getSignedJacobian(nodesXYZ,normal,jacobian);
+ computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normal,jacobian);
break;
}
case 0 :
case 3 : {
fullMatrix<double> dum;
- getSignedJacobian(nodesXYZ,dum,jacobian);
+ computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,dum,jacobian);
break;
}
@@ -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
-// straight element for regularization and scaling
-void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
+// Calculate scaled (signed) Jacobian for one element, with normal vectors to straight element for regularization
+// and scaling. Evaluation points depend on the given matrices for shape function gradients.
+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 : {
fullMatrix<double> normals(2,3);
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
- getSignedJacobian(nodesXYZ,normals,jacobian);
+ computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normals,jacobian);
break;
}
@@ -248,7 +348,7 @@ void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVe
fullMatrix<double> normal(1,3);
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
- getSignedJacobian(nodesXYZ,normal,jacobian);
+ computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normal,jacobian);
break;
}
@@ -256,7 +356,7 @@ void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVe
case 3 : {
fullMatrix<double> dum;
const double scale = 1./getPrimJac3D(nodesXYZ);
- getSignedJacobian(nodesXYZ,dum,jacobian);
+ computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,dum,jacobian);
jacobian.scale(scale);
break;
}
@@ -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
-// of line and surface Jacobians in 3D
-void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals, fullVector<double> &jacobian) const
+// Calculate (signed) Jacobian for several elements.
+// Evaluation points depend on the given matrices for shape function gradients.
+// TODO: Optimize and test 1D & 2D
+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()) {
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;
}
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 int numEl = nodesX.size2();
+ fullMatrix<double> dxdX(nJacNodes,numEl), dydX(nJacNodes,numEl), dzdX(nJacNodes,numEl);
+ gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.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);
- 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;
}
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);
+ const int numEl = nodesX.size2();
+ fullMatrix<double> dxdX(nJacNodes,numEl), dydX(nJacNodes,numEl), dzdX(nJacNodes,numEl);
+ fullMatrix<double> dxdY(nJacNodes,numEl), dydY(nJacNodes,numEl), dzdY(nJacNodes,numEl);
+ gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
+ gSMatY.mult(nodesX, dxdY); gSMatY.mult(nodesY, dydY); gSMatY.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 < 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;
}
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);
- }
+ const int numEl = nodesX.size2();
+ fullMatrix<double> dxdX(nJacNodes,numEl), dydX(nJacNodes,numEl), dzdX(nJacNodes,numEl);
+ fullMatrix<double> dxdY(nJacNodes,numEl), dydY(nJacNodes,numEl), dzdY(nJacNodes,numEl);
+ fullMatrix<double> dxdZ(nJacNodes,numEl), dydZ(nJacNodes,numEl), dzdZ(nJacNodes,numEl);
+ gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
+ gSMatY.mult(nodesX, dxdY); gSMatY.mult(nodesY, dydY); gSMatY.mult(nodesZ, dzdY);
+ gSMatZ.mult(nodesX, dxdZ); gSMatZ.mult(nodesY, dydZ); gSMatZ.mult(nodesZ, dzdZ);
+ for (int iEl = 0; iEl < numEl; iEl++)
+ 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;
}
@@ -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
-// of line and surface Jacobians in 3D
-void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals,
- fullMatrix<double> &JDJ) const
+// Calculate (signed) Jacobian and its gradients for one element, with normal vectors to straight element
+// for regularization. Evaluation points depend on the given matrices for shape function gradients.
+void JacobianBasis::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
{
switch (bezier->getDim()) {
case 0 : {
- for (int i = 0; i < numJacNodes; i++) {
+ for (int i = 0; i < nJacNodes; i++) {
for (int j = 0; j < numMapNodes; j++) {
JDJ (i,j) = 0.;
JDJ (i,j+1*numMapNodes) = 0.;
@@ -343,14 +475,14 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
}
case 1 : {
- fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3);
- gradShapeMatX.mult(nodesXYZ, dxyzdX);
- for (int i = 0; i < numJacNodes; i++) {
+ fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3);
+ gSMatX.mult(nodesXYZ, dxyzdX);
+ for (int i = 0; i < nJacNodes; 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);
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+1*numMapNodes) = dPhidX * dzdY * dxdZ - dPhidX * dxdY * dzdZ;
JDJ (i,j+2*numMapNodes) = dPhidX * dxdY * dydZ - dPhidX * dydY * dxdZ;
@@ -361,16 +493,16 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
}
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++) {
+ fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3);
+ gSMatX.mult(nodesXYZ, dxyzdX);
+ gSMatY.mult(nodesXYZ, dxyzdY);
+ for (int i = 0; i < nJacNodes; 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);
for (int j = 0; j < numMapNodes; j++) {
- const double &dPhidX = gradShapeMatX(i,j);
- const double &dPhidY = gradShapeMatY(i,j);
+ const double &dPhidX = gSMatX(i,j);
+ const double &dPhidY = gSMatY(i,j);
JDJ (i,j) =
dPhidX * dydY * dzdZ + dzdX * dPhidY * dydZ +
dPhidX * dzdY * dydZ - dydX * dPhidY * dzdZ;
@@ -389,18 +521,18 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
}
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++) {
+ fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3), dxyzdZ(nJacNodes,3);
+ gSMatX.mult(nodesXYZ, dxyzdX);
+ gSMatY.mult(nodesXYZ, dxyzdY);
+ gSMatZ.mult(nodesXYZ, dxyzdZ);
+ for (int i = 0; i < nJacNodes; 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);
for (int j = 0; j < numMapNodes; j++) {
- const double &dPhidX = gradShapeMatX(i,j);
- const double &dPhidY = gradShapeMatY(i,j);
- const double &dPhidZ = gradShapeMatZ(i,j);
+ const double &dPhidX = gSMatX(i,j);
+ const double &dPhidY = gSMatY(i,j);
+ const double &dPhidZ = gSMatZ(i,j);
JDJ (i,j) =
dPhidX * dydY * dzdZ + dzdX * dPhidY * dydZ +
dydX * dzdY * dPhidZ - dzdX * dydY * dPhidZ -
@@ -425,7 +557,7 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
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)
@@ -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)
{
int nbMonomials = (order+3)*((order+3)+1)*(2*(order+3)+1)/6 - 5;
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 31890fc5f0d69c636b706aea668b744ffb0d6f9a..e28f949958da0d1c787c447da6088be186f17645 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -16,19 +16,36 @@ class JacobianBasis {
private:
const bezierBasis *bezier;
- fullMatrix<double> lagPoints; // sampling point
fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
+ fullMatrix<double> gradShapeMatXFast, gradShapeMatYFast, gradShapeMatZFast;
fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY, primGradShapeBarycenterZ;
fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
int numJacNodes, numPrimJacNodes;
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 :
JacobianBasis(int tag);
- // get methods
+ // Get methods
inline int getNumJacNodes() const { return numJacNodes; }
+ inline int getNumJacNodesFast() const { return numJacNodesFast; }
inline int getNumMapNodes() const { return numMapNodes; }
inline int getNumPrimJacNodes() const { return numPrimJacNodes; }
inline int getNumPrimMapNodes() const { return numPrimMapNodes; }
@@ -36,22 +53,44 @@ class JacobianBasis {
inline int getNumSubNodes() const { return bezier->subDivisor.size1(); }
inline int getNumLagCoeff() const { return bezier->getNumLagCoeff(); }
- //
+ // Jacobian evaluation methods
double getPrimNormals1D(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;
- void getSignedJacobian(const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals, fullVector<double> &jacobian) const;
- void getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const;
+ inline void getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
+ getSignedJacAndGradientsGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,normals,JDJ);
+ }
+ 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,
const fullMatrix<double> &nodesXYZStraight,
fullVector<double> &lambdaJ , fullMatrix<double> &gradLambdaJ) const;
- void getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
- void getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
- void getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
-
+ inline void getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
+ getSignedJacobianGeneral(numJacNodes,gradShapeMatX,gradShapeMatY,gradShapeMatZ,nodesXYZ,jacobian);
+ }
+ 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 {
bezier->matrixLag2Bez.mult(jac,bez);
@@ -66,7 +105,7 @@ class JacobianBasis {
bezier->subDivisor.mult(bez,result);
}
- //
+ // Jacobian basis order and pyramidal basis
static int jacobianOrder(int parentType, int order);
static fullMatrix<double> generateJacMonomialsPyramid(int order);
static fullMatrix<double> generateJacPointsPyramid(int order);
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.cpp b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
index 64b2cc200b18ebf14299f20e473e1c782887502d..45db5d249c0fd6edf679a25c0a13ff8665a3d92f 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.cpp
@@ -68,8 +68,8 @@ static inline double compute_f1(double v, double barrier)
}
OptHOM::OptHOM(const std::set<MElement*> &els, std::set<MVertex*> & toFix,
- bool fixBndNodes) :
- mesh(els, toFix, fixBndNodes)
+ bool fixBndNodes, bool fastJacEval) :
+ mesh(els, toFix, fixBndNodes, fastJacEval)
{
_optimizeMetricMin = false;
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHOM.h b/contrib/HighOrderMeshOptimizer/OptHOM.h
index 19e5c42931a0645873c370afdc61f511f38cec45..debbc6a856c7858f4ad78d3419d136fdc834ae0c 100644
--- a/contrib/HighOrderMeshOptimizer/OptHOM.h
+++ b/contrib/HighOrderMeshOptimizer/OptHOM.h
@@ -44,7 +44,7 @@ class OptHOM
public:
Mesh mesh;
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
// 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
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 85efcef199e98bcb9afb543e877c4a17500cea15..4ed1a1a119ad8f8780b280324a31b6d6d62dbb55 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -35,7 +35,8 @@
#include "ParamCoord.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();
@@ -65,7 +66,7 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn
it != els.end(); ++it, ++iEl) {
_el[iEl] = *it;
const JacobianBasis *jac = _el[iEl]->getJacobianFuncSpace();
- _nBezEl[iEl] = jac->getNumJacNodes();
+ _nBezEl[iEl] = _fastJacEval ? jac->getNumJacNodesFast() : jac->getNumJacNodes();
_nNodEl[iEl] = jac->getNumMapNodes();
for (int iVEl = 0; iVEl < jac->getNumMapNodes(); iVEl++) {
MVertex *vert = _el[iEl]->getVertex(iVEl);
@@ -73,10 +74,12 @@ Mesh::Mesh(const std::set<MElement*> &els, std::set<MVertex*> &toFix, bool fixBn
_el2V[iEl].push_back(iV);
const int nPCV = _pc->nCoord(vert);
bool isFV = false;
- if (fixBndNodes)
- isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
- else
- isFV = (vert->onWhat()->dim() >= 1) && (toFix.find(vert) == toFix.end());
+ if (nPCV > 0) {
+ if (fixBndNodes)
+ isFV = (vert->onWhat()->dim() == _dim) && (toFix.find(vert) == toFix.end());
+ else
+ isFV = (vert->onWhat()->dim() >= 1) && (toFix.find(vert) == toFix.end());
+ }
if (isFV) {
int iFV = addFreeVert(vert,iV,nPCV,toFix);
_el2FV[iEl].push_back(iFV);
@@ -285,10 +288,9 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ,
std::vector<double> &gSJ)
{
const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace();
- const int &numJacNodes = jacBasis->getNumJacNodes();
- const int &numMapNodes = jacBasis->getNumMapNodes();
- fullMatrix<double> JDJ (numJacNodes,3*numMapNodes+1);
- fullMatrix<double> BDB (numJacNodes,3*numMapNodes+1);
+ const int &numJacNodes = _nBezEl[iEl];
+ const int &numMapNodes = _nNodEl[iEl];
+ fullMatrix<double> JDJ(numJacNodes,3*numMapNodes+1), BDB(numJacNodes,3*numMapNodes+1);
// Coordinates of nodes
fullMatrix<double> nodesXYZ(numMapNodes,3), normals(_dim,3);
@@ -301,11 +303,17 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ,
// Calculate Jacobian and gradients, scale if 3D (already scaled by
// 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]);
// Transform Jacobian and gradients from Lagrangian to Bezier basis
- jacBasis->lag2Bez(JDJ,BDB);
+ if (_fastJacEval)
+ BDB = JDJ;
+ else
+ jacBasis->lag2Bez(JDJ,BDB);
// Scaled jacobian
for (int l = 0; l < numJacNodes; l++) sJ [l] = BDB (l,3*numMapNodes);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index ec49b7717918ae9de808c777621c78bf4b6ec1e9..06876da687614a01987d7249b2bc954e7b2a57f2 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -41,7 +41,7 @@
class Mesh
{
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 int nVert() { return _vert.size(); }
@@ -75,6 +75,8 @@ private:
int _dim;
// Total nb. of parametric coordinates
int _nPC;
+ // Use fast Jacobian estimation?
+ bool _fastJacEval;
// List of elements
std::vector<MElement*> _el;
// Normals to 2D elements for Jacobian regularization and scaling
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 5a988b76b987d5fe35721011384fe72c5702edb6..b9ca0aaae0ce7a2173f8404a10943fe1b8097e13 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -308,8 +308,12 @@ static void optimizeConnectedBlobs
//std::ostringstream ossI1;
//ossI1 << "initial_ITER_" << i << ".msh";
//temp.mesh.writeMSH(ossI1.str().c_str());
- int success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX,
- false, samples, p.itMax, p.optPassMax);
+ int success = -1;
+ 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) {
Msg::Info("Jacobian optimization succeed, starting svd optimization");
success = temp.optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC, p.BARRIER_MAX,
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.h b/contrib/HighOrderMeshOptimizer/ParamCoord.h
index e812d78d8da34c808de6ca3d854898a6479f6bb0..9b50e50280a609c0dd36b8cb078074538d32276f 100644
--- a/contrib/HighOrderMeshOptimizer/ParamCoord.h
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.h
@@ -92,7 +92,7 @@ public:
class ParamCoordParent : public ParamCoord
{
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 SPoint3 getUvw(MVertex* vert);
virtual SPoint3 uvw2Xyz(MVertex* vert, const SPoint3 &uvw);