From 03689a7702fe83c2013f78c612beaf214e56c11f Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 19 Nov 2008 21:55:00 +0000
Subject: [PATCH] more work on high order
---
Geo/MElement.cpp | 109 +++++++++++++++++++++++++----------------------
Geo/MVertex.cpp | 23 +++++++---
2 files changed, 76 insertions(+), 56 deletions(-)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 5d340647af..7fe335040a 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -814,12 +814,52 @@ const gmshFunctionSpace* MTriangle::getFunctionSpace(int o) const
return 0;
}
-const int numSubEdges = 12;
+#define NUM_SUB_EDGES (_order)
-int MTriangleN::getNumFacesRep(){ return numSubEdges * numSubEdges; }
+int MTriangleN::getNumEdgesRep(){ return 3 * NUM_SUB_EDGES; }
+
+void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ int numSubEdges = NUM_SUB_EDGES;
+
+ n[0] = n[1] = n[2] = getFace(0).normal();
+
+ if (num < numSubEdges){
+ SPoint3 pnt1, pnt2;
+ pnt((double)num / numSubEdges, 0., 0.,pnt1);
+ pnt((double)(num + 1) / numSubEdges, 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();
+ return;
+ }
+ if (num < 2 * numSubEdges){
+ SPoint3 pnt1, pnt2;
+ num -= numSubEdges;
+ pnt(1. - (double)num / numSubEdges, (double)num / numSubEdges, 0, pnt1);
+ pnt(1. - (double)(num + 1) / numSubEdges, (double)(num + 1) / numSubEdges, 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();
+ return ;
+ }
+ {
+ SPoint3 pnt1, pnt2;
+ num -= 2 * numSubEdges;
+ pnt(0, (double)num / numSubEdges, 0,pnt1);
+ pnt(0, (double)(num + 1) / numSubEdges, 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();
+ }
+}
+
+int MTriangleN::getNumFacesRep(){ return NUM_SUB_EDGES * NUM_SUB_EDGES; }
void MTriangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
{
+ int numSubEdges = NUM_SUB_EDGES;
+
// on the first layer, we have (numSubEdges-1) * 2 + 1 triangles
// on the second layer, we have (numSubEdges-2) * 2 + 1 triangles
// on the ith layer, we have (numSubEdges-1-i) * 2 + 1 triangles
@@ -879,42 +919,6 @@ void MTriangleN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *
z[0] = pnt1.z(); z[1] = pnt2.z(); z[2] = pnt3.z();
}
-int MTriangleN::getNumEdgesRep(){ return 3 * numSubEdges; }
-
-void MTriangleN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
-{
- n[0] = n[1] = n[2] = getFace(0).normal();
- int N = getNumEdgesRep() / 3;
- if (num < N){
- SPoint3 pnt1, 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();
- return;
- }
- if (num < 2 * N){
- SPoint3 pnt1, pnt2;
- num -= N;
- 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();
- return ;
- }
- {
- SPoint3 pnt1, pnt2;
- num -= 2 * N;
- 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();
- }
-}
-
void MTriangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
{
#if !defined(HAVE_GMSH_EMBEDDED)
@@ -1033,40 +1037,43 @@ const gmshFunctionSpace* MTetrahedron::getFunctionSpace(int o) const
return 0;
}
-int MTetrahedronN::getNumEdgesRep(){ return 6 * numSubEdges; }
+int MTetrahedronN::getNumEdgesRep(){ return 6 * NUM_SUB_EDGES; }
void MTetrahedronN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
+ int numSubEdges = NUM_SUB_EDGES;
+
static double pp[4][3] = {{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
static int ed [6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
int iEdge = num / numSubEdges;
int iSubEdge = num % numSubEdges;
-
int iVertex1 = ed [iEdge][0];
int iVertex2 = ed [iEdge][1];
- double t1 = (double) iSubEdge / (double) numSubEdges;
- double u1 = pp[iVertex1][0] * (1.-t1) + pp[iVertex2][0] * t1;
- double v1 = pp[iVertex1][1] * (1.-t1) + pp[iVertex2][1] * t1;
- double w1 = pp[iVertex1][2] * (1.-t1) + pp[iVertex2][2] * t1;
+ double t1 = (double)iSubEdge / (double)numSubEdges;
+ double u1 = pp[iVertex1][0] * (1. - t1) + pp[iVertex2][0] * t1;
+ double v1 = pp[iVertex1][1] * (1. - t1) + pp[iVertex2][1] * t1;
+ double w1 = pp[iVertex1][2] * (1. - t1) + pp[iVertex2][2] * t1;
- double t2 = (double) (iSubEdge+1) / (double) numSubEdges;
- double u2 = pp[iVertex1][0] * (1.-t2) + pp[iVertex2][0] * t2;
- double v2 = pp[iVertex1][1] * (1.-t2) + pp[iVertex2][1] * t2;
- double w2 = pp[iVertex1][2] * (1.-t2) + pp[iVertex2][2] * t2;
+ double t2 = (double)(iSubEdge + 1) / (double)numSubEdges;
+ double u2 = pp[iVertex1][0] * (1. - t2) + pp[iVertex2][0] * t2;
+ double v2 = pp[iVertex1][1] * (1. - t2) + pp[iVertex2][1] * t2;
+ double w2 = pp[iVertex1][2] * (1. - t2) + pp[iVertex2][2] * t2;
SPoint3 pnt1, pnt2;
- pnt(u1,v1,w1,pnt1);
- pnt(u2,v2,w2,pnt2);
+ pnt(u1, v1, w1, pnt1);
+ pnt(u2, v2, w2, 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();
}
-int MTetrahedronN::getNumFacesRep(){ return 4 * numSubEdges * numSubEdges ; }
+int MTetrahedronN::getNumFacesRep(){ return 4 * NUM_SUB_EDGES * NUM_SUB_EDGES; }
void MTetrahedronN::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
{
+ int numSubEdges = NUM_SUB_EDGES;
+
static double pp[4][3] = {{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
static int fak [4][3] = {{0,1,2},{0,1,3},{0,2,3},{1,2,3}};
int iFace = num / (numSubEdges * numSubEdges);
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 18515fba75..0269e59a53 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -250,6 +250,7 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
if (p1.size() == 1 && p2.size() == 1){
param1 = p1[0];
param2 = p2[0];
+ return true;
}
else if (p1.size() == 1 && p2.size() == 2){
double d1 =
@@ -260,6 +261,7 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
(p1[0].x() - p2[1].y()) * (p1[0].y() - p2[1].y());
param1 = p1[0];
param2 = d2 < d1 ? p2[1] : p2[0];
+ return true;
}
else if (p2.size() == 1 && p1.size() == 2){
double d1 =
@@ -270,12 +272,20 @@ bool reparamMeshEdgeOnFace(MVertex *v1, MVertex *v2, GFace *gf,
(p2[0].x() - p1[1].y()) * (p2[0].y() - p1[1].y());
param1 = d2 < d1 ? p1[1] : p1[0];
param2 = p2[0];
+ return true;
}
+ else if(p1.size() > 1 && p2.size() > 1){
+ param1 = p1[0];
+ param2 = p2[0];
+ // shout, both vertices are on seams
+ return false;
+ }
else{
+ // brute force!
param1 = gf->parFromPoint(SPoint3(v1->x(), v1->y(), v1->z()));
param2 = gf->parFromPoint(SPoint3(v2->x(), v2->y(), v2->z()));
+ return true;
}
- return true;
}
bool reparamMeshVertexOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
@@ -296,7 +306,7 @@ bool reparamMeshVertexOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
GVertex *gv = (GVertex*)v->onWhat();
param = gv->reparamOnFace(gf, 1);
- // shout if we could be on a seam
+ // shout, we could be on a seam
std::list<GEdge*> ed = gv->edges();
for(std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); it++)
if((*it)->isSeam(gf)) return false;
@@ -307,16 +317,19 @@ bool reparamMeshVertexOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
v->getParameter(0, t);
param = ge->reparamOnFace(gf, t, 1);
- // shout if we are on a seam
+ // shout, we are on a seam
if(ge->isSeam(gf))
return false;
}
else{
double uu, vv;
- if(v->onWhat() == gf && v->getParameter(0, uu) && v->getParameter(1, vv))
+ if(v->onWhat() == gf && v->getParameter(0, uu) && v->getParameter(1, vv)){
param = SPoint2(uu, vv);
- else
+ }
+ else {
+ // brute force!
param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+ }
}
return true;
}
--
GitLab