Skip to content
Snippets Groups Projects
Commit 7ea76245 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

*** empty log message ***

parent f3109f61
No related branches found
No related tags found
No related merge requests found
......@@ -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();
......
......@@ -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){
......@@ -703,78 +711,121 @@ 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;
t1 = n2 - n1*dot(n1,n2);
t2 = n2*(dot(n1,n2)) - n1;
Pij = lt->v2 - lt->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);
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{
// //quadratic Lagrange mesh
// //-------------------------
// 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;
// t1 = n2 - n1*dot(n1,n2);
// t2 = n2*(dot(n1,n2)) - n1;
// Pij = lt->v2 - lt->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);
// 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()) );
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());
}
// }
// 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;
// t2 = n3 - n2*dot(n2,n3);
// t3 = n3*(dot(n2,n3)) - n2;
// Pij = lt->v3 - lt->v2;
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);
// 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 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());
}
// 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;
// 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);
// 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());
}
// 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(),
......
......@@ -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; }
......
......@@ -1064,9 +1064,10 @@ 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());
//printf("create face2Edges myEdges.size =%d \n", myEdges.size());
//for each actual GEdge
......@@ -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();
}
}
......@@ -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);
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);
mesh_vertices.erase(std::find(mesh_vertices.begin(), mesh_vertices.end(), vE));
}
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 , coord = %g %g %g\n",gvB->tag(),gvB->tag(), gvB->x(), gvB->y(), gvB->z());
v0 = gvB;
v1 = gvB;
//printf("set bound vertices %d %d \n",bound_vertex->tag(),bound_vertex->tag());
v0 = bound_vertex;
v1 = bound_vertex;
}
......@@ -219,7 +236,6 @@ void discreteEdge::setBoundVertices()
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];
std::map<MVertex*, SVector3>::iterator itn = _normals.find(vB);
if (itn == _normals.end()) printf("vB not found \n");
std::map<MVertex*, SVector3>::iterator itn2 = _normals.find(vE);
if (itn2 == _normals.end()) printf("vE not found \n");
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 b300,b030,b003;
SVector3 b210,b120;
double w12,w21;
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);
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;
SPoint3 X3b(.5*(v1.x()+v2.x()), .5*(v1.y()+v2.y()), .5*(v1.z()+v2.z()) );
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 PP(point.x(), point.y(), point.z());
return GPoint(PP.x(),PP.y(),PP.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());
}
// else{
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);
// //quadratic Lagrange mesh
// //-------------------------
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);
// const SVector3 n1 = _normals[vB];
// const SVector3 n2 = _normals[vE];
p = v1*f1[0] + v2*f1[1] + X3*f1[2];
return GPoint(p.x(),p.y(),p.z());
// 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;
// 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 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);
// 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());
// }
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment