diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 40aaaeb53e35885c0677a5d3a4aa043f05fe9ee8..29fb4b56cf243ba16c4b9727b2f60c6449e34d94 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -27,7 +27,6 @@ GEdgeCompound::GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound)
int N = _compound.size();
v0 = _orientation[0] ? _compound[0]->getBeginVertex() : _compound[0]->getEndVertex();
v1 = _orientation[N-1] ? _compound[N-1]->getEndVertex() : _compound[N-1]->getBeginVertex();
- //printf("vertices of compound %d -> %d\n",v0->tag(),v1->tag());
v0->addEdge(this);
v1->addEdge(this);
@@ -40,7 +39,6 @@ void GEdgeCompound::orderEdges()
std::list<GEdge*> edges ;
for (unsigned int i=0;i<_compound.size();i++){
- //printf("set compound %d for edge %d \n", tag(), _compound[i]->tag());
_compound[i]->setCompound(this);
edges.push_back(_compound[i]);
}
@@ -50,16 +48,13 @@ void GEdgeCompound::orderEdges()
for (std::list<GEdge*>::iterator it = edges.begin() ; it != edges.end() ; ++it){
GVertex *v1 = (*it)->getBeginVertex();
GVertex *v2 = (*it)->getEndVertex();
- //printf("BEFORE ORDERING segment vert=%d, %d \n", v1->tag(), v2->tag());
std::map<GVertex*,GEdge*>::iterator it1 = tempv.find(v1);
if (it1==tempv.end()) {
- //printf("insert v1=%d \n", v1->tag());
tempv.insert(std::make_pair(v1,*it));
}
else tempv.erase(it1);
std::map<GVertex*,GEdge*>::iterator it2 = tempv.find(v2);
if (it2==tempv.end()){
- //printf("insert v2=%d \n", v2->tag());
tempv.insert(std::make_pair(v2,*it));
}
else tempv.erase(it2);
@@ -97,9 +92,12 @@ void GEdgeCompound::orderEdges()
bool found = false;
for (std::list<GEdge*>::iterator it = edges.begin() ; it != edges.end() ; ++it){
GEdge *e = *it;
+ std::list<GEdge*>::iterator itp;
if (e->getBeginVertex() == last){
_c.push_back(e);
- edges.erase(it);
+ itp = it;
+ it++;
+ edges.erase(itp);
_orientation.push_back(1);
last = e->getEndVertex();
found = true;
@@ -107,7 +105,9 @@ void GEdgeCompound::orderEdges()
}
else if (e->getEndVertex() == last){
_c.push_back(e);
- edges.erase(it);
+ itp = it;
+ it++;
+ edges.erase(itp);
_orientation.push_back(0);
last = e->getBeginVertex();
found = true;
@@ -135,21 +135,11 @@ void GEdgeCompound::orderEdges()
if (_compound.size() < 2)return;
if (_orientation[0] && _compound[0]->getEndVertex() != _compound[1]->getEndVertex()
&& _compound[0]->getEndVertex() != _compound[1]->getBeginVertex()){
- // printf("coucou again\n");
for (unsigned int i = 0; i < _compound.size(); i++){
_orientation[i] = !_orientation[i] ;
}
}
-// for (int i=0;i<_compound.size();i++){
-// printf("i = %d , o %d e %d (%d,%d)\n",
-// i, (int)_orientation[i],
-// _compound[i]->tag(),
-// _compound[i]->getBeginVertex()->tag(),
-// _compound[i]->getEndVertex()->tag());
-// }
-// exit(1);
-
return;
}
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 2f8e8cb6880c5b53061d979bf64dfce618b32c09..20a69364cfdc3f33c733f4e6207014e5558f5c2a 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -16,98 +16,6 @@
#include "FunctionSpace.h"
#include "GmshDefines.h"
-struct GradNormals{
-
- gmshVector<double> g1_U;
- gmshVector<double> g1_V;
- gmshVector<double> g2_U;
- gmshVector<double> g2_V;
- gmshVector<double> g3_U;
- gmshVector<double> g3_V;
- SVector3 n1, n2, n3;
- SPoint3 v1, v2, v3;
-
-};
-
-static void normalS(gmshVector<double> &x, gmshVector<double> &res, void *data)
-{
-
- GradNormals *gn = (GradNormals*)data;
-
- //printf("g1_U= %g %g %g %g %g %g \n", gn->g1_U(0), gn->g1_U(1), gn->g1_U(2), gn->g1_U(3), gn->g1_U(4), gn->g1_U(5));
-
- double Xi[6][3];
- Xi[0][0] = gn->v1.x(); Xi[0][1] = gn->v1.y(); Xi[0][2] = gn->v1.z();
- Xi[1][0] = gn->v2.x(); Xi[1][1] = gn->v2.y(); Xi[1][2] = gn->v2.z();
- Xi[2][0] = gn->v3.x(); Xi[2][1] = gn->v3.y(); Xi[2][2] = gn->v3.z();
- Xi[3][0] = x(0); Xi[3][1] = x(1); Xi[3][2] = x(2);
- Xi[4][0] = x(3); Xi[4][1] = x(4); Xi[4][2] = x(5);
- Xi[5][0] = x(6); Xi[5][1] = x(7); Xi[5][2] = x(8);
-
- double dXdU_p1[3], dXdV_p1[3];
- double dXdU_p2[3], dXdV_p2[3];
- double dXdU_p3[3], dXdV_p3[3];
-
- dXdU_p1[0] = 0.; dXdU_p1[1] = 0.;dXdU_p1[2] = 0.;
- dXdU_p2[0] = 0.; dXdU_p2[1] = 0.;dXdU_p2[2] = 0.;
- dXdU_p3[0] = 0.; dXdU_p3[1] = 0.;dXdU_p3[2] = 0.;
- dXdV_p1[0] = 0.; dXdV_p1[1] = 0.;dXdV_p1[2] = 0.;
- dXdV_p2[0] = 0.; dXdV_p2[1] = 0.;dXdV_p2[2] = 0.;
- dXdV_p3[0] = 0.; dXdV_p3[1] = 0.;dXdV_p3[2] = 0.;
- for (int i=0; i< 6; i++){
- for ( int j =0; j<3; j++){
- dXdU_p1[j] += Xi[i][j]*gn->g1_U(i);
- dXdV_p1[j] += Xi[i][j]*gn->g1_V(i);
- dXdU_p2[j] += Xi[i][j]*gn->g2_U(i);
- dXdV_p2[j] += Xi[i][j]*gn->g2_V(i);
- dXdU_p3[j] += Xi[i][j]*gn->g3_U(i);
- dXdV_p3[j] += Xi[i][j]*gn->g3_V(i);
- }
- }
-
- SVector3 dXdU_1(dXdU_p1[0], dXdU_p1[1], dXdU_p1[2]);
- SVector3 dXdV_1(dXdV_p1[0], dXdV_p1[1], dXdV_p1[2]);
- SVector3 dXdU_2(dXdU_p2[0], dXdU_p2[1], dXdU_p2[2]);
- SVector3 dXdV_2(dXdV_p2[0], dXdV_p2[1], dXdV_p2[2]);
- SVector3 dXdU_3(dXdU_p3[0], dXdU_p3[1], dXdU_p3[2]);
- SVector3 dXdV_3(dXdV_p3[0], dXdV_p3[1], dXdV_p3[2]);
-
- SVector3 n1_UV = crossprod(dXdU_1, dXdV_1); n1_UV.normalize();
- SVector3 n2_UV = crossprod(dXdU_2, dXdV_2); n2_UV.normalize();
- SVector3 n3_UV = crossprod(dXdU_3, dXdV_3); n3_UV.normalize();
-
- res(0) = n1_UV.x() - gn->n1.x();
- res(1) = n1_UV.y() - gn->n1.y();
- res(2) = n1_UV.z() - gn->n1.z();
-
- res(3) = n2_UV.x() - gn->n2.x();
- res(4) = n2_UV.y() - gn->n2.y();
- res(5) = n2_UV.z() - gn->n2.z();
-
- res(6) = n3_UV.x() - gn->n3.x();
- res(7) = n3_UV.y() - gn->n3.y();
- res(8) = n3_UV.z() - gn->n3.z();
-
- printf("X4 =%g %g %g \n", x(0), x(1), x(2));
- printf("X5 =%g %g %g \n", x(3), x(4), x(5));
- printf("X6 =%g %g %g \n", x(6), x(7), x(8));
-
- printf(" X n1UV =%g n1 =%g \n", n1_UV.x() , gn->n1.x());
- printf(" Y n1UV =%g n1 =%g \n", n1_UV.y() , gn->n1.y());
- printf(" Z n1UV =%g n1 =%g \n", n1_UV.z() , gn->n1.z());
-
- printf(" X n2UV =%g n2 =%g \n", n2_UV.x() , gn->n2.x());
- printf(" Y n2UV =%g n2 =%g \n", n2_UV.y() , gn->n2.y());
- printf(" Z n2UV =%g n2 =%g \n", n2_UV.z() , gn->n2.z());
-
- printf(" X n3UV =%g n3 =%g \n", n3_UV.x() , gn->n3.x());
- printf(" Y n3UV =%g n3 =%g \n", n3_UV.y() , gn->n3.y());
- printf(" Z n3UV =%g n3 =%g \n", n3_UV.z() , gn->n3.z());
-
-// printf(" *** res =%g %g %g %g %g %g %g %g %g \n", res(0) , res(1) , res(2) , res(3) , res(4) , res(5) , res(6) , res(7) , res(8) );
-
-}
-
class gmshGradientBasedDiffusivity : public gmshFunction<double>
{
private:
@@ -247,7 +155,7 @@ void GFaceCompound::getBoundingEdges()
std::set<GEdge*>::iterator it = _unique.begin();
GVertex *vB = (*it)->getBeginVertex();
GVertex *vE = (*it)->getEndVertex();
- printf("boundary add edge=%d \n", (*it)->tag());
+ //printf("boundary add edge=%d \n", (*it)->tag());
_loop.push_back(*it);
_unique.erase(it);
it++;
@@ -261,7 +169,7 @@ void GFaceCompound::getBoundingEdges()
std::set<GEdge*>::iterator itp;
if ( v1 == vE ){
- printf("boundary add edge=%d \n", (*it)->tag());
+ //printf("boundary add edge=%d \n", (*it)->tag());
_loop.push_back(*it);
itp = it;
it++;
@@ -270,7 +178,7 @@ void GFaceCompound::getBoundingEdges()
i = -1;
}
else if ( v2 == vE){
- printf("boundary add edge=%d \n", (*it)->tag());
+ //printf("boundary add edge=%d \n", (*it)->tag());
_loop.push_back(*it);
itp = it;
it++;
@@ -278,6 +186,9 @@ void GFaceCompound::getBoundingEdges()
vE = v1;
i=-1;
}
+
+ if (it == _unique.end()) break;
+
}
if (vB == vE) {
@@ -695,7 +606,6 @@ GPoint GFaceCompound::point(double par1, double par2) const
getTriangle (par1, par2, <, U,V);
SPoint3 p(3, 3, 0);
if (!lt){
- // Msg::Warning("Re-Parametrized face %d --> point (%g %g) lies outside the domain", tag(), par1,par2);
GPoint gp (p.x(),p.y(),p.z(),this);
gp.setNoSuccess();
return gp;
@@ -757,145 +667,6 @@ GPoint GFaceCompound::point(double par1, double par2) const
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()) );
-
-// }
-// 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(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());
-// }
-
-// 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);
-
-// }
-
-// gmshVector<double> x(9);
-// const gmshFunctionSpace& fs = gmshFunctionSpaces::find(MSH_TRI_6);
-// double g1[6][3];
-// double g2[6][3];
-// double g3[6][3];
-// fs.df(0,0,0,g1);
-// fs.df(1,0,0,g2);
-// fs.df(0,1,0,g3);
-// gmshVector<double> dNdU_1(6);
-// gmshVector<double> dNdV_1(6);
-// for (int i=0; i<6; i++) dNdU_1(i) = g1[i][0];
-// for (int i=0; i<6; i++) dNdV_1(i) = g1[i][1];
-// gmshVector<double> dNdU_2(6);
-// gmshVector<double> dNdV_2(6);
-// for (int i=0; i<6; i++) dNdU_2(i) = g2[i][0];
-// for (int i=0; i<6; i++) dNdV_2(i) = g2[i][1];
-// gmshVector<double> dNdU_3(6);
-// gmshVector<double> dNdV_3(6);
-// for (int i=0; i<6; i++) dNdU_3(i) = g3[i][0];
-// for (int i=0; i<6; i++) dNdV_3(i) = g3[i][1];
-
-
-
-// //x4,y4,z4 //x5,y5,z5 //x6,y6,z6
-// x(0) = .5*(lt->v1.x()+lt->v2.x()); x(1) = .5*(lt->v1.y()+lt->v2.y()); x(2) = .5*(lt->v1.z()+lt->v2.z());
-// x(3) = .5*(lt->v2.x()+lt->v3.x()); x(4) = .5*(lt->v2.y()+lt->v3.y()); x(5) = .5*(lt->v2.z()+lt->v3.z());
-// x(6) = .5*(lt->v1.x()+lt->v3.x()); x(7) = .5*(lt->v1.y()+lt->v3.y()); x(8) = .5*(lt->v1.z()+lt->v3.z());
-
-// GradNormals gn = {dNdU_1, dNdV_1, dNdU_2, dNdV_2,dNdU_3, dNdV_3,
-// n1, n2, n3, lt->v1, lt->v2, lt->v3};
-// printf("AV X4 =%g %g %g \n", x(0), x(1), x(2));
-// printf("AV X5 =%g %g %g \n", x(3), x(4), x(5));
-// printf("AV X6 =%g %g %g \n", x(6), x(7), x(8));
-
-// bool success = newton_fd(normalS, x, &gn);
-// if (success) printf("succees !\n");
-
-// printf("NEW X4 =%g %g %g \n", x(0), x(1), x(2));
-// printf("NEW X5 =%g %g %g \n", x(3), x(4), x(5));
-// printf("NEW X6 =%g %g %g \n", x(6), x(7), x(8));
-// //exit(1);
-
-// SPoint3 X4(x(0), x(1), x(2));
-// SPoint3 X5(x(3), x(4), x(5));
-// SPoint3 X6(x(6), x(7), x(8));
-// double f1[6];
-// fs.f(U,V,0,f1);
-// // printf("U=%g V=%g \n", U, V);
-// // printf("X1 =%g %g %g \n", lt->v1.x(), lt->v1.y(), lt->v1.z());
-// // printf("X2 =%g %g %g \n", lt->v2.x(), lt->v2.y(), lt->v2.z());
-// // printf("X3 =%g %g %g \n", lt->v3.x(), lt->v3.y(), lt->v3.z());
-// // printf("X4 =%g %g %g \n", x(0), x(1), x(2));
-// // printf("X5 =%g %g %g \n", x(3), x(4), x(5));
-// // printf("X6 =%g %g %g \n", x(6), x(7), x(8));
-// // printf("f1=%g %g %g %g %g %g \n", f1[0], f1[1],f1[2],f1[3],f1[4],f1[5]);
-
-
-
-
-
-
-
-
}
@@ -1079,19 +850,19 @@ void GFaceCompound::printStuff() const
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");
+ 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");
diff --git a/Numeric/gmshLaplace.cpp b/Numeric/gmshLaplace.cpp
index 4c227e4a60b88adbcddb5a92d7167bf1d4ac292d..634ac96cb9dcb092f060065c6789442d2520c274 100644
--- a/Numeric/gmshLaplace.cpp
+++ b/Numeric/gmshLaplace.cpp
@@ -37,12 +37,15 @@ void gmshLaplaceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
invjac[2][2] * grads[j][2];
}
+ double pi = 3.14;
double H_x=1.0;
+ double H_y=1.0;
+ double H_z=1.0;
for (int j = 0; j < nbNodes; j++){
for (int k = 0; k <= j; k++){
m(j, k) += (H_x*Grads[j][0] * Grads[k][0] +
- Grads[j][1] * Grads[k][1] +
- Grads[j][2] * Grads[k][2]) * weight * detJ * _diff;
+ H_y*Grads[j][1] * Grads[k][1] +
+ H_z*Grads[j][2] * Grads[k][2]) * weight * detJ * _diff;
}
}
}