diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 84ec84deffde046c6065dc3306d0a743e8506645..af37c8d204fe0c69673c57f67ac259ddf45febff 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -188,10 +188,50 @@ static void myPolygon(std::vector<MElement*> &vTri, std::vector<MVertex*> &vPoly
// }
}
+bool checkCavity(std::vector<MElement*> &vTri, std::map<MVertex*, SPoint2> &vCoord) {
+
+ bool badCavity = false;
+
+ unsigned int nbV = vTri.size();
+ double a_old = 0, a_new;
+ for(unsigned int i = 0; i < nbV; ++i){
+ MTriangle *t = (MTriangle*) vTri[i];
+ SPoint2 v1 = vCoord[t->getVertex(0)];
+ SPoint2 v2 = vCoord[t->getVertex(1)];
+ SPoint2 v3 = vCoord[t->getVertex(2)];
+ double p1[2] = {v1[0],v1[1]};
+ double p2[2] = {v2[0],v2[1]};
+ double p3[2] = {v3[0],v3[1]};
+ a_new = robustPredicates::orient2d(p1, p2, p3);
+ if(i == 0) a_old=a_new;
+ if(a_new*a_old < 0.) badCavity = true;
+ a_old = a_new;
+ }
+
+ return badCavity;
+}
+
+static bool closedCavity(MVertex *v, std::vector<MElement*> &vTri){
+ std::set<MVertex *> vs;
+ for (unsigned int i = 0; i < vTri.size(); i++){
+ MElement *t = vTri[i];
+ for (int j = 0; j < t->getNumVertices(); j++){
+ MVertex *vv = t->getVertex(j);
+ if (vv != v){
+ if (vs.find(vv) == vs.end())vs.insert(vv);
+ else vs.erase(vv);
+ }
+ }
+ }
+ return vs.empty();
+}
+
void GFaceCompound::fillNeumannBCS() const
{
+ fillTris.clear();
+
//close neuman bcs
for(std::list<std::list<GEdge*> >::const_iterator iloop = _interior_loops.begin();
iloop != _interior_loops.end(); iloop++){
@@ -220,25 +260,47 @@ void GFaceCompound::fillNeumannBCS() const
for (unsigned int i= 0; i< (*ite)->getNumMeshElements(); i++){
MVertex *v0 = (*ite)->getMeshElement(i)->getVertex(0);
MVertex *v1 = (*ite)->getMeshElement(i)->getVertex(1);
- MTriangle *myTri = new MTriangle(v0,v1, c);
- fillTris.push_back(myTri);
+
+ //fillTris.push_back(new MTriangle(v0,v1, c));
+
+// MVertex *v2 = new MVertex(.5*(v0->x()+c->x()), .5*(v0->y()+c->y()), .5*(v0->z()+c->z()));
+// MVertex *v3 = new MVertex(.5*(v1->x()+c->x()), .5*(v1->y()+c->y()), .5*(v1->z()+c->z()));
+// fillTris.push_back(new MTriangle(v0,v2,v3));
+// fillTris.push_back(new MTriangle(v2,c, v3));
+// fillTris.push_back(new MTriangle(v0,v3, v1);
+
+ MVertex *v2 = new MVertex(.66*v0->x()+.33*c->x(), .66*v0->y()+.33*c->y(), .66*v0->z()+.33*c->z());
+ MVertex *v3 = new MVertex(.66*v1->x()+.33*c->x(), .66*v1->y()+.33*c->y(), .66*v1->z()+.33*c->z());
+ MVertex *v4 = new MVertex(.33*v0->x()+.66*c->x(), .33*v0->y()+.66*c->y(), .33*v0->z()+.66*c->z());
+ MVertex *v5 = new MVertex(.33*v1->x()+.66*c->x(), .33*v1->y()+.66*c->y(), .33*v1->z()+.66*c->z());
+ fillTris.push_back(new MTriangle(v0,v2,v3));
+ fillTris.push_back(new MTriangle(v2,v5,v3));
+ fillTris.push_back(new MTriangle(v2,v4,v5));
+ fillTris.push_back(new MTriangle(v4,c,v5));
+ fillTris.push_back(new MTriangle(v0,v3,v1));
+
}
}
}
}
-// FILE * ftri = fopen("fillTris.pos","w");
-// fprintf(ftri,"View \"\"{\n");
-// for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
-// MTriangle *t = (*it2);
-// fprintf(ftri,"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(),
-// 1., 1., 1.);
-// }
-// fprintf(ftri,"};\n");
-// fclose(ftri);
+ //if (fillTris.size() > 0)
+ //printf("**** Filling Neuman BCs with %d triangles \n", fillTris.size());
+
+// FILE * ftri = fopen("fillTris.pos","w");
+// fprintf(ftri,"View \"\"{\n");
+// for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+// MTriangle *t = (*it2);
+// fprintf(ftri,"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(),
+// 1., 1., 1.);
+// }
+// fprintf(ftri,"};\n");
+// fclose(ftri);
+
+// exit(1);
}
@@ -252,53 +314,6 @@ 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
@@ -309,13 +324,11 @@ bool GFaceCompound::checkOrientation(int iter) const
//Only check orientation for stl files (1 patch)
// if(_compound.size() > 1.0) return true;
- //return true;
-
std::list<GFace*>::const_iterator it = _compound.begin();
- double a_old = 0, a_new;
+ double a_old = 0.0, a_new=0.0;
bool oriented = true;
+ int count = 0;
for( ; it != _compound.end(); ++it){
- int count = 0;
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
SPoint3 v1 = coordinates[t->getVertex(0)];
@@ -326,92 +339,57 @@ bool GFaceCompound::checkOrientation(int iter) const
double p2[2] = {v3[0],v3[1]};
a_new = robustPredicates::orient2d(p0, p1, p2);
if(count == 0) a_old=a_new;
- if(a_new*a_old < 0.){// && fabs(a_new) > 1.e10){
+ if(a_new*a_old < 0.){
oriented = false;
break;
}
else{
a_old = a_new;
}
- count ++;
+ count++;
}
}
- if(!oriented && iter < 10){
- if (iter == 0)Msg::Warning("*** Parametrization is NOT 1 to 1 : applying cavity checks.");
+ int iterMax = 5;
+ if(!oriented && iter < iterMax){
+ if (iter == 0) Msg::Warning("*** Parametrization is NOT 1 to 1 : applying cavity checks.");
Msg::Warning("*** Cavity Check - iter %d -",iter);
one2OneMap();
return checkOrientation(iter+1);
}
- else if (iter < 10){
- Msg::Info("Parametrization is 1 to 1 :-)");
+ else if (iter < iterMax){
+ Msg::Debug("Parametrization is 1 to 1 :-)");
}
return oriented;
}
-bool GFaceCompound::checkCavity(std::vector<MElement*> &vTri) const
-{
-
- bool badCavity = false;
-
- unsigned int nbV = vTri.size();
- double a_old = 0, a_new;
- for(unsigned int i = 0; i < nbV; ++i){
- MTriangle *t = (MTriangle*) vTri[i];
- SPoint3 v1 = coordinates[t->getVertex(0)];
- SPoint3 v2 = coordinates[t->getVertex(1)];
- SPoint3 v3 = coordinates[t->getVertex(2)];
- double p1[2] = {v1[0],v1[1]};
- double p2[2] = {v2[0],v2[1]};
- double p3[2] = {v3[0],v3[1]};
- //printf("p1=(%g %g) p2=(%g %g) p3=(%g %g)\n",v1[0],v1[1], v2[0],v2[1], v3[0],v3[1] );
- a_new = robustPredicates::orient2d(p1, p2, p3);
- if(i == 0) a_old=a_new;
- if(a_new*a_old < 0.) badCavity = true;
- a_old = a_new;
- }
-
- return badCavity;
-}
-
-static bool closedCavity(MVertex *v, std::vector<MElement*> &vTri){
- std::set<MVertex *> vs;
- for (unsigned int i = 0; i < vTri.size(); i++){
- MElement *t = vTri[i];
- for (int j = 0; j < t->getNumVertices(); j++){
- MVertex *vv = t->getVertex(j);
- if (vv != v){
- if (vs.find(vv) == vs.end())vs.insert(vv);
- else vs.erase(vv);
- }
- }
- }
- return vs.empty();
-}
-
void GFaceCompound::one2OneMap() const
{
#if defined(HAVE_MESH)
- if(!mapv2Tri){
+ if(adjv.size() == 0){
std::vector<MTriangle*> allTri;
std::list<GFace*>::const_iterator it = _compound.begin();
for( ; it != _compound.end(); ++it){
allTri.insert(allTri.end(), (*it)->triangles.begin(), (*it)->triangles.end() );
}
buildVertexToTriangle(allTri, adjv);
- mapv2Tri=true;
}
for(v2t_cont::iterator it = adjv.begin(); it!= adjv.end(); ++it){
MVertex *v = it->first;
- std::map<MVertex*,SPoint3>::iterator itV = coordinates.find(v);
-
- std::vector<MVertex*> cavV;
+ SPoint2 p=getCoordinates(v);
std::vector<MElement*> vTri = it->second;
- bool badCavity = closedCavity(v,vTri) ? checkCavity(vTri) : false;
-
+ std::map<MVertex*,SPoint2> vCoord;
+ for (int j=0; j < vTri.size(); j++){
+ for (int k= 0; k < vTri[j]->getNumVertices(); k++){
+ MVertex *vk = vTri[j]->getVertex(k);
+ vCoord[vk] = getCoordinates(vk);
+ }
+ }
+ bool badCavity = closedCavity(v,vTri) ? checkCavity(vTri, vCoord) : false;
+
if(badCavity){
Msg::Debug("Wrong cavity around vertex %d (onwhat=%d).",
v->getNum(), v->onWhat()->dim());
@@ -419,11 +397,12 @@ void GFaceCompound::one2OneMap() const
vTri.size());
double u_cg, v_cg;
+ std::vector<MVertex*> cavV;
myPolygon(vTri, cavV);
computeCGKernelPolygon(coordinates, cavV, u_cg, v_cg);
SPoint3 p_cg(u_cg,v_cg,0);
coordinates[v] = p_cg;
-
+
}
}
#endif
@@ -433,7 +412,6 @@ bool GFaceCompound::parametrize() const
{
bool paramOK = true;
-
if(oct) return paramOK;
if(trivial()) return paramOK;
@@ -441,17 +419,29 @@ bool GFaceCompound::parametrize() const
if(allNodes.empty()) buildAllNodes();
+// // TEST MULTISCALE LAPLACE ----------------------------------
+// std::vector<MElement*> _elements;
+// for( std::list<GFace*>::const_iterator itt = _compound.begin(); itt != _compound.end(); ++itt)
+// for(unsigned int i = 0; i < (*itt)->triangles.size(); ++i)
+// _elements.push_back((*itt)->triangles[i]);
+// multiscaleLaplace multiLaplace(_elements);
+// printf("ended multiscale \n");
+// exit(1);
+ // TEST MULTISCALE LAPLACE ----------------------------------
+
//Laplace parametrization
//-----------------
if (_mapping == HARMONIC){
- Msg::Info("Parametrizing surface %d with 'harmonic map'", tag());
+ Msg::Debug("Parametrizing surface %d with 'harmonic map'", tag());
+ fillNeumannBCS();
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
}
//Conformal map parametrization
//-----------------
else if (_mapping == CONFORMAL){
- Msg::Info("Parametrizing surface %d with 'conformal map'", tag());
+ Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
+ fillNeumannBCS();
parametrize_conformal();
}
//Distance function
@@ -459,11 +449,6 @@ bool GFaceCompound::parametrize() const
//compute_distance();
buildOct();
-
- if (!checkAspectRatio()){
- paramOK = false;
- return paramOK;
- }
if (!checkOrientation(0)){
Msg::Info("Parametrization failed using standard techniques : moving to convex combination");
@@ -475,8 +460,16 @@ bool GFaceCompound::parametrize() const
buildOct();
}
+ printStuff();
+
computeNormals();
+ if (!checkAspectRatio()){
+ printf("WARNING: geom aspect ratio too high \n");
+ exit(1);
+ paramOK = false;
+ }
+
return paramOK;
}
@@ -517,8 +510,8 @@ void GFaceCompound::getBoundingEdges()
(*itf)->addFace(this);
}
std::list<GEdge*> loop;
- computeALoop(_unique,loop); //_U0);
- while(!_unique.empty()) computeALoop(_unique, loop); //_U1);
+ computeALoop(_unique,loop);
+ while(!_unique.empty()) computeALoop(_unique, loop);
//assign Derichlet BC (_U0) to bound with largest size
double maxSize = 0.0;
@@ -537,6 +530,26 @@ void GFaceCompound::getBoundingEdges()
}
+double GFaceCompound::getSizeH() const
+{
+
+ SBoundingBox3d bb;
+ SOrientedBoundingBox obbox ;
+ std::vector<SPoint3> vertices;
+ for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ SPoint3 pt((*itv)->x(),(*itv)->y(), (*itv)->z());
+ vertices.push_back(pt);
+ bb +=pt;
+ }
+
+ //obbox = SOrientedBoundingBox::buildOBB(vertices);
+ //double H = obbox.getMaxSize();
+
+ double H = norm(SVector3(bb.max(), bb.min()));
+
+ return H;
+
+}
double GFaceCompound::getSizeBB(const std::list<GEdge* > &elist) const
{
@@ -609,6 +622,7 @@ tr.low());
void GFaceCompound::getUniqueEdges(std::set<GEdge*> &_unique)
{
_unique.clear();
+
// in case the bounding edges are not given Boundary { {} };
std::multiset<GEdge*> _touched;
std::list<GFace*>::iterator it = _compound.begin();
@@ -638,7 +652,6 @@ void GFaceCompound::computeALoop(std::set<GEdge*> &_unique, std::list<GEdge*> &l
while(!_unique.empty()) {
std::set<GEdge*>::iterator it = _unique.begin();
- //printf("for face %d bound edge =%d\n", tag(), (*it)->tag());
GVertex *vB = (*it)->getBeginVertex();
GVertex *vE = (*it)->getEndVertex();
_loop.push_back(*it);
@@ -716,8 +729,6 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
_lsys = new linearSystemFull<double>;
#endif
}
- nbSplit = 0;
- checked = false;
for(std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
if(!(*it)){
@@ -731,7 +742,9 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
else if(!_V1.size()) _type = UNITCIRCLE;
else _type = SQUARE;
- mapv2Tri = false;
+ nbSplit = 0;
+ fillTris.clear();
+
}
GFaceCompound::~GFaceCompound()
@@ -765,8 +778,6 @@ static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
l.clear();
coord.clear();
- //printf("%d lines\n",e.size());
-
std::list<GEdge*>::const_iterator it = e.begin();
std::list<MLine*> temp;
double tot_length = 0;
@@ -844,20 +855,17 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
//getParameter Point
v->getParameter(0,tGlob);
- //printf("*** global param for point (%g, %g, %g ) = %g\n", v->x(), v->y(), v->z(), tGlob);
-
+
//find compound Edge
GEdgeCompound *gec = dynamic_cast<GEdgeCompound*>(v->onWhat());
- //printf("tag compound=%d, beginvertex=%d, endvertex=%d \n", gec->tag(), gec->getBeginVertex()->tag(), gec->getEndVertex()->tag());
-
+
if(gec){
//compute local parameter on Edge
gec->getLocalParameter(tGlob,iEdge,tLoc);
std::vector<GEdge*> gev = gec->getEdgesOfCompound();
GEdge *ge = gev[iEdge];
- //printf("iEdge =%d, leftV =%d, rightV=%d, and local param=%g \n", ge->tag(), ge->getBeginVertex()->tag(), ge->getEndVertex()->tag(), tLoc);
-
+
//left and right vertex of the Edge
MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
MVertex *v1 = ge->getEndVertex()->mesh_vertices[0];
@@ -865,12 +873,10 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
std::map<MVertex*,SPoint3>::iterator itR = coordinates.find(v1);
//for the Edge, find the left and right vertices of the initial 1D mesh and interpolate to find (u,v)
- //printf("number of vertices =%d for Edge =%d \n", ge->mesh_vertices.size(), ge->tag());
MVertex *vL = v0;
MVertex *vR = v1;
double tB = ge->parBounds(0).low();
double tE = ge->parBounds(0).high();
- //printf("tB=%g uv (%g %g) tE=%g (%g %g), tLoc=%g\n", tB, itL->second.x(), itL->second.y(), tE, itR->second.x(), itR->second.y(), tLoc);
int j = 0;
tL=tB;
bool found = false;
@@ -881,7 +887,6 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
Msg::Error("vertex vr %p not medgevertex \n", vR);
Msg::Exit(1);
}
- //printf("***tLoc=%g tL=%g, tR=%g L=%p (%g,%g,%g) et R= %p (%g,%g,%g)\n", tLoc, tL, tR, vL, vL->x(),vL->y(),vL->z(),vR, vR->x(), vR->y(), vR->z());
if(tLoc >= tL && tLoc <= tR){
found = true;
itR = coordinates.find(vR);
@@ -897,23 +902,16 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
tL = tR;
}
j++;
- //printf("in while j =%d, vL (%g,%g,%g), -> uv= (%g,%g)\n",j, vL->x(), vL->y(), vL->z(), itL->second.x(), itL->second.y());
}
if(!found) {
vR = v1;
tR = tE;
}
- //printf("vL (%g,%g,%g), -> uv= (%g,%g)\n",vL->x(), vL->y(), vL->z(), itL->second.x(), itL->second.y());
- //printf("vR (%g,%g,%g), -> uv= (%g,%g)\n",vR->x(), vR->y(), vR->z(), itR->second.x(), itR->second.y());
- //printf("tL:%g, tR=%g, tLoc=%g \n", tL, tR, tLoc);
//Linear interpolation between tL and tR
double uloc, vloc;
uloc = itL->second.x() + (tLoc-tL)/(tR-tL) * (itR->second.x()-itL->second.x());
vloc = itL->second.y() + (tLoc-tL)/(tR-tL) * (itR->second.y()-itL->second.y());
- //printf("uloc=%g, vloc=%g \n", uloc,vloc);
- //exit(1);
-
return SPoint2(uloc,vloc);
}
}
@@ -922,8 +920,8 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
}
void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
-{
-
+{
+
dofManager<double> myAssembler(_lsys);
simpleFunction<double> ONE(1.0);
@@ -937,20 +935,6 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
return;
}
- // TEST --------------------------------------------------
- if(step == ITERV){
- std::vector<MElement*> _elements;
- 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];
- _elements.push_back(t);
- }
- }
- // multiscaleLaplace(_elements,ordered,coords);
- }
- // TEST --------------------------------------------------
-
for(unsigned int i = 0; i < ordered.size(); i++){
MVertex *v = ordered[i];
const double theta = 2 * M_PI * coords[i];
@@ -991,46 +975,39 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
myAssembler.numberVertex(t->getVertex(2), 0, 1);
}
}
- if (tom == CONVEXCOMBINATION){
- for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
- MTriangle *t = (*it2);
- myAssembler.numberVertex(t->getVertex(0), 0, 1);
- myAssembler.numberVertex(t->getVertex(1), 0, 1);
- myAssembler.numberVertex(t->getVertex(2), 0, 1);
- }
- }
-
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ MTriangle *t = (*it2);
+ myAssembler.numberVertex(t->getVertex(0), 0, 1);
+ myAssembler.numberVertex(t->getVertex(1), 0, 1);
+ myAssembler.numberVertex(t->getVertex(2), 0, 1);
+ }
+
Msg::Debug("Creating term %d dofs numbered %d fixed",
myAssembler.sizeOfR(), myAssembler.sizeOfF());
double t1 = Cpu();
- if (tom == CONVEXCOMBINATION){
- convexCombinationTerm laplace(model(), 1, &ONE);
- it = _compound.begin();
- for( ; it != _compound.end() ; ++it){
- for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- SElement se((*it)->triangles[i]);
- laplace.addToMatrix(myAssembler, &se);
- }
- }
- for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
- SElement se((*it2));
- laplace.addToMatrix(myAssembler, &se);
- }
+ femTerm<double> *mapping;
+ if (tom == HARMONIC){
+ mapping = new laplaceTerm(0, 1, &ONE);
}
- else {
- laplaceTerm laplace(model(), 1, &ONE);
- it = _compound.begin();
- for( ; it != _compound.end() ; ++it){
- for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- SElement se((*it)->triangles[i]);
- laplace.addToMatrix(myAssembler, &se);
- }
+ else if (tom == CONVEXCOMBINATION){
+ mapping = new convexCombinationTerm(0, 1, &ONE);
+ }
+
+ it = _compound.begin();
+ for( ; it != _compound.end() ; ++it){
+ for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ SElement se((*it)->triangles[i]);
+ mapping->addToMatrix(myAssembler, &se);
}
}
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ SElement se((*it2));
+ mapping->addToMatrix(myAssembler, &se);
+ }
double t2 = Cpu();
Msg::Debug("Assembly done in %8.3f seconds", t2 - t1);
@@ -1070,23 +1047,12 @@ void GFaceCompound::parametrize_conformal() const
MVertex *v1 = ordered[0];
MVertex *v2 = ordered[1];
-// if (!_interior_loops.empty()){
-// std::vector<MVertex*> ordered2;
-// bool success2 = orderVertices(_U1, ordered2, coords);
-// v2 = ordered2[ordered.size()/4];
-// }
-// else{
-// v2 = ordered[1];
-// //v2 = ordered[ordered.size()/2];
-// }
-
- //printf("Pinned vertex %g %g %g / %g %g %g \n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z());
- //exit(1);
-
myAssembler.fixVertex(v1, 0, 1, 0);
myAssembler.fixVertex(v1, 0, 2, 0);
myAssembler.fixVertex(v2, 0, 1, 1);
myAssembler.fixVertex(v2, 0, 2, 0);
+ //printf("Pinned vertex %g %g %g / %g %g %g \n", v1->x(), v1->y(), v1->z(), v2->x(), v2->y(), v2->z());
+ //exit(1);
std::list<GFace*>::const_iterator it = _compound.begin();
for( ; it != _compound.end(); ++it){
@@ -1099,7 +1065,17 @@ void GFaceCompound::parametrize_conformal() const
myAssembler.numberVertex(t->getVertex(1), 0, 2);
myAssembler.numberVertex(t->getVertex(2), 0, 2);
}
- }
+ }
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ MTriangle *t = (*it2);
+ myAssembler.numberVertex(t->getVertex(0), 0, 1);
+ myAssembler.numberVertex(t->getVertex(1), 0, 1);
+ myAssembler.numberVertex(t->getVertex(2), 0, 1);
+ myAssembler.numberVertex(t->getVertex(0), 0, 2);
+ myAssembler.numberVertex(t->getVertex(1), 0, 2);
+ myAssembler.numberVertex(t->getVertex(2), 0, 2);
+ }
+
simpleFunction<double> ONE(1.0);
simpleFunction<double> MONE(-1.0 );
@@ -1117,6 +1093,13 @@ void GFaceCompound::parametrize_conformal() const
cross21.addToMatrix(myAssembler, &se);
}
}
+ 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);
+ }
Msg::Debug("Assembly done");
_lsys->systemSolve();
@@ -1297,7 +1280,7 @@ GPoint GFaceCompound::point(double par1, double par2) const
// return lt->gf->point(pv.x(),pv.y());
// }
- const bool LINEARMESH = false; //false
+ const bool LINEARMESH = true; //false; //false
if(LINEARMESH){
@@ -1540,7 +1523,6 @@ void GFaceCompound::getTriangle(double u, double v,
void GFaceCompound::buildOct() const
{
- printStuff();
SBoundingBox3d bb;
int count = 0;
@@ -1601,26 +1583,47 @@ void GFaceCompound::buildOct() const
nbT = count;
Octree_Arrange(oct);
}
-//Verify topological conditions for computing the harmonic map
+
bool GFaceCompound::checkTopology() const
{
+
// FIXME!!! I think those things are wrong with cross-patch reparametrization
//if ((*(_compound.begin()))->geomType() != GEntity::DiscreteSurface)return true;
bool correctTopo = true;
-
+ if(allNodes.empty()) buildAllNodes();
+
int Nb = _interior_loops.size();
int G = genusGeom() ;
- if( G != 0 || Nb < 1) correctTopo = false;
+
+ double H = getSizeH();
+ double D = H;
+ if (_interior_loops.size() > 0)
+ D = getSizeBB(_U0);
+// for(std::list<std::list<GEdge*> >::const_iterator it = _interior_loops.begin(); it != _interior_loops.end(); it++){
+// double size = getSizeBB(*it);
+// D = std::min(D, size);
+// }
+ int AR = (int) ceil(H/D);
- if (!correctTopo){
- Msg::Warning("Wrong topology: Genus=%d and N boundary=%d", G, Nb);
- nbSplit = G+2;//std::max(G+2, 2);
- Msg::Info("Split surface %d in %d parts with Mesh partitioner", tag(), nbSplit);
+ if (G != 0 || Nb < 1){
+ correctTopo = false;
+ nbSplit = 2; //std::max(G+2, 2);
+ Msg::Info("-----------------------------------------------------------");
+ Msg::Warning("Wrong topology: Genus=%d, Nb boundaries=%d, AR=%g", G, Nb, H/D);
+ Msg::Info("*** Split surface %d in %d parts with Multilevel Mesh partitioner", tag(), nbSplit);
+ }
+ else if (G == 0 && AR > 3){
+ correctTopo = false;
+ nbSplit = -2;
+ Msg::Info("-----------------------------------------------------------");
+ Msg::Warning("Wrong topology: Genus=%d, Nb boundaries=%d, AR=%d", G, Nb, AR);
+ Msg::Info("*** Split surface %d in 2 parts with Laplacian Mesh partitioner", tag());
+ }
+ else{
+ Msg::Debug("Correct topology: Genus=%d and Nb boundaries=%d, AR=%g", G, Nb, H/D);
}
- else
- Msg::Info("Correct topology: Genus=%d and N boundary=%d", G, Nb);
-
+
return correctTopo;
}
@@ -1634,7 +1637,7 @@ bool GFaceCompound::checkAspectRatio() const
bool paramOK = true;
if(allNodes.empty()) buildAllNodes();
- double limit = 1.e10;
+ double limit = 1.e12;
double areaMax = 0.0;
int nb = 0;
std::list<GFace*>::const_iterator it = _compound.begin();
@@ -1662,20 +1665,19 @@ bool GFaceCompound::checkAspectRatio() const
}
}
- if (areaMax > limit && nb > 10) {
+ if (areaMax > limit && nb > 2) {
Msg::Warning("Geometrical aspect ratio too high (1/area_2D=%g)", areaMax);
SBoundingBox3d bboxH = bounds();
- double H = norm(SVector3(bboxH.max(), bboxH.min()));
+ double H = getSizeH();
double D = getSizeBB(_U0);
double eta = H/D;
- int split = (int)ceil(eta/3);
- nbSplit = std::max(split,2);
+ int split = -2;
//printf("H=%g, D=%g split=%d nbSplit=%d \n", H, D, split, nbSplit);
- Msg::Info("Partition geometry in N=%d parts", nbSplit);
+ Msg::Info("Partition geometry in N=%d parts", std::abs(nbSplit));
paramOK = false;
}
else {
- Msg::Info("Geometrical aspect ratio is OK :-)", areaMax);
+ Msg::Debug("Geometrical aspect ratio is OK :-)");
paramOK = true;
}
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 3e4f7ea3f709c49ac7a07b25b47af7f45af08991..64809a8e3e4b57cd0afbc0367b29e8cfea7b6e50 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -75,7 +75,6 @@ class GFaceCompound : public GFace {
void compute_distance() const;
bool checkOrientation(int iter) const;
void one2OneMap() const;
- bool checkCavity(std::vector<MElement*> &vTri) const;
bool checkAspectRatio() const;
void computeNormals () const;
void getBoundingEdges();
@@ -87,6 +86,7 @@ class GFaceCompound : public GFace {
void printStuff() const;
bool trivial() const ;
linearSystem <double> *_lsys;
+ double getSizeH() const;
double getSizeBB(const std::list<GEdge* > &elist) const;
SBoundingBox3d boundEdges(const std::list<GEdge* > &elist) const;
SOrientedBoundingBox obb_boundEdges(const std::list<GEdge* > &elist) const;
@@ -101,6 +101,7 @@ class GFaceCompound : public GFace {
Range<double> parBounds(int i) const
{ return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
virtual GPoint point(double par1, double par2) const;
+ typeOfMapping getTypeOfMapping() { return _mapping;}
SPoint2 parFromPoint(const SPoint3 &p) const;
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
@@ -113,10 +114,8 @@ class GFaceCompound : public GFace {
virtual bool checkTopology() const;
bool parametrize() const ;
void coherenceNormals();
- void partitionFaceCM();
virtual std::list<GFace*> getCompounds() const {return _compound;};
mutable int nbSplit;
- mutable bool checked;
private:
typeOfIsomorphism _type;
typeOfMapping _mapping;
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 25c1cf0e75ba90de2db41b81a35b8d24b204535d..6949444b6940c2d06fe3b09f67b7909a58e4ed30 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1072,6 +1072,8 @@ void GModel::createTopologyFromMesh()
if((*it)->geomType() == GEntity::DiscreteVolume)
discRegions.push_back((discreteRegion*) *it);
+ //EMI-FIX in case of createTopology for Volumes
+ //all faces are set to the volume
for (std::vector<discreteRegion*>::iterator it = discRegions.begin();
it != discRegions.end(); it++)
(*it)->setBoundFaces();
@@ -1137,10 +1139,14 @@ void GModel::createTopologyFromMesh()
//create Topo From Faces
createTopologyFromFaces(discFaces);
+ //create
+ exportDiscreteGEOInternals();
+
}
void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
{
+
std::vector<discreteEdge*> discEdges;
for(eiter it = firstEdge(); it != lastEdge(); it++){
@@ -1154,7 +1160,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
for (std::vector<discreteFace*>::iterator it = discFaces.begin(); it != discFaces.end(); it++)
(*it)->findEdges(map_edges);
- //return if no boundary edges (torus, sphere, ..)
+ //return if no boundary edges (torus, sphere, ...)
if (map_edges.empty()) return;
// create reverse map, for each face find set of MEdges that are
diff --git a/Geo/GModel.h b/Geo/GModel.h
index b0c17356829fbe64da9f60b0538dd69d6089c7d4..b43b3d33d6e492da559f76e6d7566fdf6e34782b 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -351,6 +351,7 @@ class GModel
static int readGEO(const std::string &name);
int writeGEO(const std::string &name, bool printLabels=true);
int importGEOInternals();
+ int exportDiscreteGEOInternals();
// Fourier model
int readFourier();
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index d1fe8ce0f3dcfe73b62c28c5c478a20bf8a0cdbf..465a1ce8c69fffd8b0ba70a9d639025655d6bbb0 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -41,8 +41,76 @@ int GModel::readGEO(const std::string &name)
return GModel::current()->importGEOInternals();
}
+int GModel::exportDiscreteGEOInternals()
+{
+
+ _geo_internals = new GEO_Internals;
+
+ for(viter it = firstVertex(); it != lastVertex(); it++){
+ Vertex *v = Create_Vertex((*it)->tag(), (*it)->x(), (*it)->y(), (*it)->z(), (*it)->prescribedMeshSizeAtVertex(), 1.0);
+ Tree_Add(GModel::current()->getGEOInternals()->Points, &v);
+ }
+
+ for(eiter it = firstEdge(); it != lastEdge(); it++){
+ if((*it)->geomType() == GEntity::DiscreteCurve){
+ Curve *c = Create_Curve((*it)->tag(), MSH_SEGM_DISCRETE, 1, NULL, NULL, -1, -1, 0., 1.);
+ List_T *points = Tree2List(_geo_internals->Points);
+ GVertex *gvb = (*it)->getBeginVertex();
+ GVertex *gve = (*it)->getEndVertex();
+ c->Control_Points = List_Create(2, 1, sizeof(Vertex *));
+ for(int i = 0; i < List_Nbr(points); i++) {
+ Vertex *v;
+ List_Read(points, i, &v);
+ if (v->Num == gvb->tag()) {
+ List_Add(c->Control_Points, &v);
+ c->beg = v;
+ }
+ if (v->Num == gve->tag()) {
+ List_Add(c->Control_Points, &v);
+ c->end = v;
+ }
+ }
+ Tree_Add(GModel::current()->getGEOInternals()->Curves, &c);
+ CreateReversedCurve(c);
+ }
+ }
+
+ for(fiter it = firstFace(); it != lastFace(); it++){
+ if((*it)->geomType() == GEntity::DiscreteSurface){
+ Surface *s = Create_Surface((*it)->tag(), MSH_SURF_DISCRETE);
+ std::list<GEdge*> edges = (*it)->edges();
+ s->Generatrices = List_Create(edges.size(), 1, sizeof(Curve *));
+ List_T *curves = Tree2List(_geo_internals->Curves);
+ Curve *c;
+ for(std::list<GEdge*>::iterator ite = edges.begin(); ite != edges.end(); ite++){
+ for(int i = 0; i < List_Nbr(curves); i++) {
+ List_Read(curves, i, &c);
+ if (c->Num == (*ite)->tag()) {
+ List_Add(s->Generatrices, &c);
+ }
+ }
+ }
+ Tree_Add(GModel::current()->getGEOInternals()->Surfaces, &s);
+ }
+ }
+
+ //create Volumes from discreteRegions
+ //TODO
+
+ Msg::Debug("Geo internal model has:");
+ List_T *points = Tree2List(_geo_internals->Points);
+ List_T *curves = Tree2List(_geo_internals->Curves);
+ List_T *surfaces = Tree2List(_geo_internals->Surfaces);
+ Msg::Debug("%d Vertices", List_Nbr(points));
+ Msg::Debug("%d Edges", List_Nbr(curves));
+ Msg::Debug("%d Faces", List_Nbr(surfaces));
+
+ return 1;
+}
+
int GModel::importGEOInternals()
{
+
if(Tree_Nbr(_geo_internals->Points)) {
List_T *points = Tree2List(_geo_internals->Points);
for(int i = 0; i < List_Nbr(points); i++){
@@ -193,7 +261,7 @@ int GModel::importGEOInternals()
}
- Msg::Debug("Gmsh model imported:");
+ Msg::Debug("Gmsh model (GModel) imported:");
Msg::Debug("%d Vertices", vertices.size());
Msg::Debug("%d Edges", edges.size());
Msg::Debug("%d Faces", faces.size());
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 79142edf79e281b76fe6be3ebd632c4e6eb31669..f0643c4e6528abb8c115c075e46f1640773c345f 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2352,7 +2352,6 @@ static int Extrude_ProtudeSurface(int type, int is,
Msg::Debug("Extrude Surface %d", is);
chapeau = DuplicateSurface(ps, false);
-
chapeau->Extrude = new ExtrudeParams(COPIED_ENTITY);
chapeau->Extrude->fill(type, T0, T1, T2, A0, A1, A2, X0, X1, X2, alpha);
chapeau->Extrude->geo.Source = is; // not ps->Num: we need the sign info
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 0b1aadbdd7dff1aa9238b8c6dd4c16681c6c4817..c5a2de40e86730616e50d0dba0fe412fc79d6d34 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -418,7 +418,7 @@ GPoint discreteEdge::point(double par) const
MVertex *vB = lines[iEdge]->getVertex(0);
MVertex *vE = lines[iEdge]->getVertex(1);
- const bool LINEARMESH = false; //false; //false;
+ const bool LINEARMESH = true; //false; //false;
if (LINEARMESH){
//linear Lagrange mesh
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 09d1ca0abbb67b9b74469cb02f0efdc11b40ab4d..4d2f3fde2ee2040ee621a9939f4c8dca633d2723 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -35,6 +35,7 @@
#include "meshPartition.h"
#include "CreateFile.h"
#include "Context.h"
+#include "multiscalePartition.h"
void fourthPoint(double *p1, double *p2, double *p3, double *p4)
{
@@ -1285,8 +1286,6 @@ void meshGFace::operator() (GFace *gf)
Msg::Error("Impossible to mesh face %d", gf->tag());
}
- // gmshQMorph(gf);
-
Msg::Debug("Type %d %d triangles generated, %d internal vertices",
gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size());
}
@@ -1295,7 +1294,7 @@ bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges)
{
bool isMeshed = false;
#if defined(HAVE_SOLVER)
- //Check Topology
+
bool correctTopo = gf->checkTopology();
if (!correctTopo){
partitionAndRemesh((GFaceCompound*) gf);
@@ -1303,12 +1302,12 @@ bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges)
return isMeshed;
}
- //Parametrize
bool correctParam = gf->parametrize();
+
if (!correctParam){
- partitionAndRemesh((GFaceCompound*) gf);
- isMeshed = true;
- return isMeshed;
+ partitionAndRemesh((GFaceCompound*) gf);
+ isMeshed = true;
+ return isMeshed;
}
//Replace edges by their compounds
@@ -1316,11 +1315,9 @@ bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges)
std::list<GEdge*>::iterator it = edges.begin();
while(it != edges.end()){
if((*it)->getCompound()){
- //printf("replace edge %d by %d\n", (*it)->tag(), (*it)->getCompound()->tag() );
mySet.insert((*it)->getCompound());
}
else{
- //printf("insert edge =%d \n", (*it)->tag());
mySet.insert(*it);
}
++it;
@@ -1338,10 +1335,21 @@ void partitionAndRemesh(GFaceCompound *gf)
//Partition the mesh and createTopology for new faces
//-----------------------------------------------------
- int NF = gf->nbSplit;
std::list<GFace*> cFaces = gf->getCompounds();
- PartitionZeroGenus(cFaces, NF);
+ //PartitionMeshFace(cFaces, options);
+
+ std::vector<MElement *> elements;
+ for (std::list<GFace*>::iterator it = cFaces.begin(); it != cFaces.end(); it++)
+ for(unsigned int j = 0; j < (*it)->getNumMeshElements(); j++)
+ elements.push_back((*it)->getMeshElement(j));
+
+ typeOfPartition method;
+ if(gf->nbSplit > 0) method = MULTILEVEL;
+ else method = LAPLACIAN;
+
+ multiscalePartition *msp = new multiscalePartition(elements, abs(gf->nbSplit), method);
+ int NF = msp->getNumberOfParts();
int numv = gf->model()->maxVertexNum() + 1;
int nume = gf->model()->maxEdgeNum() + 1;
int numf = gf->model()->maxFaceNum() + 1;
@@ -1350,18 +1358,13 @@ void partitionAndRemesh(GFaceCompound *gf)
gf->model()->createTopologyFromFaces(pFaces);
- Msg::Info("-----------------------------------------------------------");
- Msg::Info("Multiscale Partition SUCCESSFULLY PERFORMED : %d parts", NF);
+ Msg::Info("*** Multiscale Partition SUCCESSFULLY PERFORMED : %d parts", NF );
CreateOutputFile("multiscalePARTS.msh", CTX::instance()->mesh.format);
- Msg::Info("-----------------------------------------------------------");
- //Remesh new faces (Compound Lines and Compound Surfaces)
+ //Remesh new faces (Compound Lines and Compound Surfaces)
//-----------------------------------------------------
-
- Msg::Info("-----------------------------------------------------------");
- Msg::Info("Parametrize Compounds");
- Msg::Info("-----------------------------------------------------------");
-
+ Msg::Info("*** Parametrize Compounds:");
+
//Parametrize Compound Lines
int NE = gf->model()->maxEdgeNum() - nume + 1;
for (int i=0; i < NE; i++){
@@ -1385,15 +1388,14 @@ void partitionAndRemesh(GFaceCompound *gf)
int num_gfc = numf + NF + i ;
f_compound.push_back(pf);
Msg::Info("Parametrize Compound Surface (%d) = %d discrete face", num_gfc, pf->tag() );
- GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound, b[0], b[1], b[2], b[3]);
+ GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound,
+ b[0], b[1], b[2], b[3], 0, gf->getTypeOfMapping() );
gf->model()->add(gfc);
gfc->parametrize();
}
- Msg::Info("-----------------------------------------------------------");
- Msg::Info("Mesh Compounds");
- Msg::Info("-----------------------------------------------------------");
+ Msg::Info("*** Mesh Compounds:");
//Mesh 1D and 2D
for (int i=0; i < NE; i++){
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index e4e1d03d4313ca0f0d951912056dbd5093a39b32..ff9821f0a8436bf6bd9480cd82794b761bea9890 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -27,7 +27,6 @@
#include "discreteEdge.h"
#include "discreteFace.h"
#include "GFaceCompound.h"
-#include "multiscalePartition.h"
#include "Context.h"
@@ -103,28 +102,7 @@ void MakeGraphDIM(const EntIter begin, const EntIter end,
*
******************************************************************************/
-bool PartitionZeroGenus(std::list<GFace*> &cFaces, int &nbParts){
-
- meshPartitionOptions options;
- options = CTX::instance()->partitionOptions;
- options.num_partitions = nbParts;
- options.partitioner = 1;//1 CHACO //2 METIS
- if ( options.partitioner == 1){
- options.global_method = 2;// 1 Multilevel-KL 2 Spectral
- options.mesh_dims[0] = nbParts;
- }
-
-// PartitionMeshFace(cFaces, options);
-
- std::vector<MElement *> elements;
- for (std::list<GFace*>::iterator it = cFaces.begin(); it != cFaces.end(); it++)
- for(unsigned int j = 0; j < (*it)->getNumMeshElements(); j++)
- elements.push_back((*it)->getMeshElement(j));
-
- multiscalePartition *msp = new multiscalePartition(elements, options);
- nbParts = msp->getNumberOfParts();
- return true;
-}
+
int PartitionMeshElements( std::vector<MElement*> &elements, meshPartitionOptions &options){
diff --git a/Mesh/meshPartition.h b/Mesh/meshPartition.h
index b64aed669fac2041f453b284293a3180930df665..9e7780981a1dc0e13784c3ea064f27e2d0799c95 100644
--- a/Mesh/meshPartition.h
+++ b/Mesh/meshPartition.h
@@ -18,6 +18,7 @@ class GFace;
class Graph;
typedef std::vector<BoElemGr> BoElemGrVec;
+typedef enum {LAPLACIAN= 0, MULTILEVEL=1} typeOfPartition;
/*******************************************************************************
*
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index dda34a8a845659a29f94ce26cafdab678cfe909c..8cf1826f0fc394cbf8781c28951e48cd4388b622 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -4,6 +4,8 @@
#include "meshPartition.h"
#include "MEdge.h"
#include "MElement.h"
+#include "multiscaleLaplace.h"
+#include "Context.h"
static void recur_connect (MVertex *v,
std::multimap<MVertex*,MEdge> &v2e,
@@ -113,38 +115,33 @@ static int getAspectRatio (std::vector<MElement *> &elements,
//double H = obbox.getMaxSize();
double H = norm(SVector3(bb.max(), bb.min()));
- double D = 0.0;
- if (boundaries.size() == 0){
- D=1.;
- }
- else{
- for (unsigned i=0; i< boundaries.size(); i++){
- std::set<MVertex*> vb;
- std::vector<MEdge> iBound = boundaries[i];
- for (unsigned j=0; j< iBound.size(); j++){
- MEdge e = iBound[j];
- vb.insert(e.getVertex(0));
- vb.insert(e.getVertex(1));
- }
- std::vector<SPoint3> vBounds;
- for (std::set<MVertex* >::iterator it = vb.begin(); it != vb.end(); it++){
- SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
- vBounds.push_back(pt);
- }
- SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds);
- D = std::max(D, obboxD.getMaxSize());
- }
- }
- int eta = (int)ceil(H/D);
+ double D = H;
+ for (unsigned i=0; i< boundaries.size(); i++){
+ std::set<MVertex*> vb;
+ std::vector<MEdge> iBound = boundaries[i];
+ for (unsigned j=0; j< iBound.size(); j++){
+ MEdge e = iBound[j];
+ vb.insert(e.getVertex(0));
+ vb.insert(e.getVertex(1));
+ }
+ std::vector<SPoint3> vBounds;
+ for (std::set<MVertex* >::iterator it = vb.begin(); it != vb.end(); it++){
+ SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
+ vBounds.push_back(pt);
+ }
+ SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds);
+ D = std::min(D, obboxD.getMaxSize());
+ }
+ int AR = (int)ceil(H/D);
- return eta;
+ return AR;
}
-static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int &eta){
+static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int &AR){
std::vector<std::vector<MEdge> > boundaries;
genus = getGenus(elements, boundaries);
- eta = getAspectRatio(elements, boundaries);
+ AR = getAspectRatio(elements, boundaries);
return;
@@ -162,11 +159,54 @@ static void partitionRegions (std::vector<MElement*> &elements,
}
+static void printLevel ( std::vector<MElement *> &elements, int recur, int region){
+
+ char fn[256];
+ sprintf(fn, "part_%d_%d.msh", recur, region );
+ double version = 2.0;
+
+ std::set<MVertex*> vs;
+ for (int i=0;i<elements.size();i++){
+ for (int j=0;j<elements[i]->getNumVertices();j++){
+ vs.insert(elements[i]->getVertex(j));
+ }
+ }
+
+ bool binary = false;
+ FILE *fp = fopen (fn, "w");
+ fprintf(fp, "$MeshFormat\n");
+ fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
+ fprintf(fp, "$EndMeshFormat\n");
+
+ fprintf(fp,"$Nodes\n%d\n",vs.size());
+ std::set<MVertex*> :: iterator it = vs.begin();
+ int index = 1;
+ for (; it != vs.end() ; ++it){
+ (*it)->setIndex(index++);
+ fprintf(fp,"%d %g %g %g\n",(*it)->getIndex(),
+ (*it)->x(),(*it)->y(),(*it)->z());
+ }
+ fprintf(fp,"$EndNodes\n",elements.size());
+
+ fprintf(fp,"$Elements\n%d\n",elements.size());
+ for (int i=0;i<elements.size();i++){
+ elements[i]->writeMSH(fp,version);
+ }
+ fprintf(fp,"$EndElements\n%d\n",elements.size());
+
+ fclose(fp);
+}
-multiscalePartition::multiscalePartition (std::vector<MElement *> &elements,
- meshPartitionOptions &_options){
+multiscalePartition::multiscalePartition (std::vector<MElement *> &elements, int nbParts, typeOfPartition method)
+{
- options = _options;
+ options = CTX::instance()->partitionOptions;
+ options.num_partitions = nbParts;
+ options.partitioner = 1; //1 CHACO, 2 METIS
+ if ( options.partitioner == 1){
+ options.global_method = 2;// 1 Multilevel-KL, 2 Spectral
+ options.mesh_dims[0] = nbParts;
+ }
partitionLevel *level = new partitionLevel;
level->elements.insert(level->elements.begin(),elements.begin(),elements.end());
@@ -175,7 +215,7 @@ multiscalePartition::multiscalePartition (std::vector<MElement *> &elements,
levels.push_back(level);
- partition(*level);
+ partition(*level, nbParts, method);
totalParts = assembleAllPartitions();
@@ -188,11 +228,19 @@ void multiscalePartition:: setNumberOfPartitions(int &nbParts)
}
}
-void multiscalePartition::partition(partitionLevel & level){
+void multiscalePartition::partition(partitionLevel & level, int nbParts, typeOfPartition method){
+
#if defined(HAVE_METIS) || defined(HAVE_CHACO)
- PartitionMeshElements(level.elements, options);
- std::vector<std::vector<MElement*> > regions(options.num_partitions);
+ if (method == LAPLACIAN){
+ multiscaleLaplace multiLaplace(level.elements, level.recur);
+ }
+ else if (method == MULTILEVEL){
+ setNumberOfPartitions(nbParts);
+ PartitionMeshElements(level.elements, options);
+ }
+
+ std::vector<std::vector<MElement*> > regions(nbParts);
partitionRegions(level.elements, regions);
level.elements.clear();
@@ -204,23 +252,31 @@ void multiscalePartition::partition(partitionLevel & level){
nextLevel->region = i;
levels.push_back(nextLevel);
- int genus, eta;
- getGenusAndRatio(regions[i], genus, eta);
+ int genus, AR;
+ getGenusAndRatio(regions[i], genus, AR);
+
if (genus < 0) {
Msg::Error("Genus partition is negative G=%d!", genus);
exit(1);
}
- if (genus != 0 || eta > 4 ){
- int nbParts = std::max(genus+2, (int)eta/4);
- setNumberOfPartitions(nbParts);
- Msg::Info("Multiscale partition, level %d region %d is %d-GENUS (ETA=%d) ---> partition %d parts",
- nextLevel->recur,nextLevel->region, genus, eta, nbParts);
-
- partition(*nextLevel);
+
+ printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
+
+ if (genus != 0 ){
+ int nbParts = 2; //std::max(genus+2,2);
+ Msg::Info("Multiscale part: level %d region %d is %d-GENUS (AR=%d) ---> MULTILEVEL partition %d parts",
+ nextLevel->recur,nextLevel->region, genus, AR, nbParts);
+ partition(*nextLevel, nbParts, MULTILEVEL);
}
+ else if (genus == 0 && AR > 3 ){
+ int nbParts = 2;
+ Msg::Info("Multiscale part: level %d region %d is ZERO-GENUS (AR=%d) ---> LAPLACIAN partition %d parts",
+ nextLevel->recur,nextLevel->region, AR, nbParts);
+ partition(*nextLevel, nbParts, LAPLACIAN);
+ }
else {
- Msg::Info("Multiscale Partition, level %d, region %d is ZERO-GENUS (ETA=%d)",
- nextLevel->recur,nextLevel->region, eta);
+ Msg::Info("Multiscale part: level %d, region %d is ZERO-GENUS (AR=%d)",
+ nextLevel->recur,nextLevel->region, AR);
}
}
diff --git a/Mesh/multiscalePartition.h b/Mesh/multiscalePartition.h
index 1014005814ab477a475b471178873807e2f0681e..3c071685426d444f2fd3eec523546b7e47a70b25 100644
--- a/Mesh/multiscalePartition.h
+++ b/Mesh/multiscalePartition.h
@@ -10,6 +10,7 @@
#include <map>
#include "linearSystemGMM.h"
#include "meshPartitionOptions.h"
+#include "meshPartition.h"
class MElement;
struct meshPartitionOptions;
@@ -24,13 +25,12 @@ class multiscalePartition{
private:
std::vector<partitionLevel*> levels;
- void partition(partitionLevel &level);
+ void partition(partitionLevel &level, int nbParts, typeOfPartition method);
int totalParts;
meshPartitionOptions options;
public:
- multiscalePartition(std::vector<MElement *> &elements,
- meshPartitionOptions &options);
+ multiscalePartition(std::vector<MElement *> &elements, int nbParts, typeOfPartition method);
int assembleAllPartitions();
void setNumberOfPartitions(int &nbParts);
int getNumberOfParts(){return totalParts;}
diff --git a/Solver/TESTCASES/AdvectionDiffusion.lua b/Solver/TESTCASES/AdvectionDiffusion.lua
index 68cb64479b5ddfa26bccf2b39fc01fd598b87269..847b2e9e2ea8be66e4756bf7ea9ed139cebd3f1f 100644
--- a/Solver/TESTCASES/AdvectionDiffusion.lua
+++ b/Solver/TESTCASES/AdvectionDiffusion.lua
@@ -1,6 +1,6 @@
model = GModel ()
model:load ('square.geo')
-model:load ('square.msh')
+model:mesh(2)
dg = dgSystemOfEquations (model)
dg:setOrder(5)
@@ -40,9 +40,9 @@ dg:L2Projection(createFunction.lua(1,'initial_condition','XYZ'))
dg:exportSolution('output/Advection_00000')
-- main loop
-for i=1,10000 do
+for i=1,10 do
norm = dg:RK44(0.001)
- if (i % 50 == 0) then
+ if (i % 1 == 0) then
print('iter',i,norm)
dg:exportSolution(string.format("output/Advection-%05d", i))
end
diff --git a/Solver/TESTCASES/Stommel.lua b/Solver/TESTCASES/Stommel.lua
index bc3c89380dabf3c86b86db405f0d251dae8d9d13..3081a4b5ca300549eeb611698a96bf4023437ce8 100644
--- a/Solver/TESTCASES/Stommel.lua
+++ b/Solver/TESTCASES/Stommel.lua
@@ -2,8 +2,12 @@
order = 3
model = GModel()
model:load ('stommel_square.msh')
+print('load..')
+
dg = dgSystemOfEquations (model)
dg:setOrder (order)
+
+
dg:setConservationLaw('ShallowWater2d')
dg:addBoundaryCondition('Wall','Wall')
dg:setup()
diff --git a/Solver/TESTCASES/WavePulse.lua b/Solver/TESTCASES/WavePulse.lua
index 0d5c1e189a8978cfddfb3439c4273b7f55021c6f..02dc142988641a9a80ad52494485aa8b004ec155 100644
--- a/Solver/TESTCASES/WavePulse.lua
+++ b/Solver/TESTCASES/WavePulse.lua
@@ -22,12 +22,13 @@ end
order = 5
print'*** Loading the mesh and the model ***'
myModel = GModel ()
-myModel:load ('square.geo')
-myModel:load ('square.msh')
+myModel:load ('box.geo')
+myModel:load ('box.msh')
print'*** Create a dg solver ***'
DG = dgSystemOfEquations (myModel)
DG:setOrder(order)
DG:setConservationLaw('WaveEquation')
+
DG:addBoundaryCondition('Border','Wall')
DG:setup()
diff --git a/Solver/TESTCASES/box.geo b/Solver/TESTCASES/box.geo
new file mode 100644
index 0000000000000000000000000000000000000000..62fd4b8e50c641a8f97b5644bcff41e57043b547
--- /dev/null
+++ b/Solver/TESTCASES/box.geo
@@ -0,0 +1,24 @@
+lc = 1./2.5;
+Point(1) = { 0.0 ,0.0 ,0., lc};
+Point(2) = { 1.0 ,0.0 ,0., lc};
+Point(3) = { 1.0 ,1.0 ,0., lc};
+Point(4) = { 0.0 ,1.0 ,0., lc};
+
+// Point(1) = { 0., 0.0 ,0.0 ,lc};
+// Point(2) = { 0., 1.0 ,0.0 , lc};
+// Point(3) = { 0., 1.0 ,1.0 , lc};
+// Point(4) = { 0., 0.0 ,1.0 , lc};
+
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+
+Line Loop(4) = {1,2,3,4};
+Plane Surface(5) = {4};
+
+Extrude {0, 0, 0.1} {
+ Surface{5};
+}
+Surface Loop(28) = {5, 14, 18, 22, 26, 27};
+Volume(29) = {28};
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index b11d0812aeb2ed6bed85e63a64258d49828a2427..031b69a463ab58cf3a40ffb43efe4bc28fa46b87 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -8,7 +8,11 @@
#include "SPoint3.h"
#include "dofManager.h"
#include "laplaceTerm.h"
+#include "convexCombinationTerm.h"
#include "linearSystemCSR.h"
+#include "robustPredicates.h"
+#include "meshGFaceOptimize.h"
+#include "GFaceCompound.h"
#ifdef HAVE_GMM
#include "linearSystemGMM.h"
@@ -18,18 +22,7 @@
#include "MTriangle.h"
#include "robustPredicates.h"
-/*
-//-- returns 0 if no intersection occurs --\\
-//-- returns 1
-p1 (1-u) + p2 u = q1 (1-t) + q2 t
-
-we search the intersection into segment q1-q2
-
-robustPredicates
-
-if ( isLeftOf (
-
-*/
+//--------------------------------------------------------------
static int intersection_segments_b (SPoint2 &p1, SPoint2 &p2,
SPoint2 &q1, SPoint2 &q2,
double x[2]){
@@ -45,8 +38,159 @@ static int intersection_segments_b (SPoint2 &p1, SPoint2 &p2,
double QP2 = robustPredicates::orient2d(Q1,Q2,P2);
}
+//--------------------------------------------------------------
+static void recur_connect (MVertex *v,
+ std::multimap<MVertex*,MEdge> &v2e,
+ std::set<MEdge,Less_Edge> &group,
+ std::set<MVertex*> &touched){
+
+ if (touched.find(v) != touched.end())return;
+
+ touched.insert(v);
+ for (std::multimap <MVertex*,MEdge>::iterator it = v2e.lower_bound(v);
+ it != v2e.upper_bound(v) ; ++it){
+ group.insert(it->second);
+ for (int i=0;i<it->second.getNumVertices();++i){
+ recur_connect (it->second.getVertex(i),v2e,group,touched);
+ }
+ }
+}
+//--------------------------------------------------------------
+static void connected_bounds (std::vector<MElement*> &elements,
+ std::vector<std::vector<MEdge> > &boundaries){
+
+ std::vector<MEdge> bEdges;
+ for(unsigned int i = 0; i < elements.size(); i++){
+ for(int j = 0; j < elements[i]->getNumEdges(); j++){
+ MEdge me = elements[i]->getEdge(j);
+ if(std::find(bEdges.begin(), bEdges.end(), me) == bEdges.end())
+ bEdges.push_back(me);
+ else
+ bEdges.erase(std::find(bEdges.begin(), bEdges.end(),me));
+ }
+ }
+
+ std::multimap<MVertex*,MEdge> v2e;
+ for (unsigned i=0;i<bEdges.size();++i){
+ for (int j=0;j<bEdges[i].getNumVertices();j++){
+ v2e.insert(std::make_pair(bEdges[i].getVertex(j),bEdges[i]));
+ }
+ }
+ while (!v2e.empty()){
+ std::set<MEdge, Less_Edge> group;
+ std::set<MVertex*> touched;
+ recur_connect (v2e.begin()->first,v2e,group,touched);
+ std::vector<MEdge> temp;
+ temp.insert(temp.begin(), group.begin(), group.end());
+ boundaries.push_back(temp);
+ for ( std::set<MVertex*>::iterator it = touched.begin() ; it != touched.end();++it)
+ v2e.erase(*it);
+ }
+ return;
+}
+//--------------------------------------------------------------
+static double getSizeBB(std::vector<MEdge> &me){
+
+ SBoundingBox3d bb ;
+ SOrientedBoundingBox obbox ;
+
+ std::vector<SPoint3> vertices;
+ for (int i=0; i< me.size(); i++){
+ MVertex *v0 = me[i].getVertex(0);
+ MVertex *v1 = me[i].getVertex(1);
+ SPoint3 pt1(v0->x(),v0->y(), v0->z());
+ vertices.push_back(pt1);
+ SPoint3 pt2(v1->x(),v1->y(), v1->z());
+ vertices.push_back(pt2);
+ bb+=pt1;
+ bb+=pt2;
+ }
+
+ //double H = norm(SVector3(bb.max(), bb.min()));
+ //printf("H=%g \n", H);
+
+ obbox = SOrientedBoundingBox::buildOBB(vertices);
+ double H = obbox.getMaxSize();
+
+ return H;
+
+}
+//--------------------------------------------------------------
+static void ordering_dirichlet(std::vector<MElement*> &elements,
+ std::vector<std::pair<MVertex*,double> > &dirichletNodes){
+
+ //finding all boundaries
+ std::vector<std::vector<MEdge> > boundaries;
+ connected_bounds(elements,boundaries);
+
+ //largest boundary is dirichlet boundary
+ std::vector<MEdge> dirichletEdges;
+ double maxSize = 0.0;
+ for(int i=0; i< boundaries.size(); i++){
+ std::vector<MEdge> iBound = boundaries[i];
+ double size = getSizeBB(iBound);
+ if (size > maxSize) {
+ dirichletEdges = iBound;
+ maxSize = size;
+ }
+ }
+
+ //ordering dirichletNodes
+ dirichletNodes.clear();
+ std::list<MEdge> temp;
+ double tot_length = 0.0;
+ for(int i= 0 ; i < dirichletEdges.size(); i++ ){
+ MVertex *v0 = dirichletEdges[i].getVertex(0);
+ MVertex *v1 = dirichletEdges[i].getVertex(1);
+ double len = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z())) ;
+ tot_length += sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()));
+ temp.push_back(dirichletEdges[i]);
+ }
+
+ dirichletNodes.push_back(std::make_pair(dirichletEdges[0].getVertex(0),0.0));
+ MVertex *current_v = dirichletEdges[0].getVertex(1);
+ temp.erase(temp.begin());
+
+ while(temp.size()){
+ bool found = false;
+ for(std::list<MEdge>::iterator itl = temp.begin(); itl != temp.end(); ++itl){
+ MVertex *v0 = itl->getVertex(0);
+ MVertex *v1 = itl->getVertex(1);
+ if(v0 == current_v){
+ found = true;
+ current_v = v1;
+ temp.erase(itl);
+ double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()));
+ double iLength = dirichletNodes[dirichletNodes.size()-1].second + (length / tot_length);
+ dirichletNodes.push_back(std::make_pair(current_v,iLength));
+ break;
+ }
+ else if(v1 == current_v){
+ found = true;
+ current_v = v0;
+ temp.erase(itl);
+ double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()));
+ double iLength = dirichletNodes[dirichletNodes.size()-1].second + (length / tot_length);
+ dirichletNodes.push_back(std::make_pair(current_v,iLength));
+ break;
+ }
+ }
+ if(!found) return ;
+ }
+
+ return;
+}
+//--------------------------------------------------------------
static int intersection_segments (SPoint2 &p1, SPoint2 &p2,
SPoint2 &q1, SPoint2 &q2,
double x[2]){
@@ -62,20 +206,22 @@ static int intersection_segments (SPoint2 &p1, SPoint2 &p2,
x[1] >= 0.0 && x[1] <= 1.);
}
-
+//--------------------------------------------------------------
static void recur_compute_centers_ (double R, double a1, double a2,
multiscaleLaplaceLevel * root ){
+ root->radius = R;
std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers = root->cut;
+ centers.clear();
multiscaleLaplaceLevel* zero = 0;
SPoint2 PL (R*cos(a1),R*sin(a1));
SPoint2 PR (R*cos(a2),R*sin(a2));
centers.push_back(std::make_pair(PL,zero));
centers.push_back(std::make_pair(PR,zero));
- for (int i=0;i<root->childeren.size();i++){
- centers.push_back(std::make_pair(root->childeren[i]->center,root->childeren[i]));
- multiscaleLaplaceLevel* m = root->childeren[i];
+ for (int i=0;i<root->children.size();i++){
+ multiscaleLaplaceLevel* m = root->children[i];
+ centers.push_back(std::make_pair(m->center,m));
m->radius = 0.0;
for (std::map<MVertex*,SPoint2>::iterator it = m->coordinates.begin();
it != m->coordinates.end() ; ++it){
@@ -84,6 +230,35 @@ static void recur_compute_centers_ (double R, double a1, double a2,
(m->center.y() - p.y())*(m->center.y() - p.y())));
}
}
+
+ //compute interior_loops
+ std::vector<std::vector<MEdge> > boundaries;
+ connected_bounds(root->elements, boundaries);
+
+ //add the center of real holes ...
+ if (root->children.size()==0 || boundaries.size()-1 != root->children.size()){
+ for (int i = 0; i < boundaries.size(); i++){
+ std::vector<MEdge> me = boundaries[i];
+ SPoint2 c(0.0,0.0);
+ double rad = 0.0;
+ for(int j= 0; j< me.size(); j++){
+ MVertex *v = me[j].getVertex(0);
+ std::map<MVertex *, SPoint2>::iterator it0 = root->coordinates.find(v);
+ c += it0->second;
+ }
+ c *= 1./((double)me.size());
+ for(int j= 0; j< me.size(); j++){
+ SPoint2 p = root->coordinates[me[j].getVertex(0)];
+ rad = std::max(rad,sqrt ((c.x() - p.x())*(c.x() - p.x())+
+ (c.y() - p.y())*(c.y() - p.y())));
+ }
+ if (std::abs(rad/root->radius) < 0.8 && abs(rad) < 0.99){
+ centers.push_back(std::make_pair(c,zero));
+ }
+ }
+ }
+
+ //sort centers
std::sort(centers.begin(),centers.end());
for (int i=1;i<centers.size()-1;i++){
@@ -96,8 +271,10 @@ static void recur_compute_centers_ (double R, double a1, double a2,
recur_compute_centers_ (m2->radius, a1, a2, m2);
}
}
-}
+
+}
+//--------------------------------------------------------------
static void recur_cut_edges_ (multiscaleLaplaceLevel * root,
std::multimap<MEdge,MElement*,Less_Edge> &e2e,
std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
@@ -145,12 +322,13 @@ static void recur_cut_edges_ (multiscaleLaplaceLevel * root,
}
}
}
-
+//--------------------------------------------------------------
static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
std::set<MVertex*> &cutVertices,
std::set<MEdge,Less_Edge> &theCut,
std::vector<MElement*> &_all){
+
std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers = root->cut;
std::vector<MElement*> newElements;
for (int i=0;i<root->elements.size();i++){
@@ -231,7 +409,7 @@ static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
}
}
}
-
+//--------------------------------------------------------------
static void recur_split_ (MElement *e,
std::multimap<MEdge,MElement*,Less_Edge> &e2e,
std::set<MElement*> &group,
@@ -248,7 +426,7 @@ static void recur_split_ (MElement *e,
}
}
}
-
+//--------------------------------------------------------------
// starting form a list of elements, returns
// lists of lists that are all simply connected
@@ -267,6 +445,7 @@ static void recur_connect (const MEdge &e,
}
}
+//--------------------------------------------------------------
static void connectedRegions (std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > ®ions)
{
@@ -287,7 +466,122 @@ static void connectedRegions (std::vector<MElement*> &elements,
e2e.erase(*it);
}
}
+//--------------------------------------------------------------
+static void recur_cut_ (double R, double a1, double a2,
+ multiscaleLaplaceLevel * root,
+ std::vector<MElement *> &left,
+ std::vector<MElement *> &right ){
+
+ SPoint2 PL (R*cos(a1),R*sin(a1));
+ SPoint2 PR (R*cos(a2),R*sin(a2));
+ std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > centers = root->cut;
+
+ double d = sqrt((PL.x()-PR.x())*(PL.x()-PR.x())+
+ (PL.y()-PR.y())*(PL.y()-PR.y()));
+ SPoint2 farLeft (0.5*(PL.x()+PR.x()) - (PR.y()-PL.y())/d ,
+ 0.5*(PL.y()+PR.y()) + (PR.x()-PL.x())/d );
+
+ for (int i=0;i<root->elements.size();i++){
+ SPoint2 pp (0,0);
+ for (int j=0; j<root->elements[i]->getNumVertices(); j++){
+ pp += root->coordinates[root->elements[i]->getVertex(j)];
+ }
+ pp *= 1./(double)root->elements[i]->getNumVertices();
+ int nbIntersect = 0;
+ for (int j=0;j<centers.size()-1;j++){
+ double x[2];
+ nbIntersect += intersection_segments (centers[j].first,centers[j+1].first,pp,farLeft,x);
+ }
+ if (root->recur != 0){
+ if (nbIntersect %2 == 0)
+ left.push_back(root->elements[i]);
+ else
+ right.push_back(root->elements[i]);
+ }
+ else{
+ if (nbIntersect %2 != 0)
+ left.push_back(root->elements[i]);
+ else
+ right.push_back(root->elements[i]);
+ }
+ }
+
+ for (int i=1;i<centers.size()-1;i++){
+ multiscaleLaplaceLevel* m1 = centers[i-1].second;
+ multiscaleLaplaceLevel* m2 = centers[i].second;
+ multiscaleLaplaceLevel* m3 = centers[i+1].second;
+ if (m2){
+ /*center of the local system is always 0,0
+ its relative position to its parent is center
+ only 2 angles have to be computed for in and out*/
+ a1 = atan2 (m2->center.y() - centers[i-1].first.y(),m2->center.x() - centers[i-1].first.x());
+ a2 = atan2 (m2->center.y() - centers[i+1].first.y(),m2->center.x() - centers[i+1].first.x());
+ recur_cut_ (m2->radius, a1, a2, m2, left, right);
+ }
+ }
+}
+
+//--------------------------------------------------------------
+static void connected_left_right (std::vector<MElement *> &left,
+ std::vector<MElement *> &right ){
+
+ //connected left
+ std::vector<std::vector<MElement*> > subRegionsL;
+ connectedRegions (left,subRegionsL);
+ int indexL=0;
+
+ if (subRegionsL.size() > 0){
+ int maxSize= subRegionsL[0].size();
+ for (int i=1;i< subRegionsL.size() ; i++){
+ int size = subRegionsL[i].size();
+ if(size > maxSize){
+ size = maxSize;
+ indexL = i;
+ }
+ }
+ }
+
+ left.clear();
+ for (int i=0;i< subRegionsL.size() ; i++){
+ if (i == indexL)
+ left.insert(left.begin(), subRegionsL[i].begin(), subRegionsL[i].end());
+ else
+ right.insert(right.begin(), subRegionsL[i].begin(), subRegionsL[i].end());
+ }
+
+ //connected right
+ std::vector<std::vector<MElement*> > subRegionsR;
+ connectedRegions (right,subRegionsR);
+ int indexR=0;
+
+ if (subRegionsR.size() > 0){
+ int maxSize= subRegionsR[0].size();
+ for (int i=1;i< subRegionsR.size() ; i++){
+ int size = subRegionsR[i].size();
+ if(size > maxSize){
+ size = maxSize;
+ indexR = i;
+ }
+ }
+ }
+
+ right.clear();
+ for (int i=0;i< subRegionsR.size() ; i++){
+ if (i == indexR)
+ right.insert(right.begin(), subRegionsR[i].begin(), subRegionsR[i].end());
+ else
+ left.insert(left.begin(), subRegionsR[i].begin(), subRegionsR[i].end());
+ }
+
+ //assign partitions
+ for (int i= 0; i< left.size(); i++)
+ left[i]->setPartition(1);
+ for (int i= 0; i< right.size(); i++)
+ right[i]->setPartition(2);
+
+}
+//--------------------------------------------------------------
static void printLevel ( const char* fn,
std::vector<MElement *> &elements,
std::map<MVertex*,SPoint2> *coordinates,
@@ -321,18 +615,117 @@ static void printLevel ( const char* fn,
fprintf(fp,"$Elements\n%d\n",elements.size());
for (int i=0;i<elements.size();i++){
- elements[i]->writeMSH(fp);
+ elements[i]->writeMSH(fp,version);
}
fprintf(fp,"$EndElements\n%d\n",elements.size());
fclose(fp);
}
+//--------------------------------------------------------------
+static double localSize(MElement *e, std::map<MVertex*,SPoint2> &solution){
+ SBoundingBox3d local;
+ for(unsigned int j = 0; j<e->getNumVertices(); ++j){
+ SPoint2 p = solution[e->getVertex(j)];
+ local += SPoint3(p.x(),p.y(),0.0);
+ }
+ return local.max().distance(local.min());
+}
+//-------------------------------------------------------------
+static void one2OneMap(std::vector<MElement *> &elements, std::map<MVertex*,SPoint2> &solution) {
+
+
+// v2t_cont adjv;
+// std::vector<MTriangle*> allTri;
+// for(int i=0; i< elements.size(); i++){
+// allTri.push_back( (MTriangle*) elements[i] );
+// }
+// buildVertexToTriangle(allTri, adjv);
+
+// for(v2t_cont::iterator it = adjv.begin(); it!= adjv.end(); ++it){
+// MVertex *v = it->first;
+// std::vector<MElement*> vTri = it->second;
+// std::map<MVertex*,SPoint2> vCoord;
+// for (int j=0; j < vTri.size(); j++){
+// for (int k= 0; k < vTri[j]->getNumVertices(); k++){
+// MVertex *vk = vTri[j]->getVertex(k);
+// vCoord[vk] = solution(vk);
+// }
+// }
+// bool badCavity = closedCavity(v,vTri) ? checkCavity(vTri, vCoord) : false;
+
+// if(badCavity){
+// Msg::Debug("Wrong cavity around vertex %d (onwhat=%d).",
+// v->getNum(), v->onWhat()->dim());
+// Msg::Debug("--> Place vertex at center of gravity of %d-Polygon kernel." ,
+// vTri.size());
+
+// double u_cg, v_cg;
+// std::vector<MVertex*> cavV;
+// myPolygon(vTri, cavV);
+// computeCGKernelPolygon(coordinates, cavV, u_cg, v_cg);
+// SPoint3 p_cg(u_cg,v_cg,0);
+// coordinates[v] = p_cg;
+
+// }
+// }
+
+ return;
+
+}
+//--------------------------------------------------------------
+static bool checkOrientation(std::vector<MElement *> &elements,
+ std::map<MVertex*,SPoint2> &solution, int iter) {
+
+ double a_old = 0, a_new;
+ bool oriented = true;
+ int count = 0;
+ for(unsigned int i = 0; i < elements.size(); ++i){
+ MElement *t = elements[i];
+ SPoint2 v1 = solution[t->getVertex(0)];
+ SPoint2 v2 = solution[t->getVertex(1)];
+ SPoint2 v3 = solution[t->getVertex(2)];
+ double p0[2] = {v1[0],v1[1]};
+ double p1[2] = {v2[0],v2[1]};
+ double p2[2] = {v3[0],v3[1]};
+ a_new = robustPredicates::orient2d(p0, p1, p2);
+ if(count == 0) a_old=a_new;
+ if(a_new*a_old < 0.){
+ oriented = false;
+ break;
+ }
+ else{
+ a_old = a_new;
+ }
+ count++;
+ }
+
+ int iterMax = 1;
+ if(!oriented && iter < iterMax){
+ if (iter == 0) Msg::Warning("*** Parametrization is NOT 1 to 1 : applying cavity checks.");
+ Msg::Warning("*** Cavity Check - iter %d -",iter);
+ one2OneMap(elements, solution);
+ return checkOrientation(elements, solution, iter+1);
+ }
+ else if (iter < iterMax){
+ Msg::Debug("Parametrization is 1 to 1 :-)");
+ }
-multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
- std::vector<MVertex*> &boundaryNodes,
- std::vector<double> &linearAbscissa) {
+ return oriented;
+
+}
+//--------------------------------------------------------------
+
+multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements, int iPart)
+{
+
+ //To go through this execute gmsh with the option -optimize_hom
+ //if (!CTX::instance()->mesh.smoothInternalEdges)return;
+
+ //Find the boundary loops
+ //Loop with the largest equivalent radius is the Dirichlet boundary
+ std::vector<std::pair<MVertex*,double> > boundaryNodes;
+ ordering_dirichlet(elements,boundaryNodes);
- if (!CTX::instance()->mesh.smoothInternalEdges)return;
#if defined(HAVE_TAUCS)
_lsys = new linearSystemCSRTaucs<double>;
#elif defined(HAVE_GMM)
@@ -342,93 +735,82 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
#else
_lsys = new linearSystemFull<double>;
#endif
-
+
+ //Assign Dirichlet BCs
root = new multiscaleLaplaceLevel;
root->elements = elements;
for(unsigned int i = 0; i < boundaryNodes.size(); i++){
- MVertex *v = boundaryNodes[i];
- const double theta = 2 * M_PI * linearAbscissa[i];
+ MVertex *v = boundaryNodes[i].first;
+ const double theta = 2 * M_PI * boundaryNodes[i].second;
root->coordinates[v] = SPoint2(cos(theta),sin(theta));
}
+
+ //Recursively parametrize
root->recur = 0;
root->region = 0;
parametrize(*root);
- std::multimap<MEdge,MElement*,Less_Edge> e2e;
- for (int i=0;i<elements.size();++i){
- for (int j=0;j<elements[i]->getNumEdges();j++){
- e2e.insert(std::make_pair(elements[i]->getEdge(j),elements[i]));
- }
- }
- std::map<MEdge,MVertex*,Less_Edge> cutEdges;
- std::set<MVertex*> cutVertices;
- recur_compute_centers_ (1.0,0, M_PI, root);
- recur_cut_edges_ (root,e2e,cutEdges,cutVertices);
-
- if (1){
- // DEBUG : EXPORT -----------------
- std::map<MEdge,MVertex*,Less_Edge>::iterator it = cutEdges.begin();
- FILE *f = fopen("points.pos","w");
- fprintf(f,"View\"\"{\n");
- for ( ; it != cutEdges.end();++it){
- fprintf(f,"SP(%g,%g,%g){1.0};\n",it->second->x(),it->second->y(),it->second->z());
- }
- fprintf(f,"};\n");
- fclose(f);
- // ------END DEBUG
- }
- std::set<MEdge,Less_Edge> theCut;
- std::vector<MElement*> _all;
- recur_cut_elements_ (root,cutEdges,cutVertices,theCut,_all);
- if (1){
- // DEBUG : EXPORT -----------------
- std::set<MEdge,Less_Edge>::iterator it = theCut.begin();
- FILE *f = fopen("edges.pos","w");
- fprintf(f,"View\"\"{\n");
- for ( ; it != theCut.end();++it){
- fprintf(f,"SL(%g,%g,%g,%g,%g,%g){1.0,1.0};\n",it->getVertex(0)->x(),it->getVertex(0)->y(),it->getVertex(0)->z(),
- it->getVertex(1)->x(),it->getVertex(1)->y(),it->getVertex(1)->z());
- }
- fprintf(f,"};\n");
- fclose(f);
- // ------END DEBUG
- }
- e2e.clear();
- for (int i=0;i<_all.size();++i){
- for (int j=0;j<_all[i]->getNumEdges();j++){
- e2e.insert(std::make_pair(_all[i]->getEdge(j),_all[i]));
- }
- }
- // std::set<MElement*> leftSet;
- // recur_split_ (_all[0],e2e,leftSet,theCut);
+ //Compute centers for the cut
+ recur_compute_centers_ (1.0, 0., M_PI, root);
+
+ //Partition the mesh in left and right
+ cut (elements, iPart);
+
+ //---- Testing other cut for partitionning ----
+ //---- cutEdges and connected_regions ----
+// std::multimap<MEdge,MElement*,Less_Edge> e2e;
+// for (int i=0;i<elements.size();++i){
+// for (int j=0;j<elements[i]->getNumEdges();j++){
+// e2e.insert(std::make_pair(elements[i]->getEdge(j),elements[i]));
+// }
+// }
+// std::map<MEdge,MVertex*,Less_Edge> cutEdges;
+// std::set<MVertex*> cutVertices;
+// recur_cut_edges_ (root,e2e,cutEdges,cutVertices);
+
+// printf("Writing points.pos \n");
+// std::map<MEdge,MVertex*,Less_Edge>::iterator ite = cutEdges.begin();
+// FILE *f1 = fopen("points.pos","w");
+// fprintf(f1,"View\"\"{\n");
+// for ( ; ite != cutEdges.end();++ite){
+// fprintf(f1,"SP(%g,%g,%g){1.0};\n",ite->second->x(),ite->second->y(),ite->second->z());
+// }
+// fprintf(f1,"};\n");
+// fclose(f1);
+
+// std::set<MEdge,Less_Edge> theCut;
+// std::vector<MElement*> _all;
+// recur_cut_elements_ (root,cutEdges,cutVertices,theCut,_all);
+
+// printf("Writing edges.pos \n");
+// std::set<MEdge,Less_Edge>::iterator itc = theCut.begin();
+// FILE *f2 = fopen("edges.pos","w");
+// fprintf(f2,"View\"\"{\n");
+// for ( ; itc != theCut.end();++itc){
+// fprintf(f2,"SL(%g,%g,%g,%g,%g,%g){1.0,1.0};\n",itc->getVertex(0)->x(),itc->getVertex(0)->y(),itc->getVertex(0)->z(),
+// itc->getVertex(1)->x(),itc->getVertex(1)->y(),itc->getVertex(1)->z());
+// }
+// fprintf(f2,"};\n");
+// fclose(f2);
+
+// e2e.clear();
+// for (int i=0;i<_all.size();++i){
+// for (int j=0;j<_all[i]->getNumEdges();j++){
+// e2e.insert(std::make_pair(_all[i]->getEdge(j),_all[i]));
+// }
+// }
+// std::set<MElement*> leftSet;
+// recur_split_ (_all[0],e2e,leftSet,theCut);
+// for (int i=0;i<_all.size();i++){
+// if(leftSet.find(_all[i]) == leftSet.end())right.push_back(_all[i]);
+// else left.push_back(_all[i]);
+// }
- std::vector<MElement*> left,right;
- /*
- for (int i=0;i<_all.size();i++){
- if(leftSet.find(_all[i]) == leftSet.end())right.push_back(_all[i]);
- else left.push_back(_all[i]);
- }
- */
- cut (left,right);
-
- printLevel ("multiscale_left.msh",left,0,1.0);
- printLevel ("multiscale_right.msh",right,0,1.0);
- printLevel ("multiscale_all.msh",_all,0,1.0);
- // FIXME !!!!!
- throw;
-}
-
-static double localSize(MElement *e, std::map<MVertex*,SPoint2> &solution){
- SBoundingBox3d local;
- for(unsigned int j = 0; j<e->getNumVertices(); ++j){
- SPoint2 p = solution[e->getVertex(j)];
- local += SPoint3(p.x(),p.y(),0.0);
- }
- return local.max().distance(local.min());
}
void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
- simpleFunction<double> ONE(1.0);
+
+
std::set<MVertex*> allNodes;
for(unsigned int i = 0; i < level.elements.size(); ++i){
MElement *e = level.elements[i];
@@ -436,41 +818,16 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
allNodes.insert(e->getVertex(j));
}
}
+
+ //Parametrize level
std::map<MVertex*,SPoint2> solution;
- // PARAMETRIZE ---------------------------------
- for (int step =0 ; step<2 ; step++){
- dofManager<double> myAssembler(_lsys);
- for(std::map<MVertex*,SPoint2>::iterator it = level.coordinates.begin();
- it != level.coordinates.end(); ++it){
- MVertex *v = it->first;
- myAssembler.fixVertex(v, 0, 1, it->second[step]);
- }
- // do the numbering
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
- MVertex *v = *itv;
- myAssembler.numberVertex(v, 0, 1);
- }
- // assemble
- laplaceTerm laplace(0, 1, &ONE);
- for(unsigned int i = 0; i < level.elements.size(); ++i){
- MElement *e = level.elements[i];
- SElement se(e);
- laplace.addToMatrix(myAssembler, &se);
- }
- // solve
- if (myAssembler.sizeOfR() != 0) _lsys->systemSolve();
- // get the values
- int count = 0;
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
- MVertex *v = *itv;
- double value = myAssembler.getDofValue(v, 0, 1);
- if (step == 0)solution[v] = SPoint2(value,0);
- else solution[v] = SPoint2(solution[v][0],value);
- }
- _lsys->clear();
+ parametrize_method(level, allNodes, solution, HARMONIC);
+ if (!checkOrientation(level.elements, solution, 0)){
+ Msg::Info("Parametrization failed using standard techniques : moving to convex combination");
+ parametrize_method(level, allNodes, solution, CONVEXCOMBINATION);
}
-
- // compute the bbox of the parametric space
+
+ //Compute the bbox of the parametric space
SBoundingBox3d bbox;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
@@ -478,7 +835,8 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
bbox += SPoint3(p.x(),p.y(),0.0);
}
double global_size = bbox.max().distance(bbox.min());
- // check elements that are too small
+
+ //Check elements that are too small
std::vector<MElement*> tooSmall, goodSize;
for(unsigned int i = 0; i < level.elements.size(); ++i){
@@ -492,14 +850,14 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
goodSize.push_back(e);
}
}
- // only keep the right elements and nodes
+
+ //Only keep the right elements and nodes
level.elements = goodSize;
std::vector<std::vector<MElement*> > regions, regions_;
connectedRegions (tooSmall,regions_);
- // remove small regions
-
+ //Remove small regions
for (int i=0;i< regions_.size() ; i++){
bool region_has_really_small_elements = false;
for (int k=0; k<regions_[i].size() ; k++){
@@ -522,16 +880,16 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
}
}
- //END PARAMETRIZE ---------------------------------
// DEBUG
char name[245];
sprintf(name,"multiscale_level%d_region%d_real.msh",level.recur, level.region);
- printLevel (name,level.elements,0,1.0);
+ printLevel (name,level.elements,0,2.0);
sprintf(name,"multiscale_level%d_region%d_param.msh",level.recur, level.region);
- printLevel (name,level.elements,&level.coordinates,1.0);
+ printLevel (name,level.elements,&level.coordinates,2.0);
// END DEBUG
- Msg::Info("%d connected regions\n",regions.size());
+ if (regions.size() >0)
+ Msg::Info("%d connected regions\n",regions.size());
for (int i=0;i< regions.size() ; i++){
// check nodes that appear both on too small and good size
@@ -567,7 +925,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
// recursively continue if tooSmall is not empty
if (!tooSmall.empty()){
Msg::Info("Multiscale Laplace, level %d region %d, %d too small",level.recur,level.region,tooSmall.size());
- level.childeren.push_back(nextLevel);
+ level.children.push_back(nextLevel);
parametrize (*nextLevel);
}
else {
@@ -577,85 +935,75 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
}
}
-static void recur_cut_ (double R, double a1, double a2,
- multiscaleLaplaceLevel * root,
- std::vector<MElement *> &left,
- std::vector<MElement *> &right ){
- std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers = root->cut;
- centers.clear();
- multiscaleLaplaceLevel* zero = 0;
+void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
+ std::set<MVertex*> &allNodes,
+ std::map<MVertex*,SPoint2> &solution,
+ typeOfMapping tom)
+{
- SPoint2 PL (R*cos(a1),R*sin(a1));
- SPoint2 PR (R*cos(a2),R*sin(a2));
- centers.push_back(std::make_pair(PL,zero));
- centers.push_back(std::make_pair(PR,zero));
- for (int i=0;i<root->childeren.size();i++){
- centers.push_back(std::make_pair(root->childeren[i]->center,root->childeren[i]));
- multiscaleLaplaceLevel* m = root->childeren[i];
- m->radius = 0.0;
- for (std::map<MVertex*,SPoint2>::iterator it = m->coordinates.begin();
- it != m->coordinates.end() ; ++it){
- const SPoint2 &p = it->second;
- m->radius = std::max(m->radius,sqrt ((m->center.x() - p.x())*(m->center.x() - p.x())+
- (m->center.y() - p.y())*(m->center.y() - p.y())));
- }
- }
- std::sort(centers.begin(),centers.end());
- /*
- printf("%d %d ",root->recur,root->region);
- for (int j=0;j<centers.size();j++){
- printf("(%g %g) ",centers[j].first.x(),centers[j].first.y());
- }
- printf("\n");
- */
- double d = sqrt((PL.x()-PR.x())*(PL.x()-PR.x())+
- (PL.y()-PR.y())*(PL.y()-PR.y()));
- SPoint2 farLeft (0.5*(PL.x()+PR.x()) - (PR.y()-PL.y())/d ,
- 0.5*(PL.y()+PR.y()) + (PR.x()-PL.x())/d );
-
- for (int i=0;i<root->elements.size();i++){
- SPoint2 pp (0,0);
- for (int j=0; j<root->elements[i]->getNumVertices(); j++){
- pp += root->coordinates[root->elements[i]->getVertex(j)];
+ solution.clear();
+ simpleFunction<double> ONE(1.0);
+
+ for (int step =0 ; step<2 ; step++){
+
+ dofManager<double> myAssembler(_lsys);
+ for(std::map<MVertex*,SPoint2>::iterator it = level.coordinates.begin();
+ it != level.coordinates.end(); ++it){
+ MVertex *v = it->first;
+ myAssembler.fixVertex(v, 0, 1, it->second[step]);
}
- pp *= 1./(double)root->elements[i]->getNumVertices();
- int nbIntersect = 0;
- for (int j=0;j<centers.size()-1;j++){
- double x[2];
- nbIntersect += intersection_segments (centers[j].first,centers[j+1].first,pp,farLeft,x);
+
+ // do the numbering
+ for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ MVertex *v = *itv;
+ myAssembler.numberVertex(v, 0, 1);
}
- if (root->recur != 0){
- if (nbIntersect %2 == 0)
- left.push_back(root->elements[i]);
- else
- right.push_back(root->elements[i]);
+
+ // assemble
+ femTerm<double> *mapping;
+ if (tom == HARMONIC){
+ mapping = new laplaceTerm(0, 1, &ONE);
}
- else{
- if (nbIntersect %2 != 0)
- left.push_back(root->elements[i]);
- else
- right.push_back(root->elements[i]);
+ else if (tom == CONVEXCOMBINATION){
+ mapping = new convexCombinationTerm(0, 1, &ONE);
}
- }
-
- for (int i=1;i<centers.size()-1;i++){
- multiscaleLaplaceLevel* m1 = centers[i-1].second;
- multiscaleLaplaceLevel* m2 = centers[i].second;
- multiscaleLaplaceLevel* m3 = centers[i+1].second;
- if (m2){
- /*center of the local system is always 0,0
- its relative position to its parent is center
- only 2 angles have to be computed for in and out*/
- a1 = atan2 (m2->center.y() - centers[i-1].first.y(),m2->center.x() - centers[i-1].first.x());
- a2 = atan2 (m2->center.y() - centers[i+1].first.y(),m2->center.x() - centers[i+1].first.x());
- recur_cut_ (m2->radius, a1, a2, m2, left, right);
+
+ for(unsigned int i = 0; i < level.elements.size(); ++i){
+ MElement *e = level.elements[i];
+ SElement se(e);
+ mapping->addToMatrix(myAssembler, &se);
}
+
+ // solve
+ if (myAssembler.sizeOfR() != 0) _lsys->systemSolve();
+
+ // get the values
+ int count = 0;
+ for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ MVertex *v = *itv;
+ double value = myAssembler.getDofValue(v, 0, 1);
+ if (step == 0)solution[v] = SPoint2(value,0);
+ else solution[v] = SPoint2(solution[v][0],value);
+ }
+ _lsys->clear();
+
}
}
+void multiscaleLaplace::cut (std::vector<MElement *> &elements, int iPart)
+{
+ std::vector<MElement*> left,right;
+ recur_cut_ (1.0, 0, M_PI, root,left,right);
+
+ connected_left_right(left, right);
+
+ elements.clear();
+ elements.insert(elements.end(),left.begin(),left.end());
+ elements.insert(elements.end(),right.begin(),right.end());
+
+ char name[256];
+ sprintf(name, "laplace_%d.msh", iPart);
+ printLevel (name,elements,0,2.0);
-void multiscaleLaplace::cut (std::vector<MElement *> &left,
- std::vector<MElement *> &right)
-{
- recur_cut_ (1.0,0, M_PI, root,left,right);
}
+
diff --git a/Solver/multiscaleLaplace.h b/Solver/multiscaleLaplace.h
index 90f662d261c658a7b957a3e244f39c8a5080165b..8b456290e2cef5c5732d68f38c5d256e574092cf 100644
--- a/Solver/multiscaleLaplace.h
+++ b/Solver/multiscaleLaplace.h
@@ -3,6 +3,7 @@
#include <vector>
#include <map>
+#include <set>
#include "SPoint2.h"
#include "linearSystem.h"
@@ -14,21 +15,26 @@ struct multiscaleLaplaceLevel {
double scale;
double radius;
int recur,region;
- std::vector<multiscaleLaplaceLevel*> childeren;
+ std::vector<multiscaleLaplaceLevel*> children;
std::vector<MElement *> elements;
std::map<MVertex*,SPoint2> coordinates;
std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > cut;
};
class multiscaleLaplace{
+public:
+ multiscaleLaplace (std::vector<MElement *> &elements, int iPart=0);
+ void cut (std::vector<MElement *> &elements,int iPart=0);
+ typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3} typeOfMapping;
+
linearSystem<double> *_lsys;
multiscaleLaplaceLevel* root;
- void parametrize (multiscaleLaplaceLevel &);
-public:
- multiscaleLaplace (std::vector<MElement *> &elements,
- std::vector<MVertex*> &boundaryNodes,
- std::vector<double> &linearAbscissa) ;
- void cut (std::vector<MElement *> &left,
- std::vector<MElement *> &right);
+ void parametrize (multiscaleLaplaceLevel &);
+ void parametrize_method (multiscaleLaplaceLevel & level,
+ std::set<MVertex*> &allNodes,
+ std::map<MVertex*,SPoint2> &solution,
+ typeOfMapping tom);
+
+
};
#endif
diff --git a/benchmarks/boolean/JAW.geo b/benchmarks/boolean/JAW.geo
index 0a222a6161aa31677b77c87f00aa9dba880705d0..2430dcbcd065309cf06c2c2f5980044964d15083 100644
--- a/benchmarks/boolean/JAW.geo
+++ b/benchmarks/boolean/JAW.geo
@@ -1,3 +1,6 @@
+Mesh.CharacteristicLengthFactor=0.15;
+
+
L = 100;
H = 30;
Z = 10;
@@ -12,3 +15,5 @@ EndFor
OCCShape("Sphere",{H-X,H/2,Z/2,R},"Fuse");
OCCShape("Torus",{L,H/2,Z/2,0,0,1,2*R,R/2},"Fuse");
+
+Compound Surface(100) = {1 ... 26};
\ No newline at end of file
diff --git a/benchmarks/extrude/sphere_boundary_layer.geo b/benchmarks/extrude/sphere_boundary_layer.geo
index 66014e14b307ddaeb5d31f7bce7a2f8ea9982c3e..f181894dd9c7d43689a3885ccd910254fa16b88d 100644
--- a/benchmarks/extrude/sphere_boundary_layer.geo
+++ b/benchmarks/extrude/sphere_boundary_layer.geo
@@ -38,7 +38,7 @@ Line Loop(27) = {-4,12,-6};
Ruled Surface(28) = {27};
Extrude {
- Surface{14:28:2}; Layers{10, 0.1}; // Recombine;
+ Surface{14:28:2}; //Layers{10, 0.1}; // Recombine;
}
Surface Loop(29) = {28,26,16,14,20,24,22,18};