diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 62802a05a72e1afb2f5c7fe1fb70010bbcfcf72c..207e81cb7cacc0e0580191dbdd585d0761670a0a 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1589,7 +1589,7 @@ int GModel::readMESH(const std::string &name)
}
std::vector<MVertex*> vertexVector;
- std::map<int, std::vector<MElement*> > elements[4];
+ std::map<int, std::vector<MElement*> > elements[5];
while(!feof(fp)) {
if(!fgets(buffer, 256, fp)) break;
@@ -1625,7 +1625,7 @@ int GModel::readMESH(const std::string &name)
for(int j = 0; j < 2; j++) n[j]--;
std::vector<MVertex*> vertices;
if(!getVertices(2, n, vertexVector, vertices)) return 0;
- elements[3][cl].push_back(new MLine(vertices));
+ elements[0][cl].push_back(new MLine(vertices));
}
}
else if(!strcmp(str, "Triangles")){
@@ -1640,7 +1640,7 @@ int GModel::readMESH(const std::string &name)
for(int j = 0; j < 3; j++) n[j]--;
std::vector<MVertex*> vertices;
if(!getVertices(3, n, vertexVector, vertices)) return 0;
- elements[0][cl].push_back(new MTriangle(vertices));
+ elements[1][cl].push_back(new MTriangle(vertices));
}
}
else if(!strcmp(str, "Quadrilaterals")) {
@@ -1655,7 +1655,7 @@ int GModel::readMESH(const std::string &name)
for(int j = 0; j < 4; j++) n[j]--;
std::vector<MVertex*> vertices;
if(!getVertices(4, n, vertexVector, vertices)) return 0;
- elements[1][cl].push_back(new MQuadrangle(vertices));
+ elements[2][cl].push_back(new MQuadrangle(vertices));
}
}
else if(!strcmp(str, "Tetrahedra")) {
@@ -1670,7 +1670,23 @@ int GModel::readMESH(const std::string &name)
for(int j = 0; j < 4; j++) n[j]--;
std::vector<MVertex*> vertices;
if(!getVertices(4, n, vertexVector, vertices)) return 0;
- elements[2][cl].push_back(new MTetrahedron(vertices));
+ elements[3][cl].push_back(new MTetrahedron(vertices));
+ }
+ }
+ else if(!strcmp(str, "Hexahedra")) {
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int nbe;
+ sscanf(buffer, "%d", &nbe);
+ Msg::Info("%d hexahedra", nbe);
+ for(int i = 0; i < nbe; i++) {
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int n[8], cl;
+ sscanf(buffer, "%d %d %d %d %d %d %d %d %d", &n[0], &n[1], &n[2], &n[3],
+ &n[4], &n[5], &n[6], &n[7], &cl);
+ for(int j = 0; j < 8; j++) n[j]--;
+ std::vector<MVertex*> vertices;
+ if(!getVertices(8, n, vertexVector, vertices)) return 0;
+ elements[4][cl].push_back(new MHexahedron(vertices));
}
}
}
@@ -1710,7 +1726,8 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
entities[i]->mesh_vertices[j]->writeMESH(fp, scalingFactor);
- int numEdges = 0, numTriangles = 0, numQuadrangles = 0, numTetrahedra = 0;
+ int numEdges = 0, numTriangles = 0, numQuadrangles = 0;
+ int numTetrahedra = 0, numHexahedra = 0;
for(eiter it = firstEdge(); it != lastEdge(); ++it){
if(saveAll || (*it)->physicals.size()){
numEdges += (*it)->lines.size();
@@ -1725,6 +1742,7 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
for(riter it = firstRegion(); it != lastRegion(); ++it){
if(saveAll || (*it)->physicals.size()){
numTetrahedra += (*it)->tetrahedra.size();
+ numHexahedra += (*it)->hexahedra.size();
}
}
@@ -1776,6 +1794,18 @@ int GModel::writeMESH(const std::string &name, int elementTagType,
}
}
}
+ if(numHexahedra){
+ fprintf(fp, " Hexahedra\n");
+ fprintf(fp, " %d\n", numHexahedra);
+ for(riter it = firstRegion(); it != lastRegion(); ++it){
+ int numPhys = (*it)->physicals.size();
+ if(saveAll || numPhys){
+ for(unsigned int i = 0; i < (*it)->hexahedra.size(); i++)
+ (*it)->hexahedra[i]->writeMESH(fp, elementTagType, (*it)->tag(),
+ numPhys ? (*it)->physicals[0] : 0);
+ }
+ }
+ }
fprintf(fp, " End\n");
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index faf16cd95963c0450bc2f9c0948ef2d3caf0904f..7ff6329166283d61d070cfdc07db67d16e403f43 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -351,79 +351,78 @@ void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv
void planarQuad_xyz2xy(double *x, double *y, double *z, double *xn, double *yn)
{
- double v1[3] = {x[1] - x[0], y[1] - y[0], z[1] - z[0]};
- double v2[3] = {x[2] - x[0], y[2] - y[0], z[2] - z[0]};
- double v3[3] = {x[3] - x[0], y[3] - y[0], z[3] - z[0]};
-
- double vx[3] = {x[1] - x[0], y[1] - y[0], z[1] - z[0]};
- double vy[3] = {x[2] - x[0], y[2] - y[0], z[2] - z[0]};
- double vz[3]; prodve(vx, vy, vz); prodve(vz, vx, vy);
-
- norme(vx); norme(vy); norme(vz);
-
- double p1P[2] = {0.0, 0.0};
- double p2P[2]; prosca(v1, vx, &p2P[0]); prosca(v1, vy, &p2P[1]);
- double p3P[2]; prosca(v2, vx, &p3P[0]); prosca(v2, vy, &p3P[1]);
- double p4P[2]; prosca(v3, vx, &p4P[0]); prosca(v3, vy, &p4P[1]);
-
- xn[0] = p1P[0];
- xn[1] = p2P[0];
- xn[2] = p3P[0];
- xn[3] = p4P[0];
- yn[0] = p1P[1];
- yn[1] = p2P[1];
- yn[2] = p3P[1];
- yn[3] = p4P[1];
+ double v1[3] = {x[1] - x[0], y[1] - y[0], z[1] - z[0]};
+ double v2[3] = {x[2] - x[0], y[2] - y[0], z[2] - z[0]};
+ double v3[3] = {x[3] - x[0], y[3] - y[0], z[3] - z[0]};
+
+ double vx[3] = {x[1] - x[0], y[1] - y[0], z[1] - z[0]};
+ double vy[3] = {x[2] - x[0], y[2] - y[0], z[2] - z[0]};
+ double vz[3]; prodve(vx, vy, vz); prodve(vz, vx, vy);
+
+ norme(vx); norme(vy); norme(vz);
+
+ double p1P[2] = {0.0, 0.0};
+ double p2P[2]; prosca(v1, vx, &p2P[0]); prosca(v1, vy, &p2P[1]);
+ double p3P[2]; prosca(v2, vx, &p3P[0]); prosca(v2, vy, &p3P[1]);
+ double p4P[2]; prosca(v3, vx, &p4P[0]); prosca(v3, vy, &p4P[1]);
+
+ xn[0] = p1P[0];
+ xn[1] = p2P[0];
+ xn[2] = p3P[0];
+ xn[3] = p4P[0];
+ yn[0] = p1P[1];
+ yn[1] = p2P[1];
+ yn[2] = p3P[1];
+ yn[3] = p4P[1];
}
-double computeInnerRadiusForQuad(double *x, double *y, int i){
-
- // parameters of the equations of the 3 edges
- double a1 = y[(4+i)%4]-y[(5+i)%4];
- double a2 = y[(5+i)%4]-y[(6+i)%4];
- double a3 = y[(6+i)%4]-y[(7+i)%4];
-
- double b1 = x[(5+i)%4]-x[(4+i)%4];
- double b2 = x[(6+i)%4]-x[(5+i)%4];
- double b3 = x[(7+i)%4]-x[(6+i)%4];
-
- double c1 = y[(5+i)%4]*x[(4+i)%4]-y[(4+i)%4]*x[(5+i)%4];
- double c2 = y[(6+i)%4]*x[(5+i)%4]-y[(5+i)%4]*x[(6+i)%4];
- double c3 = y[(7+i)%4]*x[(6+i)%4]-y[(6+i)%4]*x[(7+i)%4];
-
- // length of the 3 edges
- double l1 = sqrt(a1*a1+b1*b1);
- double l2 = sqrt(a2*a2+b2*b2);
- double l3 = sqrt(a3*a3+b3*b3);
-
- // parameters of the 2 bisectors
- double a12 = a1/l1-a2/l2;
- double a23 = a2/l2-a3/l3;
-
- double b12 = b1/l1-b2/l2;
- double b23 = b2/l2-b3/l3;
-
- double c12 = c1/l1-c2/l2;
- double c23 = c2/l2-c3/l3;
-
- // compute the coordinates of the center of the incircle,
- // that is the point where the 2 bisectors meet
- double x_s = (c12*b23-c23*b12)/(a23*b12-a12*b23);
- double y_s = 0.;
- if (b12 != 0) {
- y_s = -a12/b12*x_s-c12/b12;
- }
- else {
- y_s = -a23/b23*x_s-c23/b23;
- }
-
- // finally get the radius of the circle
- double r = (a1*x_s+b1*y_s+c1)/l1;
-
- return r;
+double computeInnerRadiusForQuad(double *x, double *y, int i)
+{
+ // parameters of the equations of the 3 edges
+ double a1 = y[(4+i)%4]-y[(5+i)%4];
+ double a2 = y[(5+i)%4]-y[(6+i)%4];
+ double a3 = y[(6+i)%4]-y[(7+i)%4];
+
+ double b1 = x[(5+i)%4]-x[(4+i)%4];
+ double b2 = x[(6+i)%4]-x[(5+i)%4];
+ double b3 = x[(7+i)%4]-x[(6+i)%4];
+
+ double c1 = y[(5+i)%4]*x[(4+i)%4]-y[(4+i)%4]*x[(5+i)%4];
+ double c2 = y[(6+i)%4]*x[(5+i)%4]-y[(5+i)%4]*x[(6+i)%4];
+ double c3 = y[(7+i)%4]*x[(6+i)%4]-y[(6+i)%4]*x[(7+i)%4];
+
+ // length of the 3 edges
+ double l1 = sqrt(a1*a1+b1*b1);
+ double l2 = sqrt(a2*a2+b2*b2);
+ double l3 = sqrt(a3*a3+b3*b3);
+
+ // parameters of the 2 bisectors
+ double a12 = a1/l1-a2/l2;
+ double a23 = a2/l2-a3/l3;
+
+ double b12 = b1/l1-b2/l2;
+ double b23 = b2/l2-b3/l3;
+
+ double c12 = c1/l1-c2/l2;
+ double c23 = c2/l2-c3/l3;
+
+ // compute the coordinates of the center of the incircle,
+ // that is the point where the 2 bisectors meet
+ double x_s = (c12*b23-c23*b12)/(a23*b12-a12*b23);
+ double y_s = 0.;
+ if (b12 != 0) {
+ y_s = -a12/b12*x_s-c12/b12;
+ }
+ else {
+ y_s = -a23/b23*x_s-c23/b23;
+ }
+
+ // finally get the radius of the circle
+ double r = (a1*x_s+b1*y_s+c1)/l1;
+
+ return r;
}
-
char float2char(float f)
{
// float normalized in [-1, 1], char in [-127, 127]
@@ -686,18 +685,18 @@ void gmshLineSearch(double (*func)(fullVector<double> &, void *), void* data,
const double TOLX = 1.0e-9;
fullVector<double> xold(x);
- const double fold = (*func)(xold,data);
+ const double fold = (*func)(xold, data);
check=0;
int n = x.size();
double norm = p.norm();
- if (norm > stpmax)p.scale(stpmax/norm);
+ if (norm > stpmax) p.scale(stpmax / norm);
double slope=0.0;
- for (i=0;i<n;i++)slope += g(i)*p(i);
+ for (i = 0; i < n; i++) slope += g(i)*p(i);
double test=0.0;
- for (i=0;i<n;i++) {
- temp=fabs(p(i))/std::max(fabs(xold(i)),1.0);
- if (temp > test) test=temp;
+ for (i = 0; i < n; i++) {
+ temp = fabs(p(i)) / std::max(fabs(xold(i)), 1.0);
+ if (temp > test) test = temp;
}
/*
for (int j=0;j<100;j++){
@@ -707,41 +706,41 @@ void gmshLineSearch(double (*func)(fullVector<double> &, void *), void* data,
}
*/
- alamin=TOLX/test;
- alam=1.0;
+ alamin = TOLX / test;
+ alam = 1.0;
while(1) {
- for (i=0;i<n;i++) x(i)=xold(i)+alam*p(i);
- f=(*func)(x,data);
+ for (i = 0; i < n; i++) x(i) = xold(i) + alam*p(i);
+ f = (*func)(x, data);
// printf("f = %g x = %g %g alam = %g p = %g %g\n",f,x(0),x(1),alam,p(0),p(1));
if (alam < alamin) {
- for (i=0;i<n;i++) x(i)=xold(i);
+ for (i = 0; i <n; i++) x(i) = xold(i);
// printf("ALERT : alam %g alamin %g\n",alam,alamin);
- check=1;
+ check = 1;
return;
}
- else if (f <= fold + ALF*alam*slope) return;
+ else if (f <= fold + ALF * alam * slope) return;
else {
if (alam == 1.0)
- tmplam = -slope/(2.0*(f-fold-slope));
+ tmplam = -slope / (2.0 * (f - fold - slope));
else {
- rhs1 = f-fold-alam*slope;
- rhs2=f2-fold2-alam2*slope;
- const double a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
- const double b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
- if (a == 0.0) tmplam = -slope/(2.0*b);
+ rhs1 = f - fold - alam * slope;
+ rhs2 = f2 - fold2 - alam2 * slope;
+ const double a = (rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
+ const double b = (-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
+ if (a == 0.0) tmplam = -slope / (2.0 * b);
else {
- const double disc=b*b-3.0*a*slope;
- if (disc<0.0) Msg::Error("Roundoff problem in gmshLineSearch.");
- else tmplam=(-b+sqrt(disc))/(3.0*a);
+ const double disc = b*b-3.0*a*slope;
+ if (disc < 0.0) Msg::Error("Roundoff problem in gmshLineSearch.");
+ else tmplam = (-b+sqrt(disc))/(3.0*a);
}
- if (tmplam>0.5*alam)
- tmplam=0.5*alam;
+ if (tmplam > 0.5 * alam)
+ tmplam = 0.5 * alam;
}
}
- alam2=alam;
+ alam2 = alam;
f2 = f;
- fold2=fold;
- alam=std::max(tmplam,0.1*alam);
+ fold2 = fold;
+ alam = std::max(tmplam, 0.1 * alam);
}
}
@@ -986,7 +985,8 @@ void changeReferential(const int direction,const SPoint3 &p,const SPoint3 &close
}
int computeDistanceRatio(const double &y, const double &yp, const double &x,
- const double &xp, double *distance, const double &r1, const double &r2)
+ const double &xp, double *distance, const double &r1,
+ const double &r2)
{
double b;
double a;