From 0f497647b651188c8fa4aceadc75d775120abc60 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Sun, 8 Dec 2013 23:42:14 +0000
Subject: [PATCH] Enforce volume positivity of 3D elements right after mesh
generation/optimization, instead of enforcing when writing to file
---
Geo/GModel.cpp | 13 +++++++++++++
Geo/GModel.h | 3 +++
Geo/MElement.cpp | 20 --------------------
3 files changed, 16 insertions(+), 20 deletions(-)
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index d2e04c5aeb..e988c6d8e4 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -538,6 +538,19 @@ int GModel::mesh(int dimension)
#endif
}
+bool GModel::setAllVolumesPositive()
+{
+ std::vector<GEntity*> entities;
+ getEntities(entities);
+ bool ok = true;
+ for(std::set<GRegion*,GEntityLessThan>::iterator itReg=regions.begin();
+ itReg != regions.end(); ++itReg) {
+ int nbEl = (*itReg)->getNumMeshElements();
+ for (int iEl=0; iEl<nbEl; ++iEl) ok = ok && (*itReg)->getMeshElement(iEl)->setVolumePositive();
+ }
+ return ok;
+}
+
int GModel::adaptMesh(std::vector<int> technique,
std::vector<simpleFunction<double>* > f,
std::vector<std::vector<double> > parameters,
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 8fc5068a07..d1aa2ec37a 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -455,6 +455,9 @@ class GModel
std::vector<std::vector<double> > parameters,
int niter, bool meshAll=false);
+ // Ensure that the Jacobian of all volume elements is positive
+ bool setAllVolumesPositive();
+
// make the mesh a high order mesh at order N
// linear is 1 if the high order points are not placed on the geometry of the model
// incomplete is 1 if incomplete basis are used
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 6457df851b..fa487a5a97 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -737,9 +737,6 @@ void MElement::writeMSH(FILE *fp, bool binary, int entity,
int type = getTypeForMSH();
if(!type) return;
- // if necessary, change the ordering of the vertices to get positive volume
- setVolumePositive();
-
std::vector<int> verts;
getVerticesIdForMSH(verts);
@@ -785,9 +782,6 @@ void MElement::writeMSH2(FILE *fp, double version, bool binary, int num,
if(!type) return;
- // if necessary, change the ordering of the vertices to get positive
- // volume
- setVolumePositive();
int n = getNumVerticesForMSH();
int par = (parentNum) ? 1 : 0;
int dom = (dom1Num) ? 2 : 0;
@@ -887,7 +881,6 @@ void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
const char *str = getStringForPOS();
if(!str) return;
- setVolumePositive();
int n = getNumVertices();
fprintf(fp, "%s(", str);
for(int i = 0; i < n; i++){
@@ -993,7 +986,6 @@ void MElement::writeSTL(FILE *fp, bool binary, double scalingFactor)
void MElement::writePLY2(FILE *fp)
{
- setVolumePositive();
fprintf(fp, "3 ");
for(int i = 0; i < getNumVertices(); i++)
fprintf(fp, " %d", getVertex(i)->getIndex() - 1);
@@ -1002,7 +994,6 @@ void MElement::writePLY2(FILE *fp)
void MElement::writeVRML(FILE *fp)
{
- setVolumePositive();
for(int i = 0; i < getNumVertices(); i++)
fprintf(fp, "%d,", getVertex(i)->getIndex() - 1);
fprintf(fp, "-1,\n");
@@ -1012,8 +1003,6 @@ void MElement::writeVTK(FILE *fp, bool binary, bool bigEndian)
{
if(!getTypeForVTK()) return;
- setVolumePositive();
-
int n = getNumVertices();
if(binary){
int verts[60];
@@ -1037,7 +1026,6 @@ void MElement::writeUNV(FILE *fp, int num, int elementary, int physical)
int type = getTypeForUNV();
if(!type) return;
- setVolumePositive();
int n = getNumVertices();
int physical_property = elementary;
int material_property = abs(physical);
@@ -1063,7 +1051,6 @@ void MElement::writeUNV(FILE *fp, int num, int elementary, int physical)
void MElement::writeMESH(FILE *fp, int elementTagType, int elementary,
int physical)
{
- setVolumePositive();
for(int i = 0; i < getNumVertices(); i++)
if (getTypeForMSH() == MSH_TET_10 && i == 8)
fprintf(fp, " %d", getVertex(9)->getIndex());
@@ -1079,8 +1066,6 @@ void MElement::writeIR3(FILE *fp, int elementTagType, int num, int elementary,
int physical)
{
int numVert = getNumVertices();
- bool ok = setVolumePositive();
- if(getDim() == 3 && !ok) Msg::Error("Element %d has zero volume", num);
fprintf(fp, "%d %d %d", num, (elementTagType == 3) ? _partition :
(elementTagType == 2) ? physical : elementary, numVert);
for(int i = 0; i < numVert; i++)
@@ -1094,7 +1079,6 @@ void MElement::writeBDF(FILE *fp, int format, int elementTagType, int elementary
const char *str = getStringForBDF();
if(!str) return;
- setVolumePositive();
int n = getNumVertices();
const char *cont[4] = {"E", "F", "G", "H"};
int ncont = 0;
@@ -1135,8 +1119,6 @@ void MElement::writeDIFF(FILE *fp, int num, bool binary, int physical_property)
const char *str = getStringForDIFF();
if(!str) return;
- setVolumePositive();
-
int n = getNumVertices();
if(binary){
// TODO
@@ -1151,7 +1133,6 @@ void MElement::writeDIFF(FILE *fp, int num, bool binary, int physical_property)
void MElement::writeINP(FILE *fp, int num)
{
- setVolumePositive();
fprintf(fp, "%d, ", num);
int n = getNumVertices();
for(int i = 0; i < n; i++){
@@ -1166,7 +1147,6 @@ void MElement::writeINP(FILE *fp, int num)
void MElement::writeSU2(FILE *fp, int num)
{
- setVolumePositive();
fprintf(fp, "%d ", getTypeForVTK());
for(int i = 0; i < getNumVertices(); i++)
fprintf(fp, "%d ", getVertexVTK(i)->getIndex() - 1);
--
GitLab