From 36f1e7d39d9ab3b4670ba106d55060060691db53 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 20 Jun 2008 12:15:44 +0000
Subject: [PATCH] High order meshes improved
---
Geo/Geo.cpp | 17 +++++-
Geo/MElement.cpp | 120 ++++++++++++++++++------------------
Geo/MElement.h | 40 ++++++------
Geo/gmshFace.cpp | 6 +-
Mesh/HighOrder.cpp | 131 +++++++++++++++++++++++++++++++++++++---
Mesh/meshGRegion.cpp | 7 +--
Numeric/FunctionSpace.h | 2 +-
Numeric/gsl_newt.cpp | 8 +--
8 files changed, 228 insertions(+), 103 deletions(-)
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index ff8ac7b3e3..933434a1a1 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: Geo.cpp,v 1.115 2008-06-20 05:51:36 geuzaine Exp $
+// $Id: Geo.cpp,v 1.116 2008-06-20 12:15:44 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -2963,6 +2963,7 @@ bool ProjectPointOnCurve(Curve *c, Vertex *v, Vertex *RES, Vertex *DER)
bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2])
{
double x[3] = { 0.5, u[0], u[1] };
+ double res[3];
Vertex vv;
int check;
SURFACE = s;
@@ -2971,6 +2972,14 @@ bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2])
double UMAX = 1.;
double VMIN = 0.;
double VMAX = 1.;
+
+ newt(x, 2, &check, projectPS);
+ vv = InterpolateSurface(s, x[1], x[2], 0, 0);
+ projectPS(2, x, res);
+
+ double resid = sqrt (res[1]*res[1]+res[2]*res[2]);
+ /*
+
double eps = 1.e-6;
int nn = 0;
while(++nn < 20) {
@@ -2981,13 +2990,15 @@ bool ProjectPointOnSurface(Surface *s, Vertex &p, double u[2])
x[1] = UMIN + (UMAX - UMIN) * ((rand() % 10000) / 10000.);
x[2] = VMIN + (VMAX - VMIN) * ((rand() % 10000) / 10000.);
}
+ */
p.Pos.X = vv.Pos.X;
p.Pos.Y = vv.Pos.Y;
p.Pos.Z = vv.Pos.Z;
u[0] = x[1];
u[1] = x[2];
- if(!check) {
- return false;
+ if (resid > 1.e-6){
+ // printf("error %12.5E\n",resid);
+ return false;
}
return true;
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index a39f774418..659348380a 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.72 2008-06-10 12:59:12 remacle Exp $
+// $Id: MElement.cpp,v 1.73 2008-06-20 12:15:44 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -633,7 +633,7 @@ void MElement::writeBDF(FILE *fp, int format, int elementary)
}
}
-void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double j[2][3])
+void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double ww, double j[2][3])
{
#if defined(HAVE_GMSH_EMBEDDED)
return;
@@ -643,21 +643,21 @@ void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double j[2][3]
if (!nf){
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, grads); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, grads); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_9).df(uu, vv, grads); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_12).df(uu, vv, grads); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_15I).df(uu, vv, grads); break;
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, ww, grads); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, ww, grads); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_9).df(uu, vv, ww, grads); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_12).df(uu, vv, ww, grads); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_15I).df(uu, vv,ww, grads); break;
default: throw;
}
}
else{
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, grads); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, grads); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_10).df(uu, vv, grads); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_15).df(uu, vv, grads); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_21).df(uu, vv, grads); break;
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).df(uu, vv, ww,grads); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).df(uu, vv, ww,grads); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_10).df(uu, vv, ww,grads); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_15).df(uu, vv, ww,grads); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_21).df(uu, vv, ww,grads); break;
default: throw;
}
}
@@ -679,7 +679,7 @@ void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double j[2][3]
#endif
}
-void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, SPoint3 &p)
+void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, double ww, SPoint3 &p)
{
#if !defined(HAVE_GMSH_EMBEDDED)
double sf[256];
@@ -687,21 +687,21 @@ void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, SPoint3 &p)
if (!nf){
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv, sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv, sf); break;
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, ww,sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, ww,sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv, ww,sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv, ww,sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv, ww,sf); break;
default: throw;
}
}
else{
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv, sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv, sf); break;
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, ww,sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, ww,sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv, ww,sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv, ww,sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv, ww,sf); break;
default: throw;
}
}
@@ -718,16 +718,16 @@ void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, SPoint3 &p)
#endif
}
-void MTetrahedron::pnt(int ord, MVertex *vs[], double uu, double vv, SPoint3 &p)
+void MTetrahedron::pnt(int ord, MVertex *vs[], double uu, double vv, double ww,SPoint3 &p)
{
#if !defined(HAVE_GMSH_EMBEDDED)
double sf[256];
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TET_4).f(uu, vv, sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TET_10).f(uu, vv, sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TET_20).f(uu, vv, sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TET_35).f(uu, vv, sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TET_56).f(uu, vv, sf); break;
+ case 1: gmshFunctionSpaces::find(MSH_TET_4).f(uu, vv, ww,sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TET_10).f(uu, vv, ww,sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TET_20).f(uu, vv, ww,sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TET_35).f(uu, vv, ww,sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TET_56).f(uu, vv, ww,sf); break;
default : throw;
}
@@ -745,18 +745,18 @@ void MTetrahedron::pnt(int ord, MVertex *vs[], double uu, double vv, SPoint3 &p)
#endif
}
-void MTetrahedron::jac(int ord, MVertex *vs[], double uu, double vv, double j[3][3])
+void MTetrahedron::jac(int ord, MVertex *vs[], double uu, double vv, double ww, double j[3][3])
{
#if defined(HAVE_GMSH_EMBEDDED)
return;
#else
double grads[256][2];
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TET_4).df(uu, vv, grads); break;
- case 2: gmshFunctionSpaces::find(MSH_TET_10).df(uu, vv, grads); break;
- case 3: gmshFunctionSpaces::find(MSH_TET_20).df(uu, vv, grads); break;
- case 4: gmshFunctionSpaces::find(MSH_TET_35).df(uu, vv, grads); break;
- case 5: gmshFunctionSpaces::find(MSH_TET_56).df(uu, vv, grads); break;
+ case 1: gmshFunctionSpaces::find(MSH_TET_4).df(uu, vv, ww,grads); break;
+ case 2: gmshFunctionSpaces::find(MSH_TET_10).df(uu, vv, ww, grads); break;
+ case 3: gmshFunctionSpaces::find(MSH_TET_20).df(uu, vv, ww, grads); break;
+ case 4: gmshFunctionSpaces::find(MSH_TET_35).df(uu, vv, ww, grads); break;
+ case 5: gmshFunctionSpaces::find(MSH_TET_56).df(uu, vv, ww, grads); break;
default: throw;
}
@@ -798,8 +798,8 @@ void MTriangle6::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
int N = getNumEdgesRep() / 3;
if (num < N){
SPoint3 pnt1, pnt2;
- pnt((double)num / N, 0., pnt1);
- pnt((double)(num + 1) / N, 0., pnt2);
+ pnt((double)num / N, 0., 0,pnt1);
+ pnt((double)(num + 1) / N, 0., 0,pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
@@ -808,8 +808,8 @@ void MTriangle6::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
if (num < 2 * N){
SPoint3 pnt1, pnt2;
num -= N;
- pnt(1. - (double)num / N, (double)num / N, pnt1);
- pnt(1. - (double)(num + 1) / N, (double)(num + 1) / N, pnt2);
+ pnt(1. - (double)num / N, (double)num / N, 0,pnt1);
+ pnt(1. - (double)(num + 1) / N, (double)(num + 1) / N, 0,pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
@@ -818,15 +818,15 @@ void MTriangle6::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
{
SPoint3 pnt1, pnt2;
num -= 2 * N;
- pnt(0, (double)num / N, pnt1);
- pnt(0, (double)(num + 1) / N, pnt2);
+ pnt(0, (double)num / N, 0,pnt1);
+ pnt(0, (double)(num + 1) / N, 0,pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
}
}
-const int numSubEdges = 6;
+const int numSubEdges = 12;
int MTriangleN::getNumFacesRep(){ return numSubEdges * numSubEdges; }
@@ -853,20 +853,20 @@ void MTriangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *
SPoint3 pnt1, pnt2, pnt3;
double J1[2][3],J2[2][3],J3[2][3];
if (ix %2 == 0){
- pnt(ix/2*d, iy*d, pnt1);
- pnt((ix/2+1)*d, iy*d, pnt2);
- pnt(ix/2*d, (iy+1)*d, pnt3);
- jac(ix/2*d, iy*d, J1);
- jac((ix/2+1)*d, iy*d, J2);
- jac(ix/2*d, (iy+1)*d, J3);
+ pnt(ix/2*d, iy*d, 0,pnt1);
+ pnt((ix/2+1)*d, iy*d, 0,pnt2);
+ pnt(ix/2*d, (iy+1)*d, 0,pnt3);
+ jac(ix/2*d, iy*d, 0,J1);
+ jac((ix/2+1)*d, iy*d, 0,J2);
+ jac(ix/2*d, (iy+1)*d, 0,J3);
}
else{
- pnt((ix/2+1)*d, iy*d, pnt1);
- pnt((ix/2+1)*d, (iy+1)*d, pnt2);
- pnt(ix/2*d, (iy+1)*d, pnt3);
- jac((ix/2+1)*d, iy*d, J1);
- jac((ix/2+1)*d, (iy+1)*d, J2);
- jac(ix/2*d, (iy+1)*d, J3);
+ pnt((ix/2+1)*d, iy*d, 0,pnt1);
+ pnt((ix/2+1)*d, (iy+1)*d, 0,pnt2);
+ pnt(ix/2*d, (iy+1)*d, 0,pnt3);
+ jac((ix/2+1)*d, iy*d, 0,J1);
+ jac((ix/2+1)*d, (iy+1)*d, 0,J2);
+ jac(ix/2*d, (iy+1)*d, 0,J3);
}
{
SVector3 d1 (J1[0][0],J1[0][1],J1[0][2]);
@@ -901,8 +901,8 @@ void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
int N = getNumEdgesRep() / 3;
if (num < N){
SPoint3 pnt1, pnt2;
- pnt((double)num / N, 0., pnt1);
- pnt((double)(num + 1) / N, 0., pnt2);
+ pnt((double)num / N, 0., 0,pnt1);
+ pnt((double)(num + 1) / N, 0., 0,pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
@@ -911,8 +911,8 @@ void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
if (num < 2 * N){
SPoint3 pnt1, pnt2;
num -= N;
- pnt(1. - (double)num / N, (double)num / N, pnt1);
- pnt(1. - (double)(num + 1) / N, (double)(num + 1) / N, pnt2);
+ pnt(1. - (double)num / N, (double)num / N, 0,pnt1);
+ pnt(1. - (double)(num + 1) / N, (double)(num + 1) / N, 0,pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
@@ -921,8 +921,8 @@ void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *
{
SPoint3 pnt1, pnt2;
num -= 2 * N;
- pnt(0, (double)num / N, pnt1);
- pnt(0, (double)(num + 1) / N, pnt2);
+ pnt(0, (double)num / N, 0,pnt1);
+ pnt(0, (double)(num + 1) / N, 0,pnt2);
x[0] = pnt1.x(); x[1] = pnt2.x();
y[0] = pnt1.y(); y[1] = pnt2.y();
z[0] = pnt1.z(); z[1] = pnt2.z();
diff --git a/Geo/MElement.h b/Geo/MElement.h
index e25d730ff5..707c9ab7c6 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -444,15 +444,15 @@ class MTriangle : public MElement {
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
- virtual void jac(int order, MVertex *verts[], double u, double v, double jac[2][3]);
- virtual void jac(double u, double v, double j[2][3])
+ virtual void jac(int order, MVertex *verts[], double u, double v, double w, double jac[2][3]);
+ virtual void jac(double u, double v, double w, double j[2][3])
{
- jac(1, 0, u, v, j);
+ jac(1, 0, u, v, w, j);
}
- virtual void pnt(int order, MVertex *verts[], double u, double v, SPoint3 &);
- virtual void pnt(double u, double v, SPoint3 &p)
+ virtual void pnt(int order, MVertex *verts[], double u, double v, double w, SPoint3 &);
+ virtual void pnt(double u, double v, double w, SPoint3 &p)
{
- pnt(1, 0, u, v, p);
+ pnt(1, 0, u, v, w, p);
}
virtual void getShapeFunction(int num, double u, double v, double w, double &s)
{
@@ -543,13 +543,13 @@ class MTriangle6 : public MTriangle {
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp;
}
- virtual void jac(double u, double v, double j[2][3])
+ virtual void jac(double u, double v, double w, double j[2][3])
{
- MTriangle::jac(2, _vs, u, v, j);
+ MTriangle::jac(2, _vs, u, v, w, j);
}
- virtual void pnt(double u, double v, SPoint3 &p)
+ virtual void pnt(double u, double v, double w, SPoint3 &p)
{
- MTriangle::pnt(2, _vs, u, v, p);
+ MTriangle::pnt(2, _vs, u, v, w, p);
}
};
@@ -616,13 +616,13 @@ class MTriangleN : public MTriangle {
inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
_vs = inv;
}
- virtual void jac(double u, double v, double j[2][3])
+ virtual void jac(double u, double v, double w, double j[2][3])
{
- MTriangle::jac(_order, &(*(_vs.begin())), u, v, j);
+ MTriangle::jac(_order, &(*(_vs.begin())), u, v, w, j);
}
- virtual void pnt(double u, double v, SPoint3 &p)
+ virtual void pnt(double u, double v, double w, SPoint3 &p)
{
- MTriangle::pnt(_order, &(*(_vs.begin())), u, v, p);
+ MTriangle::pnt(_order, &(*(_vs.begin())), u, v, w, p);
}
};
@@ -933,8 +933,8 @@ class MTetrahedron : public MElement {
default : s = 0.; break;
}
}
- virtual void jac(int order, MVertex *verts[], double u, double v, double jac[3][3]);
- virtual void pnt(int order, MVertex *verts[], double u, double v, SPoint3 &);
+ virtual void jac(int order, MVertex *verts[], double u, double v, double w, double jac[3][3]);
+ virtual void pnt(int order, MVertex *verts[], double u, double v, double w, SPoint3 &);
virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])
{
switch(num) {
@@ -1081,13 +1081,13 @@ class MTetrahedronN : public MTetrahedron {
inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
_vs = inv;
}
- virtual void jac(double u, double v, double j[3][3])
+ virtual void jac(double u, double v, double w, double j[3][3])
{
- MTetrahedron::jac(_order, &(*(_vs.begin())), u, v, j);
+ MTetrahedron::jac(_order, &(*(_vs.begin())), u, v, w, j);
}
- virtual void pnt(double u, double v, SPoint3 &p)
+ virtual void pnt(double u, double v, double w, SPoint3 &p)
{
- MTetrahedron::pnt(_order, &(*(_vs.begin())), u, v, p);
+ MTetrahedron::pnt(_order, &(*(_vs.begin())), u, v, w, p);
}
};
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 9897d25ba6..4c0f220c05 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.62 2008-06-20 11:50:00 geuzaine Exp $
+// $Id: gmshFace.cpp,v 1.63 2008-06-20 12:15:44 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -207,7 +207,9 @@ GPoint gmshFace::closestPoint(const SPoint3 & qp, const double initialGuess[2])
v.Pos.X = qp.x();
v.Pos.Y = qp.y();
v.Pos.Z = qp.z();
- ProjectPointOnSurface(s, v, u);
+ bool result = ProjectPointOnSurface(s, v, u);
+ if (!result)
+ printf("Project Point on surface %d \n",result);
return GPoint(v.Pos.X, v.Pos.Y, v.Pos.Z, this, u);
}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 70eb301412..5d5096205b 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.30 2008-06-10 08:37:34 remacle Exp $
+// $Id: HighOrder.cpp,v 1.31 2008-06-20 12:15:44 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -292,7 +292,7 @@ void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
relax /= 2.0;
if (relax < 1.e-2)break;
}
- if (relax < 1.e-2)printf("failure in computing equidistant parameters %g\n",relax);
+ if (relax < 1.e-2)Msg::Warning("failure in computing equidistant parameters %g",relax);
}
std::vector<MVertex*> temp;
for(int j = 0; j < nPts; j++){
@@ -401,6 +401,114 @@ void getEdgeVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &ve,
}
}
+void getFaceVertices(GFace *gf,
+ MElement *incomplete,
+ MElement *ele,
+ std::vector<MVertex*> &vf,
+ faceContainer &faceVertices,
+ bool linear, int nPts = 1){
+ Double_Matrix points;
+ int start = 0;
+
+ switch (nPts){
+ case 2 :
+ points = gmshFunctionSpaces::find(MSH_TRI_10).points;
+ start = 9;
+ break;
+ case 3 :
+ points = gmshFunctionSpaces::find(MSH_TRI_15).points;
+ start = 12;
+ break;
+ case 4 :
+ points = gmshFunctionSpaces::find(MSH_TRI_21).points;
+ start = 15;
+ break;
+ default :
+ // do nothing (e.g. for quad faces)
+ break;
+ }
+
+ for(int i = 0; i < ele->getNumFaces(); i++){
+ MFace face = ele->getFace(i);
+ std::vector<MVertex*> p;
+ face.getOrderedVertices(p);
+ if(faceVertices.count(p)){
+ vf.insert(vf.end(), faceVertices[p].begin(), faceVertices[p].end());
+ }
+ else{
+ SPoint2 pts[20];
+ bool reparamOK = true;
+ if(!linear &&
+ gf->geomType() != GEntity::DiscreteSurface &&
+ gf->geomType() != GEntity::BoundaryLayerSurface){
+ for (int k=0;k<incomplete->getNumVertices(); k++){
+ reparamOK &= reparamOnFace(incomplete->getVertex(k), gf, pts[k]);
+ }
+ }
+ if(face.getNumVertices() == 3){ // triangles
+ for(int k = start ; k < points.size1() ; k++){
+ MVertex *v;
+ const double t1 = points(k, 0);
+ const double t2 = points(k, 1);
+ if(!reparamOK || linear || gf->geomType() == GEntity::DiscreteSurface){
+ SPoint3 pc = face.interpolate(t1, t2);
+ v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
+ }
+ else{
+ double X(0),Y(0),Z(0),GUESS[2]={0,0};
+ for (int j=0; j<incomplete->getNumVertices(); j++){
+ double sf ; incomplete->getShapeFunction(j,t1,t2,0,sf);
+ MVertex *vt = incomplete->getVertex(j);
+ X += sf * vt->x();
+ Y += sf * vt->y();
+ Z += sf * vt->z();
+ GUESS[0] += sf * pts[j][0];
+ GUESS[1] += sf * pts[j][1];
+ }
+ GPoint gp = gf->closestPoint(SPoint3(X,Y,Z),GUESS);
+ // printf("%g %g %g -- %g %g %g\n",X,Y,Z,gp.x(),gp.y(),gp.z());
+ v = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
+ }
+ faceVertices[p].push_back(v);
+ gf->mesh_vertices.push_back(v);
+ vf.push_back(v);
+ }
+ }
+ else if(face.getNumVertices() == 4){ // quadrangles
+ for(int j = 0; j < nPts; j++){
+ for(int k = 0; k < nPts; k++){
+ MVertex *v;
+ // parameters are between -1 and 1
+ double t1 = 2. * (double)(j + 1) / (nPts + 1) - 1.;
+ double t2 = 2. * (double)(k + 1) / (nPts + 1) - 1.;
+ if(!reparamOK || linear || gf->geomType() == GEntity::DiscreteSurface){
+ SPoint3 pc = face.interpolate(t1, t2);
+ v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
+ }
+ else{
+ double uc = 0.25 * ((1 - t1) * (1 - t2) * pts[0][0] +
+ (1 + t1) * (1 - t2) * pts[1][0] +
+ (1 + t1) * (1 + t2) * pts[2][0] +
+ (1 - t1) * (1 + t2) * pts[3][0]);
+ double vc = 0.25 * ((1 - t1) * (1 - t2) * pts[0][1] +
+ (1 + t1) * (1 - t2) * pts[1][1] +
+ (1 + t1) * (1 + t2) * pts[2][1] +
+ (1 - t1) * (1 + t2) * pts[3][1]);
+ GPoint pc = gf->point(uc, vc);
+ v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, uc, vc);
+ }
+ faceVertices[p].push_back(v);
+ gf->mesh_vertices.push_back(v);
+ vf.push_back(v);
+ }
+ }
+ }
+ else throw; // not tri or quad
+ }
+ }
+}
+
+
void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
faceContainer &faceVertices, bool linear, int nPts = 1)
{
@@ -582,7 +690,14 @@ void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
ve, nPts + 1));
}
else{
- getFaceVertices(gf, t, vf, faceVertices, linear, nPts);
+ if (gf->geomType() == GEntity::Plane){
+ getFaceVertices(gf, t, vf, faceVertices, linear, nPts);
+ }
+ else{
+ MTriangleN incpl (t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ ve, nPts + 1);
+ getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts);
+ }
ve.insert(ve.end(), vf.begin(), vf.end());
triangles2.push_back
(new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
@@ -798,7 +913,7 @@ static double mesh_functional_distorsion(MTriangle *t, double u, double v)
{
// compute uncurved element jacobian d_u x and d_v x
double mat[2][3];
- t->jac(1, 0, 0, 0, mat);
+ t->jac(1, 0, 0, 0, 0,mat);
double v1[3] = {mat[0][0], mat[0][1], mat[0][2]};
double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
double normal1[3];
@@ -806,7 +921,7 @@ static double mesh_functional_distorsion(MTriangle *t, double u, double v)
double nn = sqrt(DSQR(normal1[0]) + DSQR(normal1[1]) + DSQR(normal1[2]));
// compute uncurved element jacobian d_u x and d_v x
- t->jac(u, v, mat);
+ t->jac(u, v, 0,mat);
double v1b[3] = {mat[0][0], mat[0][1], mat[0][2]};
double v2b[3] = {mat[1][0], mat[1][1], mat[1][2]};
double normal[3];
@@ -825,7 +940,7 @@ void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
{
double mat[2][3];
int n = 3;
- t->jac(1, 0, 0, 0, mat);
+ t->jac(1, 0, 0, 0, 0,mat);
double v1[3] = {mat[0][0], mat[0][1], mat[0][2]};
double v2[3] = {mat[1][0], mat[1][1], mat[1][2]};
double normal1[3], normal[3];
@@ -833,7 +948,7 @@ void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
double nn = sqrt(DSQR(normal1[0]) + DSQR(normal1[1]) + DSQR(normal1[2]));
for(int i = 0; i < n; i++){
for(int k = 0; k < n - i; k++){
- t->jac((double)i / (n - 1), (double)k / (n - 1), mat);
+ t->jac((double)i / (n - 1), (double)k / (n - 1), 0,mat);
double v1b[3] = {mat[0][0], mat[0][1], mat[0][2]};
double v2b[3] = {mat[1][0], mat[1][1], mat[1][2]};
prodve(v1b, v2b, normal);
@@ -1256,7 +1371,7 @@ void printJacobians(GModel *m, const char *nm)
SPoint3 pt;
double u = (double)i / (n - 1);
double v = (double)k / (n - 1);
- t->pnt(u,v, pt);
+ t->pnt(u,v,0, pt);
D[i][k] = mesh_functional_distorsion (t,u,v);
X[i][k] = pt.x();
Y[i][k] = pt.y();
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 09b80552a0..3a434049a4 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegion.cpp,v 1.53 2008-06-20 11:50:00 geuzaine Exp $
+// $Id: meshGRegion.cpp,v 1.54 2008-06-20 12:15:44 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -194,10 +194,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
v[j]->x(), v[j]->y(), v[j]->z(),gp.x(),gp.y(),gp.z());
// To be safe, we should ensure that this mesh motion does not lead to an invalid mesh !!!!
- v[j]->x() = gp.x();
- v[j]->y() = gp.y();
- v[j]->z() = gp.z();
- v1b = new MFaceVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf,gp.u(),gp.v());
+ v1b = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
}
else{
v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(),gf);
diff --git a/Numeric/FunctionSpace.h b/Numeric/FunctionSpace.h
index b568ea9972..267347278d 100644
--- a/Numeric/FunctionSpace.h
+++ b/Numeric/FunctionSpace.h
@@ -80,7 +80,7 @@ struct gmshFunctionSpace
}
}
}
- inline void df(double u, double v, double grads[][2]) const
+ inline void df(double u, double v, double w,double grads[][2]) const
{
double p[256][2];
computePows(u, v, p);
diff --git a/Numeric/gsl_newt.cpp b/Numeric/gsl_newt.cpp
index b95233f572..7194d47cdd 100644
--- a/Numeric/gsl_newt.cpp
+++ b/Numeric/gsl_newt.cpp
@@ -1,4 +1,4 @@
-// $Id: gsl_newt.cpp,v 1.19 2008-05-04 08:31:16 geuzaine Exp $
+// $Id: gsl_newt.cpp,v 1.20 2008-06-20 12:15:44 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -37,7 +37,7 @@
#include <gsl/gsl_linalg.h>
#define MAX_DIM_NEWT 10
-#define MAXITER 10000
+#define MAXITER 100
#define PREC 1.e-8
int nrdim;
@@ -101,8 +101,8 @@ void newt(double x[], int n, int *check,
status = gsl_multiroot_fsolver_iterate(s);
// Msg::Info("status %d %d %d %lf %lf\n",
// status,n,iter,gsl_vector_get(s->x,0),gsl_vector_get(s->x,1));
- if(status)
- break; // solver problem
+ if(status)
+ break; // solver problem
status = gsl_multiroot_test_residual(s->f, n * PREC);
}
while(status == GSL_CONTINUE && iter < MAXITER);
--
GitLab