From e3eea8908b19e5e0544a1b393c6ffc33c3c80d41 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Mon, 19 Oct 2009 15:35:25 +0000
Subject: [PATCH] Fixed bug conformal map
---
Geo/GFaceCompound.cpp | 123 ++++++++++++++++++++++++---------------
Geo/GFaceCompound.h | 1 +
Geo/GModel.cpp | 19 +++---
Geo/MElement.h | 2 +-
Geo/discreteEdge.cpp | 1 +
Mesh/meshGFace.cpp | 77 +++++++++++++++---------
Mesh/meshGFace.h | 2 +-
Mesh/meshPartition.cpp | 7 +--
Solver/crossConfTerm.h | 8 ++-
benchmarks/2d/square.geo | 2 +-
10 files changed, 152 insertions(+), 90 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 7f75d2627f..f67dc6119e 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -29,6 +29,7 @@
#include "Context.h"
#include "discreteFace.h"
+
static void fixEdgeToValue(GEdge *ed, double value, dofManager<double, double> &myAssembler)
{
myAssembler.fixVertex(ed->getBeginVertex()->mesh_vertices[0], 0, 1, value);
@@ -193,6 +194,55 @@ bool GFaceCompound::trivial() const
return false;
}
+void GFaceCompound::partitionFaceCM()
+{
+
+ double CMu = 0.0;
+ double sumArea = 0.0;
+
+ std::list<GFace*>::const_iterator it = _compound.begin();
+ for( ; it != _compound.end() ; ++it){
+ for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ MTriangle *t = (*it)->triangles[i];
+ std::map<MVertex*,SPoint3>::const_iterator it0 = coordinates.find(t->getVertex(0));
+ std::map<MVertex*,SPoint3>::const_iterator it1 = coordinates.find(t->getVertex(1));
+ std::map<MVertex*,SPoint3>::const_iterator it2 = coordinates.find(t->getVertex(2));
+ double q0[3] = {it0->second.x(), it0->second.y(), 0.0};
+ double q1[3] = {it1->second.x(), it1->second.y(), 0.0};
+ double q2[3] = {it2->second.x(), it2->second.y(), 0.0};
+ double area = 1/fabs(triangle_area(q0, q1, q2));
+ double cg_u = (q0[0]+q1[0]+q2[0])/3.;
+ CMu += cg_u*area;
+ sumArea += area;
+ }
+ CMu /= sumArea;
+ }
+
+ //printf("min size partition =%d \n", (int)allNodes.size()/2);
+ model()->setMinPartitionSize((int)allNodes.size()/2);
+ model()->setMaxPartitionSize((int)allNodes.size()/2+1);
+
+ it = _compound.begin();
+ for( ; it != _compound.end() ; ++it){
+ for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ MTriangle *t = (*it)->triangles[i];
+ std::map<MVertex*,SPoint3>::const_iterator it0 = coordinates.find(t->getVertex(0));
+ std::map<MVertex*,SPoint3>::const_iterator it1 = coordinates.find(t->getVertex(1));
+ std::map<MVertex*,SPoint3>::const_iterator it2 = coordinates.find(t->getVertex(2));
+ double cg_u = (it0->second.x()+it1->second.x()+it2->second.x())/3.;
+ if (cg_u <= CMu)t->setPartition(1);
+ else t->setPartition(2);
+ }
+ }
+
+ model()->recomputeMeshPartitions();
+
+ CreateOutputFile("toto.msh", CTX::instance()->mesh.format);
+ Msg::Exit(1);
+
+ return;
+}
+
//check if the discrete harmonic map is correct
//by checking that all the mapped triangles have
//the same normal orientation
@@ -339,25 +389,29 @@ void GFaceCompound::parametrize() const
if(!oct){
coordinates.clear();
+
+ if(allNodes.empty()) buildAllNodes();
//Laplace parametrization
//-----------------
if (_mapping == HARMONIC){
+ Msg::Info("Parametrizing surface %d with 'harmonic map'", tag());
parametrize(ITERU);
parametrize(ITERV);
}
//Conformal map parametrization
//-----------------
else if (_mapping == CONFORMAL){
+ Msg::Info("Parametrizing surface %d with 'conformal map'", tag());
parametrize_conformal();
}
//Distance function
//-----------------
//compute_distance();
-
+
checkOrientation();
-
buildOct();
+
computeNormals();
}
}
@@ -722,8 +776,7 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
void GFaceCompound::parametrize(iterationStep step) const
{
- Msg::Info("Parametrizing Surface %d coordinate (ITER %d)", tag(), step);
-
+
dofManager<double, double> myAssembler(_lsys);
simpleFunction<double> ONE(1.0);
@@ -796,13 +849,10 @@ void GFaceCompound::parametrize(iterationStep step) const
_lsys->systemSolve();
Msg::Debug("System solved");
-
- if(allNodes.empty()) buildAllNodes();
-
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value = myAssembler.getDofValue(v, 0, 1);
- std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
+ std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
if(itf == coordinates.end()){
SPoint3 p(0, 0, 0);
p[step] = value;
@@ -812,15 +862,14 @@ void GFaceCompound::parametrize(iterationStep step) const
itf->second[step]= value;
}
}
+
_lsys->clear();
+
}
void GFaceCompound::parametrize_conformal() const
{
- if(oct) return;
- Msg::Info("Parametrizing Surface %d", tag());
-
dofManager<double, double> myAssembler(_lsys);
std::vector<MVertex*> ordered;
@@ -885,28 +934,17 @@ void GFaceCompound::parametrize_conformal() const
_lsys->systemSolve();
Msg::Debug("System solved");
- it = _compound.begin();
- std::set<MVertex *>allNodes;
- for( ; it != _compound.end() ; ++it){
- for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- MTriangle *t = (*it)->triangles[i];
- for(int j = 0; j < 3; j++){
- allNodes.insert(t->getVertex(j));
- }
- }
- }
-
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value1 = myAssembler.getDofValue(v,0,1);
double value2 = myAssembler.getDofValue(v,0,2);
coordinates[v] = SPoint3(value1,value2,0.0);
}
+ // printStuff();
+ //exit(1);
- checkOrientation();
- computeNormals();
- buildOct();
_lsys->clear();
+
}
void GFaceCompound::compute_distance() const
@@ -946,18 +984,7 @@ void GFaceCompound::compute_distance() const
Msg::Info("Distance Computation: Assembly done");
_lsys->systemSolve();
Msg::Info("Distance Computation: System solved");
-
- it = _compound.begin();
- std::set<MVertex *>allNodes;
- for( ; it != _compound.end() ; ++it){
- for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- MTriangle *t = (*it)->triangles[i];
- for(int j = 0; j < 3; j++){
- allNodes.insert(t->getVertex(j));
- }
- }
- }
-
+
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value = std::min(0.9999, myAssembler.getDofValue(v, 0, 1));
@@ -1314,7 +1341,7 @@ bool GFaceCompound::checkAspectRatio() const
{
parametrize();
-
+
if (!_checkedAR){
bool paramOK = true;
@@ -1325,12 +1352,16 @@ bool GFaceCompound::checkAspectRatio() const
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
- std::map<MVertex*,SPoint3>::const_iterator it0 = coordinates.find(t->getVertex(0));
- std::map<MVertex*,SPoint3>::const_iterator it1 = coordinates.find(t->getVertex(1));
- std::map<MVertex*,SPoint3>::const_iterator it2 = coordinates.find(t->getVertex(2));
- double p0[3] = {t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()};
- double p1[3] = {t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()};
- double p2[3] = {t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z()};
+ std::vector<MVertex *> v(3);
+ for(int k = 0; k < 3; k++){
+ v[k] = t->getVertex(k);
+ }
+ std::map<MVertex*,SPoint3>::const_iterator it0 = coordinates.find(v[0]);
+ std::map<MVertex*,SPoint3>::const_iterator it1 = coordinates.find(v[1]);
+ std::map<MVertex*,SPoint3>::const_iterator it2 = coordinates.find(v[2]);
+ double p0[3] = {v[0]->x(), v[0]->y(), v[0]->z()};
+ double p1[3] = {v[1]->x(), v[1]->y(), v[1]->z()};
+ double p2[3] = {v[2]->x(), v[2]->y(), v[2]->z()};
double a_3D = fabs(triangle_area(p0, p1, p2));
double q0[3] = {it0->second.x(), it0->second.y(), 0.0};
double q1[3] = {it1->second.x(), it1->second.y(), 0.0};
@@ -1340,7 +1371,7 @@ bool GFaceCompound::checkAspectRatio() const
}
}
- if (areaMax > 1.e12) {
+ if (areaMax > 1.e15) {
nbSplit = 2;
printf("--> Geometrical aspect ratio too high (1/area_2D=%g)\n", areaMax);
printf("--> Partition geometry in N=%d parts\n", nbSplit);
@@ -1350,7 +1381,7 @@ bool GFaceCompound::checkAspectRatio() const
Msg::Info("Geometrical aspect ratio (1/area_2D=%g)", areaMax);
paramOK = true;
}
-
+
_checkedAR = true;
_paramOK = paramOK;
}
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index d675d23ee1..d7fb7e3e94 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -100,6 +100,7 @@ class GFaceCompound : public GFace {
virtual double curvatureMax(const SPoint2 ¶m) const;
virtual int genusGeom () const;
virtual bool checkTopology() const;
+ void partitionFaceCM();
virtual std::list<GFace*> getCompounds() const {return _compound;};
mutable int nbSplit;
mutable bool _checkedAR;
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 35db317475..54ca87f0e6 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -723,7 +723,7 @@ void GModel::recomputeMeshPartitions()
for(unsigned int i = 0; i < entities.size(); i++){
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
int part = entities[i]->getMeshElement(j)->getPartition();
- if(part) meshPartitions.insert(part);
+ if(part) meshPartitions.insert(part);
}
}
}
@@ -1056,9 +1056,10 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discreteFaces)
{
std::vector<discreteEdge*> discreteEdges;
- for(eiter it = firstEdge(); it != lastEdge(); it++)
+ for(eiter it = firstEdge(); it != lastEdge(); it++){
if((*it)->geomType() == GEntity::DiscreteCurve)
discreteEdges.push_back((discreteEdge*) *it);
+ }
int initSizeEdges = discreteEdges.size();
// find boundary edges of each face and put them in a map_edges that
@@ -1075,7 +1076,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discreteFaces)
std::map<int, std::vector<int> > face2Edges;
while (!map_edges.empty()){
-
+
std::vector<MEdge> myEdges;
std::vector<int> tagFaces = map_edges.begin()->second;
myEdges.push_back(map_edges.begin()->first);
@@ -1111,18 +1112,20 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discreteFaces)
}
}
else {
+ //TODO EMI: check if edge correct by comparing two sets of MEdges
for(unsigned int i = 0; i < myEdges.size(); i++){
if (myEdges[i].getVertex(0)->onWhat()->dim() == 1) {
- int tagEdge = myEdges[i].getVertex(0)->onWhat()->tag();
+ GEntity *ge = myEdges[i].getVertex(0)->onWhat();
+ int tagEdge = ge->tag();
std::vector<int>::iterator itv = std::find(tagEdges.begin(),
tagEdges.end(), tagEdge);
- if (itv == tagEdges.end())
+ int nbElems = ge->getNumMeshElements();
+ if (itv == tagEdges.end() )// && myEdges.size() == nbElems)
tagEdges.push_back(tagEdge);
}
}
}
- for (std::vector<int>::iterator itFace = tagFaces.begin();
- itFace != tagFaces.end(); itFace++) {
+ for (std::vector<int>::iterator itFace = tagFaces.begin(); itFace != tagFaces.end(); itFace++) {
std::map<int, std::vector<int> >::iterator it = face2Edges.find(*itFace);
if (it == face2Edges.end()){
std::vector<int> allEdges;
@@ -1180,6 +1183,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discreteFaces)
}
discreteEdge *e = new discreteEdge(this, num, 0, 0);
+ printf("create new discrete edge =%d \n", num);
add(e);
discreteEdges.push_back(e);
@@ -1221,6 +1225,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discreteFaces)
// set boundary edges for each face
for (std::vector<discreteFace*>::iterator it = discreteFaces.begin(); it != discreteFaces.end(); it++){
+ //printf("set boundary edge for face =%d \n", (*it)->tag());
std::map<int, std::vector<int> >::iterator ite = face2Edges.find((*it)->tag());
if (ite != face2Edges.end()){
std::vector<int> myEdges = ite->second;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 46107c3bb2..a35759b733 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -64,7 +64,7 @@ class MElement
// get/set the partition to which the element belongs
virtual int getPartition(){ return _partition; }
- virtual void setPartition(int num){ _partition = (short)num; }
+ virtual void setPartition(int num){_partition = (short)num; }
// get/set the visibility flag
virtual char getVisibility();
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 207b05f30d..9893626554 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -163,6 +163,7 @@ void discreteEdge::setBoundVertices()
}
if(!existVertex){
GVertex *gvB = new discreteVertex(model(), model()->maxVertexNum()+1);
+ printf("create new discreteVertex =%d for edge =%d\n", gvB->tag(), this->tag());
bound_vertices.push_back(gvB);
vE->setEntity(gvB);
gvB->mesh_vertices.push_back(vE);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index fe926b5df1..bbb520d44f 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -33,6 +33,8 @@
#include "meshGEdge.h"
#include "meshPartitionOptions.h"
#include "meshPartition.h"
+#include "CreateFile.h"
+#include "Context.h"
void fourthPoint(double *p1, double *p2, double *p3, double *p4)
{
@@ -243,7 +245,6 @@ static bool recoverEdge(BDS_Mesh *m, GEdge *ge,
// Builds An initial triangular mesh that respects the boundaries of
// the domain, including embedded points and surfaces
-
static bool meshGenerator(GFace *gf, int RECUR_ITER,
bool repairSelfIntersecting1dMesh,
bool debug = true)
@@ -257,7 +258,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
if(gf->geomType() == GEntity::CompoundSurface){
bool correct = gf->checkTopology();
if (!correct){
- partitionAndRemesh((GFaceCompound*) gf);
+ partitionAndRemesh((GFaceCompound*) gf, 1);
return true;
}
std::set<GEdge*> mySet;
@@ -336,7 +337,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
bool paramOK = true;
paramOK = reparamMeshVertexOnFace(here, gf, param);
if(gf->geomType() == GEntity::CompoundSurface && !paramOK ){
- partitionAndRemesh((GFaceCompound*)gf);
+ partitionAndRemesh((GFaceCompound*)gf, 2);
return true;
}
BDS_Point *pp = m->add_point(count, param[0], param[1], gf);
@@ -1307,9 +1308,18 @@ void meshGFace::operator() (GFace *gf)
gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size());
}
-void partitionAndRemesh(GFaceCompound *gf)
+void partitionAndRemesh(GFaceCompound *gf, int algo)
{
+
#if defined(HAVE_CHACO) || defined(HAVE_METIS)
+
+ //Create list of faces to partition
+ //----------------------------------
+// std::list<GFace*> cFaces = gf->getCompounds();
+// GModel *tmp_model = new GModel();
+// for(std::list<GFace*>::iterator it = cFaces.begin(); it != cFaces.end(); it++)
+// tmp_model->add(*it);
+
//Partition the mesh and createTopology for new faces
//-----------------------------------------------------
int N = gf->nbSplit;
@@ -1317,38 +1327,40 @@ void partitionAndRemesh(GFaceCompound *gf)
options = CTX::instance()->partitionOptions;
options.num_partitions = N;
options.partitioner = 2; //METIS
- options.algorithm = 1;
- int ier = PartitionMesh(gf->model(), options);
+ options.algorithm = 1 ;
+
+ //if (algo == 1)
+ PartitionMesh(gf->model(), options);
+// else if (algo ==2)
+// gf->partitionFaceCM();
+
int numv = gf->model()->maxVertexNum() + 1;
int nume = gf->model()->maxEdgeNum() + 1;
int numf = gf->model()->maxFaceNum() + 1;
std::vector<discreteFace*> pFaces;
createPartitionFaces(gf->model(), gf, N, pFaces);
- gf->model()->createTopologyFromMesh(); //Faces(pFaces);
- //gf->model()->createTopologyFromFaces(pFaces);
- //CreateOutputFile(CTX::instance()->outputFileName, CTX::instance()->mesh.format);
- //Msg::Exit(1);
-
+ gf->model()->createTopologyFromFaces(pFaces);
+
//Remesh new faces (Compound Lines and Compound Surfaces)
//-----------------------------------------------------
+ //Remeshing discrete Edges
int numEdges = gf->model()->maxEdgeNum() - nume + 1;
- printf("*** Created %d discretEdges \n", numEdges);
+ printf("*** Created %d discreteEdges \n", numEdges);
for (int i=0; i < numEdges; i++){
std::vector<GEdge*>e_compound;
GEdge *pe = gf->model()->getEdgeByTag(nume+i);//partition edge
int num_gec = gf->model()->maxEdgeNum() + 1;
e_compound.push_back(pe);
printf("*** Remeshing discreteEdge %d with CompoundLine %d\n", pe->tag(), num_gec);
- GEdge *gec = new GEdgeCompound(gf->model(), num_gec, e_compound);
-
+ GEdge *gec = new GEdgeCompound(gf->model(), num_gec, e_compound);
meshGEdge mge;
mge(gec);//meshing 1D
-
gf->model()->remove(pe);
}
+ //Removing discrete Vertex
int numVert = gf->model()->maxVertexNum() - numv + 1;
printf("*** Created %d discreteVert \n", numVert);
for (int i=0; i < numEdges; i++){
@@ -1356,13 +1368,13 @@ void partitionAndRemesh(GFaceCompound *gf)
gf->model()->remove(pv);
}
+ //Remeshing discrete Face
std::list<GEdge*> b[4];
- std::set<MVertex*> allNod;
-
+ std::set<MVertex*> allNod;
for (int i=0; i < N; i++){
std::list<GFace*> f_compound;
GFace *pf = gf->model()->getFaceByTag(numf+i);//partition face
- //discreteFace *pf = pFaces[i];//partition face
+ //GFace *pf = pFaces[i];//partition face
int num_gfc = numf + N + i ;
f_compound.push_back(pf);
@@ -1371,17 +1383,18 @@ void partitionAndRemesh(GFaceCompound *gf)
meshGFace mgf;
mgf(gfc);//meshing 2D
- for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
- MTriangle *t = gfc->triangles[j];
- std::vector<MVertex *> v(3);
- for(int k = 0; k < 3; k++){
- v[k] = t->getVertex(k);
- allNod.insert(v[k]);
- }
- gf->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
+ for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
+ MTriangle *t = gfc->triangles[j];
+ std::vector<MVertex *> v(3);
+ for(int k = 0; k < 3; k++){
+ v[k] = t->getVertex(k);
+ allNod.insert(v[k]);
}
+ gf->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
+ }
- gf->model()->remove(pf);
+ gf->model()->remove(pf);
+
}
//Put new mesh in a new discreteFace
@@ -1411,6 +1424,16 @@ void partitionAndRemesh(GFaceCompound *gf)
printf("*** Mesh of surface %d done by assembly remeshed faces\n", gf->tag());
gf->meshStatistics.status = GFace::DONE;
+
+// for(std::list<GFace*>::iterator it = cFaces.begin(); it != cFaces.end(); it++)
+// tmp_model->remove(*it);
+// delete tmp_model;
+
+ //CreateOutputFile("toto.msh", CTX::instance()->mesh.format);
+ //Msg::Exit(1);
+
+#else
+ return;
#endif
}
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index 1d7a09d36e..9cd5a2bdf6 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -45,6 +45,6 @@ void findTransfiniteCorners(GFace *gf, std::vector<MVertex*> &corners);
int MeshTransfiniteSurface(GFace *gf);
int MeshExtrudedSurface(GFace *gf, std::set<std::pair<MVertex*, MVertex*> >
*constrainedEdges=0);
-void partitionAndRemesh(GFaceCompound *gf);
+void partitionAndRemesh(GFaceCompound *gf, int algo);
#endif
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index f5d599c546..7b4266876c 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -1009,18 +1009,17 @@ int CreatePartitionBoundaries(GModel *model)
}
void createPartitionFaces(GModel *model, GFaceCompound *gf, int N,
- std::vector<discreteFace*> &discrFaces)
+ std::vector<discreteFace*> &discreteFaces)
{
printf("---> CreateTopologyFromPartition for Compound Face %d \n", gf->tag());
- std::vector<discreteFace*> discreteFaces;//delete this
-
+
// Compound is partitioned in N discrete faces
//--------------------------------------------
std::vector<std::set<MVertex*> > allNodes;
int numMax = model->maxFaceNum() + 1;
for( int i =0; i < N; i++){
- //printf("*** Created discreteFace %d \n", numMax+i);
+ printf("*** Created discreteFace %d \n", numMax+i);
discreteFace *face = new discreteFace(model, numMax+i);
discreteFaces.push_back(face);
model->add(face);//delete this
diff --git a/Solver/crossConfTerm.h b/Solver/crossConfTerm.h
index c8d313d7cf..2b0488321c 100644
--- a/Solver/crossConfTerm.h
+++ b/Solver/crossConfTerm.h
@@ -34,11 +34,13 @@ class crossConfTerm : public femTerm<double, double> {
}
Dof getLocalDofR(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldR);
+ return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ Dof::createTypeWithTwoInts(0, _iFieldR));
}
Dof getLocalDofC(SElement *se, int iRow) const
{
- return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldC);
+ return Dof(se->getMeshElement()->getVertex(iRow)->getNum(),
+ Dof::createTypeWithTwoInts(0, _iFieldC));
}
virtual void elementMatrix(SElement *se, fullMatrix<double> &m) const
{
@@ -82,7 +84,7 @@ class crossConfTerm : public femTerm<double, double> {
}
for (int j = 0; j < nbNodes; j++)
for (int k = 0; k < j; k++)
- m(k, j) = -1. * m(j, k);
+ m(k, j) = -1.* m(j, k);
}
};
diff --git a/benchmarks/2d/square.geo b/benchmarks/2d/square.geo
index 6475d2fad2..88dbb82010 100644
--- a/benchmarks/2d/square.geo
+++ b/benchmarks/2d/square.geo
@@ -11,5 +11,5 @@ Line(4) = {1, 2};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(10) = {5};
-Compound Surface(11)={10};
+Compound Surface(-11)={10};
--
GitLab