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

pp

parent a3198234
No related branches found
No related tags found
No related merge requests found
...@@ -255,8 +255,8 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve, ...@@ -255,8 +255,8 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
US[u0 < u1 ? j + 1 : nPts - 1 - (j + 1)]); US[u0 < u1 ? j + 1 : nPts - 1 - (j + 1)]);
if (displ2D || displ3D){ if (displ2D || displ3D){
SPoint3 pc2 = edge.interpolate(t); SPoint3 pc2 = edge.interpolate(t);
if (displ2D)displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z())); if(displ2D) displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
if (displ3D)displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z())); if(displ3D) displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
} }
} }
temp.push_back(v); temp.push_back(v);
...@@ -317,8 +317,8 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve, ...@@ -317,8 +317,8 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, US[j + 1], VS[j + 1]); v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, US[j + 1], VS[j + 1]);
if (displ2D || displ3D){ if (displ2D || displ3D){
SPoint3 pc2 = edge.interpolate(t); SPoint3 pc2 = edge.interpolate(t);
if (displ3D) displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z())); if(displ3D) displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
if (displ2D) displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z())); if(displ2D) displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
} }
} }
temp.push_back(v); temp.push_back(v);
...@@ -373,6 +373,7 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele, ...@@ -373,6 +373,7 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
if(gf->geomType() == GEntity::DiscreteSurface || if(gf->geomType() == GEntity::DiscreteSurface ||
gf->geomType() == GEntity::BoundaryLayerSurface) gf->geomType() == GEntity::BoundaryLayerSurface)
linear = true; linear = true;
for(int i = 0; i < ele->getNumFaces(); i++){ for(int i = 0; i < ele->getNumFaces(); i++){
MFace face = ele->getFace(i); MFace face = ele->getFace(i);
faceContainer::iterator fIter = faceVertices.find(face); faceContainer::iterator fIter = faceVertices.find(face);
...@@ -387,8 +388,8 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele, ...@@ -387,8 +388,8 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
for(int k = 0; k < incomplete->getNumVertices(); k++) for(int k = 0; k < incomplete->getNumVertices(); k++)
reparamOK &= reparamMeshVertexOnFace(incomplete->getVertex(k), gf, pts[k]); reparamOK &= reparamMeshVertexOnFace(incomplete->getVertex(k), gf, pts[k]);
} }
int start = face.getNumVertices()*(nPts+1); int start = face.getNumVertices() * (nPts + 1);
const fullMatrix<double> &points = ele->getFunctionSpace(nPts+1)->points; const fullMatrix<double> &points = ele->getFunctionSpace(nPts + 1)->points;
for(int k = start; k < points.size1(); k++){ for(int k = start; k < points.size1(); k++){
MVertex *v; MVertex *v;
const double t1 = points(k, 0); const double t1 = points(k, 0);
...@@ -425,8 +426,8 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele, ...@@ -425,8 +426,8 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
} }
if(displ3D || displ2D){ if(displ3D || displ2D){
SPoint3 pc2 = face.interpolate(t1, t2); SPoint3 pc2 = face.interpolate(t1, t2);
if(displ3D)displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z())); if(displ3D) displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
if(displ2D)displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z())); if(displ2D) displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
} }
} }
// should be expensive -> induces a new search each time // should be expensive -> induces a new search each time
...@@ -485,21 +486,21 @@ static void reorientQuadPoints(std::vector<MVertex*> &vtcs, int orientation, ...@@ -485,21 +486,21 @@ static void reorientQuadPoints(std::vector<MVertex*> &vtcs, int orientation,
else{ else{
int i1,i2,i3,i4; int i1,i2,i3,i4;
if (!swap){ if (!swap){
if (orientation == 0){i1 = 0; i2=1; i3=2; i4 = 3;} if (orientation == 0){ i1 = 0; i2 = 1; i3 = 2; i4 = 3; }
else if (orientation == 1){i1 = 3; i2=0; i3=1; i4 = 2;} else if (orientation == 1){ i1 = 3; i2 = 0; i3 = 1; i4 = 2; }
else if (orientation == 2){i1 = 2; i2=3; i3=0; i4 = 1;} else if (orientation == 2){ i1 = 2; i2 = 3; i3 = 0; i4 = 1; }
else if (orientation == 3){i1 = 1; i2=2; i3=3; i4 = 0;} else if (orientation == 3){ i1 = 1; i2 = 2; i3 = 3; i4 = 0; }
} }
else{ else{
if (orientation == 0){i1 = 0; i2=3; i3=2; i4 = 1;} if (orientation == 0){ i1 = 0; i2 = 3; i3 = 2; i4 = 1; }
else if (orientation == 3){i1 = 3; i2=2; i3=1; i4 = 0;} else if (orientation == 3){ i1 = 3; i2 = 2; i3 = 1; i4 = 0; }
else if (orientation == 2){i1 = 2; i2=1; i3=0; i4 = 3;} else if (orientation == 2){ i1 = 2; i2 = 1; i3 = 0; i4 = 3; }
else if (orientation == 1){i1 = 1; i2=0; i3=3; i4 = 2;} else if (orientation == 1){ i1 = 1; i2 = 0; i3 = 3; i4 = 2; }
} }
int indices[4] = {i1,i2,i3,i4}; int indices[4] = {i1, i2, i3, i4};
for (int i=0;i<4;i++)tmp[i] = vtcs[start+indices[i]]; for (int i = 0; i < 4; i++) tmp[i] = vtcs[start + indices[i]];
for (int i=0;i<4;i++)vtcs[start+i] = tmp[i]; for (int i = 0; i < 4; i++) vtcs[start + i] = tmp[i];
// EDGES // EDGES
for (int iEdge=0;iEdge<4;iEdge++){ for (int iEdge=0;iEdge<4;iEdge++){
...@@ -530,10 +531,10 @@ static void reorientQuadPoints(std::vector<MVertex*> &vtcs, int orientation, ...@@ -530,10 +531,10 @@ static void reorientQuadPoints(std::vector<MVertex*> &vtcs, int orientation,
else if (p1 == 0 && p2 == 3){ else if (p1 == 0 && p2 == 3){
for (int i = 4+4*nbP-1; i >= 4+3*nbP; i--) tmp[index++] = vtcs[i]; for (int i = 4+4*nbP-1; i >= 4+3*nbP; i--) tmp[index++] = vtcs[i];
} }
else Msg::Error("ouuls"); else Msg::Error("Something wrong in reorientQuadPoints");
} }
for (int i=0;i<index;i++)vtcs[start+4+i] = tmp[i]; for (int i = 0; i < index; i++)vtcs[start+4+i] = tmp[i];
start += (4+index); start += (4 + index);
} }
order -= 2; order -= 2;
...@@ -651,7 +652,6 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v ...@@ -651,7 +652,6 @@ static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &v
} }
else if(face.getNumVertices() == 4){ // quad face else if(face.getNumVertices() == 4){ // quad face
for (int k = start; k < points.size1(); k++) { for (int k = start; k < points.size1(); k++) {
// parameters are between -1 and 1
double t1 = points(k, 0); double t1 = points(k, 0);
double t2 = points(k, 1); double t2 = points(k, 1);
SPoint3 pc = face.interpolate(t1, t2); SPoint3 pc = face.interpolate(t1, t2);
...@@ -867,7 +867,7 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices, ...@@ -867,7 +867,7 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
t->getVertex(3), ve, nPts + 1); t->getVertex(3), ve, nPts + 1);
getRegionVertices(gr, &incpl, t, vr, linear, nPts); getRegionVertices(gr, &incpl, t, vr, linear, nPts);
ve.insert(ve.end(), vr.begin(), vr.end()); ve.insert(ve.end(), vr.begin(), vr.end());
MTetrahedron* n = new MTetrahedronN(t->getVertex(0), t->getVertex(1), MTetrahedron *n = new MTetrahedronN(t->getVertex(0), t->getVertex(1),
t->getVertex(2), t->getVertex(3), ve, nPts + 1); t->getVertex(2), t->getVertex(3), ve, nPts + 1);
tetrahedra2.push_back(n); tetrahedra2.push_back(n);
} }
...@@ -910,7 +910,7 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices, ...@@ -910,7 +910,7 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
h->getVertex(6), h->getVertex(7), ve, nPts + 1); h->getVertex(6), h->getVertex(7), ve, nPts + 1);
getRegionVertices(gr, &incpl, h, vr, linear, nPts); getRegionVertices(gr, &incpl, h, vr, linear, nPts);
ve.insert(ve.end(), vr.begin(), vr.end()); ve.insert(ve.end(), vr.begin(), vr.end());
MHexahedron* n = new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2), MHexahedron *n = new MHexahedronN(h->getVertex(0), h->getVertex(1), h->getVertex(2),
h->getVertex(3), h->getVertex(4), h->getVertex(5), h->getVertex(3), h->getVertex(4), h->getVertex(5),
h->getVertex(6), h->getVertex(7), ve, nPts + 1); h->getVertex(6), h->getVertex(7), ve, nPts + 1);
hexahedra2.push_back(n); hexahedra2.push_back(n);
...@@ -1024,7 +1024,7 @@ void SetOrder1(GModel *m) ...@@ -1024,7 +1024,7 @@ void SetOrder1(GModel *m)
} }
void checkHighOrderTriangles(const char* cc, GModel *m, void checkHighOrderTriangles(const char* cc, GModel *m,
std::vector<MElement*> &bad, double &minJGlob) std::vector<MElement*> &bad, double &minJGlob)
{ {
bad.clear(); bad.clear();
minJGlob = 1.0; minJGlob = 1.0;
...@@ -1055,10 +1055,11 @@ void checkHighOrderTriangles(const char* cc, GModel *m, ...@@ -1055,10 +1055,11 @@ void checkHighOrderTriangles(const char* cc, GModel *m,
else if (disto < 0.2) nbfair++; else if (disto < 0.2) nbfair++;
} }
} }
if(!count) return;
if (minJGlob > 0) if (minJGlob > 0)
Msg::Info("%s : Worst Face Distorsion Mapping %g Gamma %g Nb elem. (0<d<0.2) = %d", Msg::Info("%s : Worst Face Distorsion Mapping %g Gamma %g Nb elem. (0<d<0.2) = %d",
cc, minJGlob, minGGlob,nbfair ); cc, minJGlob, minGGlob,nbfair );
else else
Msg::Warning("%s : Worst Face Distorsion Mapping %g (%d negative jacobians) " Msg::Warning("%s : Worst Face Distorsion Mapping %g (%d negative jacobians) "
"Worst Gamma %g Avg Smoothness %g", cc, minJGlob, bad.size(), "Worst Gamma %g Avg Smoothness %g", cc, minJGlob, bad.size(),
minGGlob, avg / (count ? count : 1)); minGGlob, avg / (count ? count : 1));
...@@ -1085,6 +1086,7 @@ static void checkHighOrderTetrahedron(const char* cc, GModel *m, ...@@ -1085,6 +1086,7 @@ static void checkHighOrderTetrahedron(const char* cc, GModel *m,
else if (disto < 0.2) nbfair++; else if (disto < 0.2) nbfair++;
} }
} }
if(!count) return;
if (minJGlob < 0) if (minJGlob < 0)
Msg::Warning("%s : Worst Tetrahedron Smoothness %g Gamma %g NbFair = %d NbBad = %d", Msg::Warning("%s : Worst Tetrahedron Smoothness %g Gamma %g NbFair = %d NbBad = %d",
cc, minJGlob, minGGlob, nbfair, bad.size()); cc, minJGlob, minGGlob, nbfair, bad.size());
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment