Skip to content
Snippets Groups Projects
Commit 9fd86cf9 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

pp

parent acb56d98
Branches
Tags
No related merge requests found
...@@ -122,11 +122,12 @@ void MElement::scaledJacRange(double &jmin, double &jmax, GEntity *ge) ...@@ -122,11 +122,12 @@ void MElement::scaledJacRange(double &jmin, double &jmax, GEntity *ge)
getNodesCoord(nodesXYZ); getNodesCoord(nodesXYZ);
fullVector<double> SJi(numJacNodes), Bi(numJacNodes); fullVector<double> SJi(numJacNodes), Bi(numJacNodes);
jac->getScaledJacobian(nodesXYZ,SJi); jac->getScaledJacobian(nodesXYZ,SJi);
if (ge && (ge->dim() == 2) && ge->haveParametrization()) { // If parametrized surface entity provided... if (ge && (ge->dim() == 2) && ge->haveParametrization()) {
SVector3 geoNorm(0.,0.,0.); // ... correct Jacobian sign with geometrical normal // If parametrized surface entity provided...
SVector3 geoNorm(0.,0.,0.);
// ... correct Jacobian sign with geometrical normal
for (int i=0; i<jac->getNumPrimMapNodes(); i++) { for (int i=0; i<jac->getNumPrimMapNodes(); i++) {
MVertex *vert = getVertex(i); MVertex *vert = getVertex(i);
GEntity *vge = vert->onWhat();
if (vert->onWhat() == ge) { if (vert->onWhat() == ge) {
double u, v; double u, v;
vert->getParameter(0,u); vert->getParameter(0,u);
...@@ -134,13 +135,15 @@ void MElement::scaledJacRange(double &jmin, double &jmax, GEntity *ge) ...@@ -134,13 +135,15 @@ void MElement::scaledJacRange(double &jmin, double &jmax, GEntity *ge)
geoNorm += ((GFace*)ge)->normal(SPoint2(u,v)); geoNorm += ((GFace*)ge)->normal(SPoint2(u,v));
} }
} }
if (geoNorm.normSq() == 0.) { // If no vertex on surface or average is zero, take normal at barycenter if (geoNorm.normSq() == 0.) {
// If no vertex on surface or average is zero, take normal at barycenter
SPoint2 param = ((GFace*)ge)->parFromPoint(barycenter(true),false); SPoint2 param = ((GFace*)ge)->parFromPoint(barycenter(true),false);
geoNorm = ((GFace*)ge)->normal(param); geoNorm = ((GFace*)ge)->normal(param);
} }
fullMatrix<double> elNorm(1,3); fullMatrix<double> elNorm(1,3);
jac->getPrimNormal2D(nodesXYZ,elNorm); jac->getPrimNormal2D(nodesXYZ,elNorm);
const double scal = geoNorm(0)*elNorm(0,0)+geoNorm(1)*elNorm(0,1)+geoNorm(2)*elNorm(0,2); const double scal = geoNorm(0) * elNorm(0,0) + geoNorm(1) * elNorm(0,1) +
geoNorm(2) * elNorm(0,2);
if (scal < 0.) SJi.scale(-1.); if (scal < 0.) SJi.scale(-1.);
} }
jac->lag2Bez(SJi,Bi); jac->lag2Bez(SJi,Bi);
......
...@@ -580,9 +580,9 @@ static void FindConnectedRegions(std::vector<GRegion*> &delaunay, ...@@ -580,9 +580,9 @@ static void FindConnectedRegions(std::vector<GRegion*> &delaunay,
if (!nbVolumes)return; if (!nbVolumes)return;
while (delaunay.size()){ while (delaunay.size()){
std::set<GRegion*> oneDomain; std::set<GRegion*> oneDomain;
std::stack<GRegion*> _stack; std::stack<GRegion*> _stack;
GRegion *r = delaunay[0]; GRegion *r = delaunay[0];
_stack.push(r); _stack.push(r);
while(!_stack.empty()){ while(!_stack.empty()){
r = _stack.top(); r = _stack.top();
_stack.pop(); _stack.pop();
...@@ -605,7 +605,7 @@ static void FindConnectedRegions(std::vector<GRegion*> &delaunay, ...@@ -605,7 +605,7 @@ static void FindConnectedRegions(std::vector<GRegion*> &delaunay,
delaunay=temp1; delaunay=temp1;
} }
Msg::Info("Delaunay Meshing %d volumes with %d connected components", Msg::Info("Delaunay Meshing %d volumes with %d connected components",
nbVolumes,connected.size()); nbVolumes,connected.size());
} }
static void Mesh3D(GModel *m) static void Mesh3D(GModel *m)
...@@ -635,7 +635,7 @@ static void Mesh3D(GModel *m) ...@@ -635,7 +635,7 @@ static void Mesh3D(GModel *m)
FindConnectedRegions(delaunay, connected); FindConnectedRegions(delaunay, connected);
// remove quads elements for volumes that are recombined // remove quads elements for volumes that are recombined
// pragma OMP ICI ?? // pragma OMP ICI ??
for(unsigned int i = 0; i < connected.size(); i++){ for(unsigned int i = 0; i < connected.size(); i++){
for(unsigned j=0;j<connected[i].size();j++){ for(unsigned j=0;j<connected[i].size();j++){
GRegion *gr = connected[i][j]; GRegion *gr = connected[i][j];
...@@ -647,7 +647,7 @@ static void Mesh3D(GModel *m) ...@@ -647,7 +647,7 @@ static void Mesh3D(GModel *m)
} }
} }
// pragma OMP ICI ?? // pragma OMP ICI ??
for(unsigned int i = 0; i < connected.size(); i++){ for(unsigned int i = 0; i < connected.size(); i++){
MeshDelaunayVolume(connected[i]); MeshDelaunayVolume(connected[i]);
...@@ -808,6 +808,7 @@ void GenerateMesh(GModel *m, int ask) ...@@ -808,6 +808,7 @@ void GenerateMesh(GModel *m, int ask)
p.BARRIER_MIN = CTX::instance()->mesh.hoThresholdMin; p.BARRIER_MIN = CTX::instance()->mesh.hoThresholdMin;
p.BARRIER_MAX = CTX::instance()->mesh.hoThresholdMax; p.BARRIER_MAX = CTX::instance()->mesh.hoThresholdMax;
p.dim = GModel::current()->getDim(); p.dim = GModel::current()->getDim();
//p.optPrimSurfMesh = true;
HighOrderMeshOptimizer(GModel::current(), p); HighOrderMeshOptimizer(GModel::current(), p);
} }
#else #else
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment