From 7ea7624509a6961e2b92e32209ada854154b3d76 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-Fran=C3=A7ois=20Remacle=20=28students=29?=
<jean-francois.remacle@uclouvain.be>
Date: Thu, 18 Jun 2009 12:50:43 +0000
Subject: [PATCH] *** empty log message ***
---
Geo/GEdgeCompound.cpp | 12 +--
Geo/GFaceCompound.cpp | 208 ++++++++++++++++++++++++++++--------------
Geo/GFaceCompound.h | 2 +
Geo/GModel.cpp | 39 ++++----
Geo/discreteEdge.cpp | 171 +++++++++++++++++++++-------------
5 files changed, 276 insertions(+), 156 deletions(-)
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index dddca595df..b4c0c0c6a5 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -16,12 +16,12 @@ GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound)
: GEdge(m, tag, 0 , 0), _compound(compound)
{
- for (std::vector<GEdge*>::iterator it = compound.begin(); it != compound.end(); ++it){
- if (!(*it)) {
- Msg::Error("Incorrect edge in compound line %d\n", tag);
- Msg::Exit(1);
- }
- }
+// for (std::vector<GEdge*>::iterator it = compound.begin(); it != compound.end(); ++it){
+// if (!(*it)) {
+// Msg::Error("Incorrect edge in compound line %d\n", tag);
+// Msg::Exit(1);
+// }
+// }
orderEdges ();
int N = _compound.size();
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 778472437c..86c11ebc90 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -618,6 +618,14 @@ void GFaceCompound::parametrize(iterationStep step) const
}
}
+
+void GFaceCompound::computeNormals (std::map<MVertex*,SVector3> &normals) const
+{
+
+ computeNormals ();
+ normals = _normals;
+
+}
void GFaceCompound::computeNormals () const
{
_normals.clear();
@@ -691,7 +699,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
return lt->gf->point(pv.x(),pv.y());
}
- const int LINEARMESH = false;
+ const bool LINEARMESH = false;
if (LINEARMESH){
@@ -702,79 +710,122 @@ GPoint GFaceCompound::point(double par1, double par2) const
}
else{
-
- //quadratic Lagrange mesh
+
+ //curved PN triangle
//-------------------------
-
+
const SVector3 n1 = _normals[lt->tri->getVertex(0)];
const SVector3 n2 = _normals[lt->tri->getVertex(1)];
const SVector3 n3 = _normals[lt->tri->getVertex(2)];
- SVector3 t1, t2, t3, Pij;
+ SVector3 b300,b030,b003;
+ SVector3 b210,b120,b021,b012,b102,b201,E,VV,b111;
+ double w12,w21,w23,w32,w31,w13;
+
+ b300 = lt->v1;
+ b030 = lt->v2;
+ b003 = lt->v3;
+ w12 = dot(lt->v2 - lt->v1, n1);
+ w21 = dot(lt->v1 - lt->v2, n2);
+ w23 = dot(lt->v3 - lt->v2, n2);
+ w32 = dot(lt->v2 - lt->v3, n3);
+ w31 = dot(lt->v1 - lt->v3, n3);
+ w13 = dot(lt->v3 - lt->v1, n1);
+ b210 = (2*lt->v1 + lt->v2 -w12*n1)*0.333;
+ b120 = (2*lt->v2 + lt->v1-w21*n2)*0.333;
+ b021 = (2*lt->v2 + lt->v3-w23*n2)*0.333;
+ b012 = (2*lt->v3 + lt->v2-w32*n3)*0.333;
+ b102 = (2*lt->v3 + lt->v1-w31*n3)*0.333;
+ b201 = (2*lt->v1 + lt->v3-w13*n1)*0.333;
+ E=(b210+b120+b021+b012+b102+b201)*0.16667;
+ VV=(lt->v1+lt->v2+lt->v3)*0.333;
+ b111=E+(E-VV)*0.5;
+
+ double W = 1-U-V;
+ SVector3 point = b300*W*W*W+b030*U*U*U+b003*V*V*V+
+ b210*3.*W*W*U+b120*3.*W*U*U+b201*3.*W*W*V+
+ b021*3.*U*U*V+b102*3.*W*V*V+b012*3.*U*V*V+
+ b111*6.*W*U*V;
+
+ SPoint3 PP(point.x(), point.y(), point.z());
+ return GPoint(PP.x(),PP.y(),PP.z(),this,par);
+
+ }
+ // else{
- t1 = n2 - n1*dot(n1,n2);
- t2 = n2*(dot(n1,n2)) - n1;
- Pij = lt->v2 - lt->v1;
+// //quadratic Lagrange mesh
+// //-------------------------
- double a = dot(t1,t2)/dot(t1,t1);
- double b = dot(Pij,t1)/dot(t1,t2);
- double ap = dot(t1,t2)/dot(t2,t2);
- double bp = dot(Pij,t2)/dot(t1,t2);
+// const SVector3 n1 = _normals[lt->tri->getVertex(0)];
+// const SVector3 n2 = _normals[lt->tri->getVertex(1)];
+// const SVector3 n3 = _normals[lt->tri->getVertex(2)];
- SPoint3 X4;
- if (dot(n1,n2)/(norm(n1)*norm(n2)) > 0.9999){
- X4.setPosition(.5*(lt->v1.x()+lt->v2.x()), .5*(lt->v1.y()+lt->v2.y()), .5*(lt->v1.z()+lt->v2.z()) );
-
- }
- else{
- SVector3 XX4 = (.5*((ap*bp*t2)-(a*bp*t1)) ) *.25 + (Pij + .5*((a*b*t1) - (ap*bp*t2)))*.5 + lt->v1;
- X4.setPosition(XX4.x(), XX4.y(), XX4.z());
- }
-
- t2 = n3 - n2*dot(n2,n3);
- t3 = n3*(dot(n2,n3)) - n2;
- Pij = lt->v3 - lt->v2;
+// SVector3 t1, t2, t3, Pij;
- a = dot(t2,t3)/dot(t2,t2);
- b = dot(Pij,t2)/dot(t2,t3);
- ap = dot(t2,t3)/dot(t3,t3);
- bp = dot(Pij,t3)/dot(t2,t3);
+// t1 = n2 - n1*dot(n1,n2);
+// t2 = n2*(dot(n1,n2)) - n1;
+// Pij = lt->v2 - lt->v1;
- SPoint3 X5;
- if (dot(n2,n3)/(norm(n2)*norm(n3)) > 0.9999){
- X5.setPosition(.5*(lt->v2.x()+lt->v3.x()), .5*(lt->v2.y()+lt->v3.y()), .5*(lt->v2.z()+lt->v3.z()));
- }
- else{
- SVector3 XX5 = (.5*((ap*bp*t3)-(a*bp*t2)) ) *.35 + (Pij + .5*((a*b*t2) - (ap*bp*t3)))*.5 + lt->v2;
- X5.setPosition(XX5.x(), XX5.y(), XX5.z());
- }
+// double a = dot(t1,t2)/dot(t1,t1);
+// double b = dot(Pij,t1)/dot(t1,t2);
+// double ap = dot(t1,t2)/dot(t2,t2);
+// double bp = dot(Pij,t2)/dot(t1,t2);
- t3 = n1 - n3*dot(n3,n1);
- t1 = n1*(dot(n3,n1)) - n3;
- Pij = lt->v1 - lt->v3;
+// SPoint3 X4;
+// if (dot(n1,n2)/(norm(n1)*norm(n2)) > 0.9999){
+// printf("close zero \n");
+// X4.setPosition(.5*(lt->v1.x()+lt->v2.x()), .5*(lt->v1.y()+lt->v2.y()), .5*(lt->v1.z()+lt->v2.z()) );
+
+// }
+// else{
+// SVector3 XX4 = (.5*((ap*bp*t2)-(a*bp*t1)) ) *.25 + (Pij + .5*((a*b*t1) - (ap*bp*t2))) *0.5 + lt->v1;
+// X4.setPosition(XX4.x(), XX4.y(), XX4.z());
+// }
+
+// t2 = n3 - n2*dot(n2,n3);
+// t3 = n3*(dot(n2,n3)) - n2;
+// Pij = lt->v3 - lt->v2;
- a = dot(t3,t1)/dot(t3,t3);
- b = dot(Pij,t3)/dot(t3,t1);
- ap = dot(t3,t1)/dot(t1,t1);
- bp = dot(Pij,t1)/dot(t3,t1);
+// a = dot(t2,t3)/dot(t2,t2);
+// b = dot(Pij,t2)/dot(t2,t3);
+// ap = dot(t2,t3)/dot(t3,t3);
+// bp = dot(Pij,t3)/dot(t2,t3);
- SPoint3 X6;
- if (dot(n1,n3)/(norm(n1)*norm(n3)) > 0.9999){
- X6.setPosition(.5*(lt->v1.x()+lt->v3.x()), .5*(lt->v1.y()+lt->v3.y()), .5*(lt->v1.z()+lt->v3.z()) );
- }
- else{
- SVector3 XX6 = (.5*((ap*bp*t1)-(a*bp*t3)) ) *.15 + (Pij + .5*((a*b*t3) - (ap*bp*t1)))*.5 + lt->v3;
- X6.setPosition(XX6.x(), XX6.y(), XX6.z());
- }
+// SPoint3 X5;
+// if (dot(n2,n3)/(norm(n2)*norm(n3)) > 0.9999){
+// X5.setPosition(.5*(lt->v2.x()+lt->v3.x()), .5*(lt->v2.y()+lt->v3.y()), .5*(lt->v2.z()+lt->v3.z()));
+// }
+// else{
+// SVector3 XX5 = (.5*((ap*bp*t3)-(a*bp*t2)) ) *.35 + (Pij + .5*((a*b*t2) - (ap*bp*t3)))*.5 + lt->v2;
+// X5.setPosition(XX5.x(), XX5.y(), XX5.z());
+// }
+
+// t3 = n1 - n3*dot(n3,n1);
+// t1 = n1*(dot(n3,n1)) - n3;
+// Pij = lt->v1 - lt->v3;
+
+// a = dot(t3,t1)/dot(t3,t3);
+// b = dot(Pij,t3)/dot(t3,t1);
+// ap = dot(t3,t1)/dot(t1,t1);
+// bp = dot(Pij,t1)/dot(t3,t1);
+
+// SPoint3 X6;
+// if (dot(n1,n3)/(norm(n1)*norm(n3)) > 0.9999){
+// X6.setPosition(.5*(lt->v1.x()+lt->v3.x()), .5*(lt->v1.y()+lt->v3.y()), .5*(lt->v1.z()+lt->v3.z()) );
+// }
+// else{
+// SVector3 XX6 = (.5*((ap*bp*t1)-(a*bp*t3)) ) *.15 + (Pij + .5*((a*b*t3) - (ap*bp*t1)))*.5 + lt->v3;
+// X6.setPosition(XX6.x(), XX6.y(), XX6.z());
+// }
- const gmshFunctionSpace& fs = gmshFunctionSpaces::find(MSH_TRI_6);
- double f1[6];
- fs.f(U,V,0,f1);
- p = lt->v1*f1[0] + lt->v2*f1[1] + lt->v3*f1[2] + X4*f1[3] + X5*f1[4] + X6*f1[5];
- return GPoint(p.x(),p.y(),p.z(),this,par);
+// const gmshFunctionSpace& fs = gmshFunctionSpaces::find(MSH_TRI_6);
+// double f1[6];
+// fs.f(U,V,0,f1);
+// p = lt->v1*f1[0] + lt->v2*f1[1] + lt->v3*f1[2] + X4*f1[3] + X5*f1[4] + X6*f1[5];
+// return GPoint(p.x(),p.y(),p.z(),this,par);
- }
+// }
// gmshVector<double> x(9);
// const gmshFunctionSpace& fs = gmshFunctionSpaces::find(MSH_TRI_6);
@@ -1020,12 +1071,28 @@ void GFaceCompound::printStuff() const
{
std::list<GFace*> :: const_iterator it = _compound.begin();
- FILE * uvx = fopen("UVX.pos","w");
- FILE * uvy = fopen("UVY.pos","w");
- FILE * uvz = fopen("UVZ.pos","w");
- FILE * xyzu = fopen("XYZU.pos","w");
- FILE * xyzv = fopen("XYZV.pos","w");
- FILE * xyzc = fopen("XYZC.pos","w");
+ char name1[256], name2[256], name3[256];
+ char name4[256], name5[256], name6[256];
+ //sprintf(name1, "UVX-%d.pos", (*it)->tag());
+ //sprintf(name2, "UVY-%d.pos", (*it)->tag());
+ //sprintf(name3, "UVZ-%d.pos", (*it)->tag());
+ //sprintf(name4, "XYZU-%d.pos", (*it)->tag());
+ //sprintf(name5, "XYZV-%d.pos", (*it)->tag());
+ //sprintf(name6, "XYZC-%d.pos", (*it)->tag());
+
+ sprintf(name1, "UVX.pos");
+ sprintf(name2, "UVY.pos");
+ sprintf(name3, "UVZ.pos");
+ sprintf(name4, "XYZU.pos");
+ sprintf(name5, "XYZV.pos");
+ sprintf(name6, "XYZC.pos");
+
+ FILE * uvx = fopen(name1,"w");
+ FILE * uvy = fopen(name2,"w");
+ FILE * uvz = fopen(name3,"w");
+ FILE * xyzu = fopen(name4,"w");
+ FILE * xyzv = fopen(name5,"w");
+ FILE * xyzc = fopen(name6,"w");
fprintf(uvx,"View \"\"{\n");
fprintf(uvy,"View \"\"{\n");
@@ -1043,13 +1110,18 @@ void GFaceCompound::printStuff() const
coordinates.find(t->getVertex(1));
std::map<MVertex*,SPoint3>::const_iterator it2 =
coordinates.find(t->getVertex(2));
- fprintf(xyzu,"VT(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
+ fprintf(xyzu,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
- (35.*it0->second.x()-t->getVertex(0)->x()), -t->getVertex(0)->y(), (35.*it0->second.y()-t->getVertex(0)->z()),
- (35.*it1->second.x()-t->getVertex(1)->x()), -t->getVertex(1)->y(), (35.*it1->second.y()-t->getVertex(1)->z()),
- (35.*it2->second.x()-t->getVertex(2)->x()), -t->getVertex(2)->y(), (35.*it2->second.y()-t->getVertex(2)->z()));
+ it0->second.x(),it1->second.x(),it2->second.x());
+// fprintf(xyzu,"VT(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
+// t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
+// t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
+// t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
+// (35.*it0->second.x()-t->getVertex(0)->x()), -t->getVertex(0)->y(), (35.*it0->second.y()-t->getVertex(0)->z()),
+// (35.*it1->second.x()-t->getVertex(1)->x()), -t->getVertex(1)->y(), (35.*it1->second.y()-t->getVertex(1)->z()),
+// (35.*it2->second.x()-t->getVertex(2)->x()), -t->getVertex(2)->y(), (35.*it2->second.y()-t->getVertex(2)->z()));
fprintf(xyzv,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 7bfa33c257..f1a17fc106 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -45,6 +45,7 @@ class Octree;
class GFaceCompound : public GFace {
public:
typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
+ void computeNormals (std::map<MVertex*,SVector3> &normals) const;
protected:
std::list<GFace*> _compound;
std::list<GEdge*> _U0, _U1, _V0, _V1;
@@ -79,6 +80,7 @@ public:
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return 0; }
+
virtual SPoint2 getCoordinates(MVertex *v) const;
virtual bool buildRepresentationCross(){ return false; }
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 2e12a708b8..cf2a6afd60 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1033,7 +1033,7 @@ void GModel::createTopologyFromMesh()
(*edge)->getEndVertex()->mesh_vertices[0] == myEdges[0].getVertex(1) ) ||
( (*edge)->getBeginVertex()->mesh_vertices[0] == myEdges[0].getVertex(1) &&
(*edge)->getEndVertex()->mesh_vertices[0] == myEdges[0].getVertex(0) )){
- // printf("**********add tagedge =%d \n", (*edge)->tag());
+ //printf("**********add tagedge =%d \n", (*edge)->tag());
tagEdges.push_back((*edge)->tag());
}
}
@@ -1064,12 +1064,13 @@ void GModel::createTopologyFromMesh()
}
face2Edges.insert(std::make_pair(*itFace, tagEdges));
}
+ //printf(" !!! END discrete edges already exist %d \n", myEdges.size());
}
else{
- printf("create face2Edges myEdges.size =%d \n", myEdges.size());
-
-
- //for each actual GEdge
+ //printf("create face2Edges myEdges.size =%d \n", myEdges.size());
+
+
+ //for each actual GEdge
while (! myEdges.empty()) {
std::vector<MEdge> myLines;
myLines.clear();
@@ -1088,17 +1089,17 @@ void GModel::createTopologyFromMesh()
for (std::vector<MEdge>::iterator it = myEdges.begin() ; it != myEdges.end(); it++){
MVertex *v1 = (*it).getVertex(0);
MVertex *v2 = (*it).getVertex(1);
- printf("mline %d %d size=%d\n", v1->getNum(), v2->getNum(), myEdges.size());
+ //printf("mline %d %d size=%d\n", v1->getNum(), v2->getNum(), myEdges.size());
if ( v1 == vE ){
- printf("->v1 = vE push back this mline \n");
+ //printf("->v1 = vE push back this mline \n");
myLines.push_back(*it);
myEdges.erase(it);
vE = v2;
i = -1;
}
else if ( v2 == vE){
- printf("->v2 = VE push back this mline \n");
+ //printf("->v2 = VE push back this mline \n");
myLines.push_back(*it);
myEdges.erase(it);
vE = v1;
@@ -1110,26 +1111,26 @@ void GModel::createTopologyFromMesh()
printf("end Edges \n");
if (vB == vE) {
- printf("vB = ve = \n");
+ //printf("vB = ve = \n");
break;
}
if (myEdges.empty()) break;
- printf("not found VB=%d vE=%d\n", vB->getNum(), vE->getNum());
+ //printf("not found VB=%d vE=%d\n", vB->getNum(), vE->getNum());
MVertex *temp = vB;
vB = vE;
vE = temp;
- printf("not found VB=%d vE=%d\n", vB->getNum(), vE->getNum());
+ //printf("not found VB=%d vE=%d\n", vB->getNum(), vE->getNum());
}
- printf("************ CANDIDATE NEW EDGE with num =%d\n", num);
- for (std::vector<MEdge>::iterator it = myLines.begin() ; it != myLines.end() ; ++it){
- MVertex *v1 = (*it).getVertex(0);
- MVertex *v2 = (*it).getVertex(1);
- printf("Line %d %d \n", v1->getNum(), v2->getNum());
- }
+// printf("************ CANDIDATE NEW EDGE with num =%d\n", num);
+// for (std::vector<MEdge>::iterator it = myLines.begin() ; it != myLines.end() ; ++it){
+// MVertex *v1 = (*it).getVertex(0);
+// MVertex *v2 = (*it).getVertex(1);
+// printf("Line %d %d \n", v1->getNum(), v2->getNum());
+// }
discreteEdge *e = new discreteEdge(this, num, 0, 0);
add(e);
Dedges.push_back(e);
@@ -1146,9 +1147,7 @@ void GModel::createTopologyFromMesh()
for (std::vector<int>::iterator itFace = tagFaces.begin(); itFace != tagFaces.end(); itFace++) {
GFace *dFace = getFaceByTag(abs(*itFace));
- printf("face =%d \n", dFace->tag());
for (std::list<MVertex*>::iterator itv = all_vertices.begin(); itv != all_vertices.end(); itv++) {
- printf("vertrx=%d \n", (*itv)->getNum());
std::vector<MVertex*>::iterator itve = std::find(dFace->mesh_vertices.begin(), dFace->mesh_vertices.end(), *itv) ;
if (itve != dFace->mesh_vertices.end()) dFace->mesh_vertices.erase(itve);
(*itv)->setEntity(e);
@@ -1232,5 +1231,7 @@ void GModel::createTopologyFromMesh()
(*edge)->parametrize();
}
+
+
}
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 68e16ac4f2..361ca618ab 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -16,6 +16,7 @@
#include "MHexahedron.h"
#include "MPyramid.h"
+
#include <vector>
#include <list>
@@ -157,13 +158,13 @@ void discreteEdge::setBoundVertices()
{
if (boundv.size()==2){
- //printf("Found a regular open Curve \n");
+ printf("Found a regular open Curve \n");
std::vector<GVertex*> bound_vertices;
for (std::map<MVertex*,MLine*>::const_iterator iter = boundv.begin(); iter != boundv.end(); iter++){
MVertex* vE = (iter)->first;
bool existVertex = false;
for(GModel::viter point = model()->firstVertex(); point != model()->lastVertex(); point++){
- //printf("Discrete point =%d %d bound vE=%d %d\n",(*point)->tag(), (*point)->mesh_vertices[0]->getNum(), vE->getNum(), vE->getIndex());
+ //printf("Discrete point =%d %d bound vE=%d\n",(*point)->tag(), (*point)->mesh_vertices[0]->getNum(), vE->getNum());
if ((*point)->mesh_vertices[0]->getNum() == vE->getNum()){
bound_vertices.push_back((*point));
existVertex = true;
@@ -190,18 +191,34 @@ void discreteEdge::setBoundVertices()
}
else if (boundv.size()==0){
//printf("Found a closed Curve add vertex \n");
-
+ GVertex* bound_vertex;
std::vector<MLine*>::const_iterator it = lines.begin();
MVertex* vE = (*it)->getVertex(0);
- GVertex *gvB = new discreteVertex(model(),vE->getNum());
- gvB->mesh_vertices.push_back(vE);
- gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
- model()->add(gvB);
- mesh_vertices.erase(std::find(mesh_vertices.begin(), mesh_vertices.end(), vE));
+ bool existVertex = false;
+
+ for(GModel::viter point = model()->firstVertex(); point != model()->lastVertex(); point++){
+ //printf("Discrete point =%d %d bound vE=%d \n",(*point)->tag(), (*point)->mesh_vertices[0]->getNum(), vE->getNum());
+ if ((*point)->mesh_vertices[0]->getNum() == vE->getNum()){
+ bound_vertex = (*point);
+ existVertex = true;
+ break;
+ }
+ }
+ if(!existVertex){
+
+ GVertex *gvB = new discreteVertex(model(),vE->getNum());
+ bound_vertex = gvB;
+ gvB->mesh_vertices.push_back(vE);
+ gvB->points.push_back(new MPoint(gvB->mesh_vertices.back()));
+ model()->add(gvB);
+ }
- ///printf("set bound vertices %d %d , coord = %g %g %g\n",gvB->tag(),gvB->tag(), gvB->x(), gvB->y(), gvB->z());
- v0 = gvB;
- v1 = gvB;
+ std::vector<MVertex*>::iterator itve = std::find(mesh_vertices.begin(), mesh_vertices.end(), vE) ;
+ if (itve != mesh_vertices.end()) mesh_vertices.erase(itve);
+
+ //printf("set bound vertices %d %d \n",bound_vertex->tag(),bound_vertex->tag());
+ v0 = bound_vertex;
+ v1 = bound_vertex;
}
@@ -217,8 +234,7 @@ void discreteEdge::setBoundVertices()
_pars[0]=0 _pars[1]=1 _pars[2]=2 _pars[N+1]=N+1
*/
void discreteEdge::parametrize()
-{
-
+{
for (int i=0;i<lines.size()+1;i++){
_pars.push_back(i);
@@ -339,11 +355,21 @@ void discreteEdge::parametrize()
// }
+ computeNormals();
+
}
void discreteEdge::computeNormals () const
{
- printf("!!!!!!! dans compute normals l_face size =%d \n", l_faces.size());
+
+// _normals.clear();
+// for(std::list<GFace*>::const_iterator iFace = l_faces.begin(); iFace != l_faces.end(); ++iFace){
+// GFaceCompound* myFace;
+// myFace = (GFaceCompound*) *iFace;
+// myFace->computeNormals(_normals);
+// printf("coucou \n");
+// }
+
_normals.clear();
double J[3][3];
@@ -365,8 +391,8 @@ void discreteEdge::computeNormals () const
std::map<MVertex*,SVector3>::iterator itn = _normals.begin();
for ( ; itn != _normals.end() ; ++itn){
itn->second.normalize();
- printf("NUM = %d xx = %g %g %g \n", itn->first->getNum(), itn->first->x(), itn->first->y(), itn->first->z() ) ;
- printf("normal =%g %g %g \n", itn->second.x(), itn->second.y(), itn->second.z());
+ //printf("NUM = %d xx = %g %g %g \n", itn->first->getNum(), itn->first->x(), itn->first->y(), itn->first->z() ) ;
+ //printf("normal =%g %g %g \n", itn->second.x(), itn->second.y(), itn->second.z());
}
printf("********* normal size = %d \n", _normals.size());
@@ -398,7 +424,7 @@ GPoint discreteEdge::point(double par) const
MVertex *vB = lines[iEdge]->getVertex(0);
MVertex *vE = lines[iEdge]->getVertex(1);
- const int LINEARMESH = false;
+ const bool LINEARMESH = false;
if (LINEARMESH){
@@ -408,68 +434,87 @@ GPoint discreteEdge::point(double par) const
y = vB->y() + tLoc * (vE->y()- vB->y());
z = vB->z() + tLoc * (vE->z()- vB->z());
- //printf("dans point par=%d iEdge =%d, tLoc=%d \n", par, iEdge, tLoc);
- //printf("Discrete Edge POint par=%g, x= %g, y = %g, z = %g \n", par, x,y,z);
-
return GPoint(x,y,z);
}
else{
-
- //quadratic Lagrange mesh
+
+ //curved PN triangle
//-------------------------
- computeNormals();
const SVector3 n1 = _normals[vB];
const SVector3 n2 = _normals[vE];
+
+ SPoint3 v1(vB->x(),vB->y(),vB->z());
+ SPoint3 v2(vE->x(),vE->y(),vE->z());
+
+ SVector3 b300,b030,b003;
+ SVector3 b210,b120;
+ double w12,w21;
- std::map<MVertex*, SVector3>::iterator itn = _normals.find(vB);
- if (itn == _normals.end()) printf("vB not found \n");
+ b300 = v1;
+ b030 = v2;
+ w12 = dot(v2 - v1, n1);
+ w21 = dot(v1 - v2, n2);
+ b210 = (2*v1 + v2 -w12*n1)*0.333;
+ b120 = (2*v2 + v1-w21*n2)*0.333;
- std::map<MVertex*, SVector3>::iterator itn2 = _normals.find(vE);
- if (itn2 == _normals.end()) printf("vE not found \n");
+ double U = tLoc;
+ double W = 1-U;
+ SVector3 point = b300*W*W*W+b030*U*U*U+b210*3.*W*W*U+b120*3.*W*U*U;
- SPoint3 v1(vB->x(),vB->y(),vB->z());
- SPoint3 v2(vE->x(),vE->y(),vE->z());
+ SPoint3 PP(point.x(), point.y(), point.z());
+ return GPoint(PP.x(),PP.y(),PP.z());
+
+ }
+// else{
+
+// //quadratic Lagrange mesh
+// //-------------------------
+
+// const SVector3 n1 = _normals[vB];
+// const SVector3 n2 = _normals[vE];
+
+// SPoint3 v1(vB->x(),vB->y(),vB->z());
+// SPoint3 v2(vE->x(),vE->y(),vE->z());
- SVector3 t1, t2, Pij;
- t1 = n2 - n1*dot(n1,n2);
- t2 = n2*(dot(n1,n2)) - n1;
- Pij = v2 - v1;
+// SVector3 t1, t2, Pij;
+// t1 = n2 - n1*dot(n1,n2);
+// t2 = n2*(dot(n1,n2)) - n1;
+// Pij = v2 - v1;
- double a = dot(t1,t2)/dot(t1,t1);
- double b = dot(Pij,t1)/dot(t1,t2);
- double ap = dot(t1,t2)/dot(t2,t2);
- double bp = dot(Pij,t2)/dot(t1,t2);
+// double a = dot(t1,t2)/dot(t1,t1);
+// double b = dot(Pij,t1)/dot(t1,t2);
+// double ap = dot(t1,t2)/dot(t2,t2);
+// double bp = dot(Pij,t2)/dot(t1,t2);
- SPoint3 X3b(.5*(v1.x()+v2.x()), .5*(v1.y()+v2.y()), .5*(v1.z()+v2.z()) );
-
- SPoint3 X3;
- if (dot(n1,n2)/(norm(n1)*norm(n2)) > 0.9999){
- X3.setPosition(.5*(v1.x()+v2.x()), .5*(v1.y()+v2.y()), .5*(v1.z()+v2.z()) );
- }
- else{
- SVector3 XX3 = (.5*((ap*bp*t2)-(a*bp*t1)) ) *.25 + (Pij + .5*((a*b*t1) - (ap*bp*t2)))*.5 + v1;
- X3.setPosition(XX3.x(), XX3.y(), XX3.z());
- }
+// SPoint3 X3;
+// if (dot(n1,n2)/(norm(n1)*norm(n2)) > 0.9999){
+// printf("coucou normals \n");
+// X3.setPosition(.5*(v1.x()+v2.x()), .5*(v1.y()+v2.y()), .5*(v1.z()+v2.z()) );
+// }
+// else{
+// SVector3 XX3 = (.5*((ap*bp*t2)-(a*bp*t1)) ) *.25 + (Pij + .5*((a*b*t1) - (ap*bp*t2)))*.5 + v1;
+// X3.setPosition(XX3.x(), XX3.y(), XX3.z());
+// }
- printf("normal x1 num=%d coord = %g %g %g \n" , vB->getNum(), v1.x(), v1.y(), v1.z());
- printf("normal n1 = %g %g %g \n", n1.x(), n1.y(), n1.z());
- printf("normal x2 num = %d coord = %g %g %g \n", vE->getNum(), v2.x(), v2.y(), v2.z());
- printf("normal n2 = %g %g %g \n", n2.x(), n2.y(), n2.z());
- printf("normal x3 = %g %g %g \n", X3.x(), X3.y(), X3.z());
- printf("normal x3b = %g %g %g \n", X3b.x(), X3b.y(), X3b.z());
- //exit(1);
+// // printf("normal x1 num=%d coord = %g %g %g \n" , vB->getNum(), v1.x(), v1.y(), v1.z());
+// // printf("normal n1 = %g %g %g \n", n1.x(), n1.y(), n1.z());
+// // printf("normal x2 num = %d coord = %g %g %g \n", vE->getNum(), v2.x(), v2.y(), v2.z());
+// // printf("normal n2 = %g %g %g \n", n2.x(), n2.y(), n2.z());
+// // printf("normal x3 = %g %g %g \n", X3.x(), X3.y(), X3.z());
+// // printf("normal x3b = %g %g %g \n", X3b.x(), X3b.y(), X3b.z());
+// //exit(1);
- const gmshFunctionSpace& fs = gmshFunctionSpaces::find(MSH_LIN_3);
- double f1[3];
- SPoint3 p(0, 0, 0);
- double U = 2*tLoc -1;
- fs.f(U,0,0,f1);
-
- p = v1*f1[0] + v2*f1[1] + X3*f1[2];
- return GPoint(p.x(),p.y(),p.z());
+// const gmshFunctionSpace& fs = gmshFunctionSpaces::find(MSH_LIN_3);
+// double f1[3];
+// SPoint3 p(0, 0, 0);
+// double U = 2*tLoc -1;
+// fs.f(U,0,0,f1);
+
+// p = v1*f1[0] + v2*f1[1] + X3*f1[2];
+// return GPoint(p.x(),p.y(),p.z());
- }
+// }
}
--
GitLab