From c292e8e0084dd9524fb7e45889225523c63ec61d Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Mon, 20 Sep 2010 10:42:36 +0000
Subject: [PATCH] Improved conformal map + first part nonlinear conformal map
---
Common/DefaultOptions.h | 2 +-
Geo/GFaceCompound.cpp | 219 ++++++++++++++++++++--------------------
Geo/GFaceCompound.h | 2 +-
Geo/GModelIO_Mesh.cpp | 4 -
Plugin/Distance.cpp | 1 -
Solver/crossConfTerm.h | 32 +++++-
Solver/eigenSolver.cpp | 22 ++--
Solver/laplaceTerm.h | 32 +++++-
8 files changed, 176 insertions(+), 138 deletions(-)
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 1c7470f07e..fee1a8ccb3 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1253,7 +1253,7 @@ StringXNumber MeshOptions_Number[] = {
{ F|O, "SmoothRatio" , opt_mesh_smooth_ratio , 1.8 ,
"Ratio between mesh sizes at vertices of a same edeg (used in BAMG)" },
{ F|O, "SubdivisionAlgorithm" , opt_mesh_algo_subdivide , 0 ,
- "Mesh subdivision algorithm (0=none, 1=all quadrangles, 2=all hexahedra)" },
+ "Mesh subdivision algorithm (0=none, 1=all quadrangles, 2=all hexahedra, 3=all quadrangles (Blossom))" },
{ F|O, "SurfaceEdges" , opt_mesh_surfaces_edges , 1. ,
"Display edges of surface mesh?" },
{ F|O, "SurfaceFaces" , opt_mesh_surfaces_faces , 0. ,
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index f8b280dc50..7ba5cdc134 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -512,7 +512,7 @@ bool GFaceCompound::checkOverlap(std::vector<MVertex*> &ordered) const
}
if ( !has_no_overlap ) {
- Msg::Warning("$$$ Overlap for compound face %d", this->tag());
+ Msg::Debug("Overlap for compound face %d", this->tag());
}
return has_no_overlap;
@@ -559,7 +559,7 @@ bool GFaceCompound::checkOrientation(int iter) const
int iterMax = 15;
if(!oriented && iter < iterMax){
- if (iter == 0) Msg::Info("--- Parametrization is flipping : applying cavity checks.");
+ if (iter == 0) Msg::Warning("--- Parametrization is flipping : applying cavity checks.");
Msg::Debug("--- Cavity Check - iter %d -",iter);
one2OneMap();
return checkOrientation(iter+1);
@@ -652,13 +652,17 @@ bool GFaceCompound::parametrize() const
Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
bool noOverlap = parametrize_conformal_spectral() ;
- if (!noOverlap) noOverlap = parametrize_conformal();
+ if (!noOverlap){
+ Msg::Warning("!!! Overlap: parametrization switched to 'FE conformal' map");
+ noOverlap = parametrize_conformal();
+ }
// if (!noOverlap) {
- // noOverlap = parametrize_conformal_nonLinear() ;
- // exit(1);
+ // Msg::Warning("!!! Overlap: parametrization switched to 'nonLIN conformal' map");
+ // noOverlap = parametrize_conformal_nonLinear() ;
+ // exit(1);
// }
if ( !noOverlap || !checkOrientation(0) ){
- Msg::Warning("$$$ Parametrization switched to harmonic map");
+ Msg::Warning("$$$ Overlap: parametrization switched to 'harmonic' map");
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
}
@@ -667,7 +671,7 @@ bool GFaceCompound::parametrize() const
buildOct();
if (!checkOrientation(0)){
- Msg::Info("### Parametrization switched to convex combination map");
+ Msg::Info("### Flipping: parametrization switched to convex combination map");
coordinates.clear();
Octree_Delete(oct);
fillNeumannBCS();
@@ -1059,11 +1063,12 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
return SPoint2(0, 0);
}
-void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom, double alpha) const
+void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
{
dofManager<double> myAssembler(_lsys);
+
if(_type == UNITCIRCLE){
// maps the boundary onto a circle
std::vector<MVertex*> ordered;
@@ -1175,21 +1180,13 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom, double al
double value;
myAssembler.getDofValue(v, 0, 1, value);
std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
-
- //combination convex and harmonic
- double valNEW ;
- if (alpha != 0.0)
- valNEW = alpha*itf->second[step] + (1-alpha)*value ;
- else
- valNEW = value;
-
if(itf == coordinates.end()){
SPoint3 p(0, 0, 0);
- p[step] = valNEW;
+ p[step] = value;
coordinates[v] = p;
}
else{
- itf->second[step]= valNEW;
+ itf->second[step]= value;
}
}
@@ -1199,11 +1196,11 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom, double al
bool GFaceCompound::parametrize_conformal_nonLinear() const
{
- Msg::Info("Conformal parametrization of type: nonLinear");
buildOct(); //use this to print conformal map that is overlapping
+ //exit(1);
+ //--create dofManager
dofManager<double> myAssembler(_lsys);
- _lsys->clear();
//--- first compute mapping harmonic
parametrize(ITERU,HARMONIC);
@@ -1222,12 +1219,12 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
int nb = ordered.size();
//--fix vertex
- MVertex *v1 = ordered[0];
- MVertex *v2 = ordered[(int)ceil((double)ordered.size()/2.)];
- myAssembler.fixVertex(v1, 0, 1, 1.);
- myAssembler.fixVertex(v1, 0, 2, 0.);
- myAssembler.fixVertex(v2, 0, 1, -1.);
- myAssembler.fixVertex(v2, 0, 2, 0.);
+ // MVertex *v1 = ordered[0];
+ // MVertex *v2 = ordered[(int)ceil((double)ordered.size()/2.)];
+ // myAssembler.fixVertex(v1, 0, 1, 1.);
+ // myAssembler.fixVertex(v1, 0, 2, 0.);
+ // myAssembler.fixVertex(v2, 0, 1, -1.);
+ // myAssembler.fixVertex(v2, 0, 2, 0.);
//--Assemble linear system
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
@@ -1240,46 +1237,48 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
myAssembler.numberVertex(v, 0, 1);
myAssembler.numberVertex(v, 0, 2);
}
- //ghost vertex for lagrange multiplier
MVertex *lag = new MVertex(0.,0.,0.);
- myAssembler.numberVertex(lag, 0, 3);
+ myAssembler.numberVertex(lag, 0, 3);//ghost vertex for lagrange multiplier
- std::vector<MElement *> allElems;
- laplaceTerm laplace1(model(), 1, ONE, &myAssembler);
- laplaceTerm laplace2(model(), 2, ONE, &myAssembler);
- crossConfTerm cross12(model(), 1, 2, ONE, &myAssembler);
- crossConfTerm cross21(model(), 2, 1, MONE, &myAssembler);
- std::list<GFace*>::const_iterator it = _compound.begin();
- for( ; it != _compound.end() ; ++it){
- for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- SElement se((*it)->triangles[i]);
- laplace1.addToMatrix(myAssembler, &se);
- laplace2.addToMatrix(myAssembler, &se);
- cross12.addToMatrix(myAssembler, &se);
- cross21.addToMatrix(myAssembler, &se);
- allElems.push_back((MElement*)((*it)->triangles[i]));
- }
- }
- for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
- SElement se((*it2));
- laplace1.addToMatrix(myAssembler, &se);
- laplace2.addToMatrix(myAssembler, &se);
- cross12.addToMatrix(myAssembler, &se);
- cross21.addToMatrix(myAssembler, &se);
- allElems.push_back((MElement*)(*it2));
- }
-
- groupOfElements gr(allElems);
- laplace1.addToRightHandSide(myAssembler, gr);
- laplace2.addToRightHandSide(myAssembler, gr);
- cross12.addToRightHandSide(myAssembler, gr);
- cross21.addToRightHandSide(myAssembler, gr);
//--- newton Loop
int nbNewton = 1;
- double lam = 0.1;
+ double lambda = 1.e-3;
for (int iNewton = 0; iNewton < nbNewton; iNewton++){
+ //-- assemble conformal matrix
+ std::vector<MElement *> allElems;
+ laplaceTerm laplace1(model(), 1, ONE, &coordinates);
+ laplaceTerm laplace2(model(), 2, ONE, &coordinates);
+ crossConfTerm cross12(model(), 1, 2, ONE, &coordinates);
+ crossConfTerm cross21(model(), 2, 1, MONE, &coordinates);
+ std::list<GFace*>::const_iterator it = _compound.begin();
+ for( ; it != _compound.end() ; ++it){
+ for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ SElement se((*it)->triangles[i]);
+ laplace1.addToMatrix(myAssembler, &se);
+ laplace2.addToMatrix(myAssembler, &se);
+ cross12.addToMatrix(myAssembler, &se);
+ cross21.addToMatrix(myAssembler, &se);
+ allElems.push_back((MElement*)((*it)->triangles[i]));
+ }
+ }
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ SElement se((*it2));
+ laplace1.addToMatrix(myAssembler, &se);
+ laplace2.addToMatrix(myAssembler, &se);
+ cross12.addToMatrix(myAssembler, &se);
+ cross21.addToMatrix(myAssembler, &se);
+ allElems.push_back((MElement*)(*it2));
+ }
+
+ groupOfElements gr(allElems);
+ laplace1.addToRightHandSide(myAssembler, gr);
+ laplace2.addToRightHandSide(myAssembler, gr);
+ cross12.addToRightHandSide(myAssembler, gr);
+ cross21.addToRightHandSide(myAssembler, gr);
+
+ //-- compute all boundary angles
double sumTheta = 0.0;
for (int i=0; i< nb; i++){
MVertex *prev = (i!=0) ? ordered[i-1] : ordered[nb-1];
@@ -1297,14 +1296,14 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
double a = norm(va), b=norm(vb), c=norm(vc);
SVector3 n = crossprod(va, vb);
- //--- compute theta
+ //- compute theta
double theta = acos((a*a+b*b-c*c)/(2*a*b)); //in rad
double sign = 1.;
if (n.z() < 0.0) {sign = -1.; theta = 2*M_PI-theta;}
//printf("theta in deg =%g \n", theta*180./M_PI);
sumTheta +=theta;
- //--- compute grad_uv theta
+ //- compute grad_uv theta
//dTheta/du1 sign*acos((a^2+b^2 -c^2)/(2*a*b))
//a=sqrt((u2-u1)^2+(v2-v1)^2))
//b=sqrt((u3-u2)^2+(v3-v2)^2))
@@ -1319,34 +1318,34 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
double dTdu3=sign*((u1 + u2 - 2*u3)/(sqrt(SQU(u1 - u2) + SQU(v1 - v2))*sqrt(SQU(u2 - u3) + SQU(v2 - v3))) - ((2*u2 - 2*u3)*(SQU(u1 - u2) + SQU(u1 - u3) + SQU(u2 - u3) + SQU(v1 - v2) + SQU(v1 - v3) + SQU(v2 - v3)))/(4*sqrt(SQU(u1 - u2) + SQU(v1 - v2))*pow((SQU(u2 - u3) + SQU(v2 - v3)),3/2)))/sqrt(1 - SQU(SQU(u1 - u2) + SQU(u1 - u3) + SQU(u2 - u3) + SQU(v1 - v2) + SQU(v1 - v3) + SQU(v2 - v3))/(4*(SQU(u1 - u2) + SQU(v1 - v2))*(SQU(u2 - u3) + SQU(v2 - v3))));
double dTdv3=((v1 + v2 - 2*v3)/(sqrt(SQU(u1 - u2) + SQU(v1 - v2))*sqrt(SQU(u2 - u3) + SQU(v2 - v3))) - ((2*v2 - 2*v3)*(SQU(u1 - u2) + SQU(u1 - u3) + SQU(u2 - u3) + SQU(v1 - v2) + SQU(v1 - v3) + SQU(v2 - v3)))/(4*sqrt(SQU(u1 - u2) + SQU(v1 - v2))*pow((SQU(u2 - u3) + SQU(v2 - v3)),3/2)))/sqrt(1 - SQU(SQU(u1 - u2) + SQU(u1 - u3) + SQU(u2 - u3) + SQU(v1 - v2) + SQU(v1 - v3) + SQU(v2 - v3))/(4*(SQU(u1 - u2) + SQU(v1 - v2))*(SQU(u2 - u3) + SQU(v2 - v3))));
- //---assemble constraint terms
- // myAssembler.assemble(lag, 0, 3, prev, 0, 1, dTdu1);
- // myAssembler.assemble(lag, 0, 3, prev, 0, 2, dTdv1);
- // myAssembler.assemble(lag, 0, 3, curr, 0, 1, dTdu2);
- // myAssembler.assemble(lag, 0, 3, curr, 0, 2, dTdv2);
- // myAssembler.assemble(lag, 0, 3, next, 0, 1, dTdu3);
- // myAssembler.assemble(lag, 0, 3, next, 0, 2, dTdv3);
-
- // myAssembler.assemble(prev, 0, 1, lag, 0, 3, dTdu1);
- // myAssembler.assemble(prev, 0, 2, lag, 0, 3, dTdv1);
- // myAssembler.assemble(curr, 0, 1, lag, 0, 3, dTdu2);
- // myAssembler.assemble(curr, 0, 2, lag, 0, 3, dTdv2);
- // myAssembler.assemble(next, 0, 1, lag, 0, 3, dTdu3);
- // myAssembler.assemble(next, 0, 2, lag, 0, 3, dTdv3);
-
- // myAssembler.assemble(prev, 0, 1, -dTdu1);
- // myAssembler.assemble(prev, 0, 2, -dTdv1);
- // myAssembler.assemble(curr, 0, 1, -dTdu2);
- // myAssembler.assemble(curr, 0, 2, -dTdv2);
- // myAssembler.assemble(next, 0, 1, -dTdu3);
- // myAssembler.assemble(next, 0, 2, -dTdv3);
+ //- assemble constraint terms
+ myAssembler.assemble(lag, 0, 3, prev, 0, 1, lambda*dTdu1);
+ myAssembler.assemble(lag, 0, 3, prev, 0, 2, lambda*dTdv1);
+ myAssembler.assemble(lag, 0, 3, curr, 0, 1, lambda*dTdu2);
+ myAssembler.assemble(lag, 0, 3, curr, 0, 2, lambda*dTdv2);
+ myAssembler.assemble(lag, 0, 3, next, 0, 1, lambda*dTdu3);
+ myAssembler.assemble(lag, 0, 3, next, 0, 2, lambda*dTdv3);
+
+ myAssembler.assemble(prev, 0, 1, lag, 0, 3, lambda*dTdu1);
+ myAssembler.assemble(prev, 0, 2, lag, 0, 3, lambda*dTdv1);
+ myAssembler.assemble(curr, 0, 1, lag, 0, 3, lambda*dTdu2);
+ myAssembler.assemble(curr, 0, 2, lag, 0, 3, lambda*dTdv2);
+ myAssembler.assemble(next, 0, 1, lag, 0, 3, lambda*dTdu3);
+ myAssembler.assemble(next, 0, 2, lag, 0, 3, lambda*dTdv3);
+
+ myAssembler.assemble(prev, 0, 1, lambda*dTdu1);
+ myAssembler.assemble(prev, 0, 2, lambda*dTdv1);
+ myAssembler.assemble(curr, 0, 1, lambda*dTdu2);
+ myAssembler.assemble(curr, 0, 2, lambda*dTdv2);
+ myAssembler.assemble(next, 0, 1, lambda*dTdu3);
+ myAssembler.assemble(next, 0, 2, lambda*dTdv3);
}
//--- compute constraint
double G = fabs(sumTheta - (nb-2)*M_PI);
- printf("constraint G = %g \n", G);
- myAssembler.assemble(lag, 0, 3, lag, 0, 3, 1.0);//change this
+ printf("Sum of angles G = %g \n", G);
+ myAssembler.assemble(lag, 0, 3, lag, 0, 3, 0.0);//change this
//--- assemble rhs of linear system
myAssembler.assemble(lag, 0, 3, G);
@@ -1364,11 +1363,24 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
myAssembler.getDofValue(v, 0, 2, value2);
coordinates[v] = SPoint3(value1,value2,0.0);
}
- double lambda;
- myAssembler.getDofValue(lag, 0, 3, lambda);
- printf("lambda=%g \n", lambda);
+ //-- update newton
+ //U=U+DU
+ //lambda = lambda +dLambda
+ for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ MVertex *v = *itv;
+ double du,dv;
+ myAssembler.getDofValue(v, 0, 1, du);
+ myAssembler.getDofValue(v, 0, 2, dv);
+ SPoint3 old = coordinates[v];
+ coordinates[v] = SPoint3(old.x()+du,old.y()+dv,0.0);
+ }
+ double dLambda;
+ myAssembler.getDofValue(lag, 0, 3, dLambda);
+ lambda = lambda+dLambda;
+ printf("dLambda=%g \n", dLambda);
- //update newton
+ _lsys->zeroMatrix();
+ _lsys->zeroRightHandSide();
}//end Newton
@@ -1378,7 +1390,7 @@ bool GFaceCompound::parametrize_conformal_nonLinear() const
bool GFaceCompound::parametrize_conformal_spectral() const
{
- Msg::Info("Conformal parametrization of type: spectral");
+
#if !defined(HAVE_PETSC) && !defined(HAVE_SLEPC)
{
@@ -1463,36 +1475,20 @@ bool GFaceCompound::parametrize_conformal_spectral() const
//-------------------------------
myAssembler.setCurrentMatrix("B");
- double small = 0.0;
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
- MVertex *v = *itv;
- if (std::find(ordered.begin(), ordered.end(), v) == ordered.end() ){
- myAssembler.assemble(v, 0, 1, v, 0, 1, small);
- myAssembler.assemble(v, 0, 2, v, 0, 2, small);
- }
- else{
- myAssembler.assemble(v, 0, 1, v, 0, 1, 1.0);
- myAssembler.assemble(v, 0, 2, v, 0, 2, 1.0);
- }
- }
- for(std::set<MVertex *>::iterator itv = fillNodes.begin(); itv !=fillNodes.end() ; ++itv){
- MVertex *v = *itv;
- myAssembler.assemble(v, 0, 1, v, 0, 1, small);
- myAssembler.assemble(v, 0, 2, v, 0, 2, small);
- }
-
//mettre max NC contraintes par bord
int NB = ordered.size();
- int NC = std::min(200,NB);
+ int NC = std::min(50,NB);
int jump = (int) NB/NC;
int nbLoop = (int) NB/jump ;
//printf("nb bound nodes=%d jump =%d \n", NB, jump);
for (int i = 0; i< nbLoop; i++){
MVertex *v1 = ordered[i*jump];
+ myAssembler.assemble(v1, 0, 1, v1, 0, 1, 1.0);
+ myAssembler.assemble(v1, 0, 2, v1, 0, 2, 1.0);
for (int j = 0; j< nbLoop; j++){
MVertex *v2 = ordered[j*jump];
- myAssembler.assemble(v1, 0, 1, v2, 0, 1, -1./NB);
- myAssembler.assemble(v1, 0, 2, v2, 0, 2, -1./NB);
+ myAssembler.assemble(v1, 0, 1, v2, 0, 1, -1./NC);
+ myAssembler.assemble(v1, 0, 2, v2, 0, 2, -1./NC);
}
}
@@ -1533,8 +1529,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
}
bool GFaceCompound::parametrize_conformal() const
{
- Msg::Info("Conformal parametrization of type: finite element");
-
+
dofManager<double> myAssembler(_lsys);
std::vector<MVertex*> ordered;
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 2d190ef295..26f77b6f58 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -76,7 +76,7 @@ class GFaceCompound : public GFace {
mutable std::set<MVertex*> fillNodes;
void buildOct() const ;
void buildAllNodes() const;
- void parametrize(iterationStep, typeOfMapping, double alpha=0.) const;
+ void parametrize(iterationStep, typeOfMapping) const;
bool parametrize_conformal() const;
bool parametrize_conformal_spectral() const;
bool parametrize_conformal_nonLinear() const;
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index fdc739c452..ec8584a60f 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -2601,7 +2601,6 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
Msg::Warning("VTK reader only accepts float or double datasets");
return 0;
}
- printf("POINTS=%d \n", numVertices);
Msg::Info("Reading %d points", numVertices);
std::vector<MVertex*> vertices(numVertices);
for(int i = 0 ; i < numVertices; i++){
@@ -2621,7 +2620,6 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
else{
if(fscanf(fp, "%lf %lf %lf", &xyz[0], &xyz[1], &xyz[2]) != 3) return 0;
}
- printf("Point(%d)={%g, %g, %g}; \n", i+1, xyz[0], xyz[1], xyz[2]);
vertices[i] = new MVertex(xyz[0], xyz[1], xyz[2]);
}
@@ -2721,7 +2719,6 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
int iLine = 1;
int num = 0;
for (int k= 0; k < numElements; k++){
- printf("*** new line =%d \n", iLine);
physicals[1][iLine][1] = "centerline";
fgets(line, sizeof(line), fp);
v0=(int)strtol(line, &pEnd, 10);//ignore firt line
@@ -2730,7 +2727,6 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
while(1){
v1 = strtol(p, &pEnd, 10);
if (p == pEnd ) break;
- printf("Line()={%d,%d};\n", vertices[v0]->getNum(),vertices[v1]->getNum());
elements[1][iLine].push_back(new MLine(vertices[v0],vertices[v1]));
p=pEnd;
v0=v1;
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index 392b40e8c6..507dda91b2 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -185,7 +185,6 @@ PView *GMSH_DistancePlugin::execute(PView *v)
for(unsigned int i = 0; i < _entities.size()-1; i++)
numnodes += _entities[i]->mesh_vertices.size();
int totNodes=numnodes + _entities[_entities.size()-1]->mesh_vertices.size();
- printf("%d\n",totNodes);
int order=ge->getMeshElement(0)->getPolynomialOrder();
int totNumNodes = totNodes+ge->getNumMeshElements()*integrationPointTetra[order-1];
diff --git a/Solver/crossConfTerm.h b/Solver/crossConfTerm.h
index dfbae205a4..011e5306fd 100644
--- a/Solver/crossConfTerm.h
+++ b/Solver/crossConfTerm.h
@@ -19,11 +19,13 @@ class crossConfTerm : public femTerm<double> {
const simpleFunction<double> *_diffusivity;
const int _iFieldR;
int _iFieldC;
+ std::map<MVertex*, SPoint3> *_coordView;
public:
crossConfTerm(GModel *gm, int iFieldR, int iFieldC,
- simpleFunction<double> *diffusivity, dofManager<double> *dofView=NULL)
+ simpleFunction<double> *diffusivity,
+ std::map<MVertex*, SPoint3> *coord=NULL)
: femTerm<double>(gm), _diffusivity(diffusivity), _iFieldR(iFieldR),
- _iFieldC(iFieldC) {}
+ _iFieldC(iFieldC), _coordView(coord) {}
virtual int sizeOfR(SElement *se) const
{
@@ -89,8 +91,30 @@ class crossConfTerm : public femTerm<double> {
}
void elementVector(SElement *se, fullVector<double> &m) const
{
- //adding here rhs
- printf("implment rhs cross term here\n");
+
+ MElement *e = se->getMeshElement();
+ int nbNodes = e->getNumVertices();
+
+ fullMatrix<double> *mat;
+ mat = new fullMatrix<double>(nbNodes,nbNodes);
+ elementMatrix(se, *mat);
+
+ fullVector<double> val(nbNodes);
+ val.scale(0.);
+ for (int i = 0; i < nbNodes; i++){
+ std::map<MVertex*, SPoint3>::iterator it = _coordView->find(e->getVertex(i));
+ SPoint3 UV = it->second;
+ if (_iFieldC == 1) val(i) = UV.x();
+ else if (_iFieldC == 2) val(i) = UV.y();
+ }
+
+ m.scale(0.);
+ for (int i = 0; i < nbNodes; i++)
+ for (int j = 0; j < nbNodes; j++)
+ m(i) += -(*mat)(i,j)*val(j);
+
+
+
}
};
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index dac0349778..8c8f08564a 100644
--- a/Solver/eigenSolver.cpp
+++ b/Solver/eigenSolver.cpp
@@ -75,15 +75,15 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
// print info
const EPSType type;
_try(EPSGetType(eps, &type));
- Msg::Info("SLEPc solution method: %s", type);
+ Msg::Debug("SLEPc solution method: %s", type);
PetscInt nev;
_try(EPSGetDimensions(eps, &nev, PETSC_NULL, PETSC_NULL));
- Msg::Info("SLEPc number of requested eigenvalues: %d", nev);
+ Msg::Debug("SLEPc number of requested eigenvalues: %d", nev);
PetscReal tol;
PetscInt maxit;
_try(EPSGetTolerances(eps, &tol, &maxit));
- Msg::Info("SLEPc stopping condition: tol=%g, maxit=%d", tol, maxit);
+ Msg::Debug("SLEPc stopping condition: tol=%g, maxit=%d", tol, maxit);
// solve
Msg::Info("SLEPc solving...");
@@ -97,7 +97,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
_try(EPSGetConvergedReason(eps, &reason));
if(reason == EPS_CONVERGED_TOL){
double t2 = Cpu();
- Msg::Info("SLEPc converged in %d iterations (%g s)", its, t2-t1);
+ Msg::Debug("SLEPc converged in %d iterations (%g s)", its, t2-t1);
}
else if(reason == EPS_DIVERGED_ITS)
Msg::Error("SLEPc diverged after %d iterations", its);
@@ -109,7 +109,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
// get number of converged approximate eigenpairs
PetscInt nconv;
_try(EPSGetConverged(eps, &nconv));
- Msg::Info("SLEPc number of converged eigenpairs: %d", nconv);
+ Msg::Debug("SLEPc number of converged eigenpairs: %d", nconv);
// ignore additional eigenvalues if we get more than what we asked
if(nconv > nev) nconv = nev;
@@ -118,8 +118,8 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
Vec xr, xi;
_try(MatGetVecs(A, PETSC_NULL, &xr));
_try(MatGetVecs(A, PETSC_NULL, &xi));
- Msg::Info(" Re[EigenValue] Im[EigenValue]"
- " Relative error");
+ Msg::Debug(" Re[EigenValue] Im[EigenValue]"
+ " Relative error");
for (int i = 0; i < nconv; i++){
PetscScalar kr, ki;
_try(EPSGetEigenpair(eps, i, &kr, &ki, xr, xi));
@@ -132,8 +132,8 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
PetscReal re = kr;
PetscReal im = ki;
#endif
- Msg::Info("EIG %03d %s%.16e %s%.16e %3.6e",
- i, (re < 0) ? "" : " ", re, (im < 0) ? "" : " ", im, error);
+ Msg::Debug("EIG %03d %s%.16e %s%.16e %3.6e",
+ i, (re < 0) ? "" : " ", re, (im < 0) ? "" : " ", im, error);
// store eigenvalues and eigenvectors
_eigenValues.push_back(std::complex<double>(re, im));
@@ -158,11 +158,11 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
_try(SlepcFinalize());
if(reason == EPS_CONVERGED_TOL){
- Msg::Info("SLEPc done");
+ Msg::Debug("SLEPc done");
return true;
}
else{
- Msg::Info("SLEPc failed");
+ Msg::Warning("SLEPc failed");
return false;
}
diff --git a/Solver/laplaceTerm.h b/Solver/laplaceTerm.h
index e7a1ceeeb1..7c8f29430b 100644
--- a/Solver/laplaceTerm.h
+++ b/Solver/laplaceTerm.h
@@ -10,13 +10,37 @@
// \nabla \cdot k \nabla U
class laplaceTerm : public helmholtzTerm<double> {
+ protected:
+ std::map<MVertex*, SPoint3> *_coordView;
+ const int _iField;
public:
- laplaceTerm(GModel *gm, int iField, simpleFunction<double> *k, dofManager<double> *dofView=NULL)
- : helmholtzTerm<double>(gm, iField, iField, k, 0) {}
+ laplaceTerm(GModel *gm, int iField, simpleFunction<double> *k,
+ std::map<MVertex*, SPoint3> *coord=NULL)
+ : helmholtzTerm<double>(gm, iField, iField, k, 0), _iField(iField), _coordView(coord) {}
void elementVector(SElement *se, fullVector<double> &m) const
{
- //adding here rhs
- printf("implment rhs laplace term here\n");
+
+ MElement *e = se->getMeshElement();
+ int nbNodes = e->getNumVertices();
+
+ fullMatrix<double> *mat;
+ mat = new fullMatrix<double>(nbNodes,nbNodes);
+ elementMatrix(se, *mat);
+
+ fullVector<double> val(nbNodes);
+ val.scale(0.);
+ for (int i = 0; i < nbNodes; i++){
+ std::map<MVertex*, SPoint3>::iterator it = _coordView->find(e->getVertex(i));
+ SPoint3 UV = it->second;
+ if (_iField == 1) val(i) = UV.x();
+ else if (_iField == 2) val(i) = UV.y();
+ }
+
+ m.scale(0.);
+ for (int i = 0; i < nbNodes; i++)
+ for (int j = 0; j < nbNodes; j++)
+ m(i) += -(*mat)(i,j)*val(j);
+
}
};
--
GitLab