From bd3997c7534ebd48a739ab7eaca16cbfefa1e38d Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Thu, 11 Mar 2010 15:09:21 +0000
Subject: [PATCH]
---
Geo/CMakeLists.txt | 1 +
Geo/GFaceCompound.cpp | 2 +-
Geo/GModel.cpp | 59 +------
Geo/discreteFace.cpp | 5 +-
Mesh/BackgroundMesh.cpp | 289 +++++++++++++++++++++++++++++++++++
Mesh/BackgroundMesh.h | 32 ++++
Mesh/highOrderSmoother.cpp | 29 ++--
Mesh/meshGFace.cpp | 16 +-
Mesh/meshGFaceBDS.cpp | 61 ++------
Mesh/meshGFaceLloyd.cpp | 14 +-
Numeric/DivideAndConquer.cpp | 13 +-
Numeric/DivideAndConquer.h | 8 +-
Solver/dgAlgorithm.cpp | 6 -
Solver/dgAlgorithm.h | 6 +
Solver/dgGroupOfElements.cpp | 9 ++
Solver/dgGroupOfElements.h | 3 +
Solver/helmholtzTerm.h | 2 +-
Solver/linearSystemCSR.cpp | 3 +-
Solver/multiscaleLaplace.cpp | 63 +++++---
Solver/multiscaleLaplace.h | 1 +
benchmarks/boolean/JAW.geo | 4 +-
benchmarks/boolean/TORUS.geo | 4 +-
benchmarks/step/capot.geo | 4 +-
23 files changed, 473 insertions(+), 161 deletions(-)
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 300a633e0a..f8b178b70e 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -28,6 +28,7 @@ set(SRC
MVertex.cpp
MFace.cpp
MElement.cpp
+ MElementOctree.cpp
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp
MHexahedron.cpp MPrism.cpp MPyramid.cpp MElementCut.cpp
MZone.cpp MZoneBoundary.cpp
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index cc843299dd..762018d202 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -790,7 +790,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
_lsysb->setGmres(1);
_lsys = _lsysb;
-#elif defined(_HAVE_TAUCS)
+#elif defined(HAVE_TAUCS)
_lsys = new linearSystemCSRTaucs<double>;
#else
_lsys = new linearSystemFull<double>;
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index b995197777..4f8059902d 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -17,6 +17,7 @@
#include "MPrism.h"
#include "MPyramid.h"
#include "MElementCut.h"
+#include "MElementOctree.h"
#include "discreteRegion.h"
#include "discreteFace.h"
#include "discreteEdge.h"
@@ -574,63 +575,12 @@ int GModel::getNumMeshElements(unsigned c[5])
return 0;
}
-static void MElementBB(void *a, double *min, double *max)
-{
- MElement *e = (MElement*)a;
- MVertex *v = e->getVertex(0);
- min[0] = max[0] = v->x();
- min[1] = max[1] = v->y();
- min[2] = max[2] = v->z();
- for(int i = 1; i < e->getNumVertices(); i++){
- v = e->getVertex(i);
- min[0] = std::min(min[0], v->x()); max[0] = std::max(max[0], v->x());
- min[1] = std::min(min[1], v->y()); max[1] = std::max(max[1], v->y());
- min[2] = std::min(min[2], v->z()); max[2] = std::max(max[2], v->z());
- }
-}
-
-static void MElementCentroid(void *a, double *x)
-{
- MElement *e = (MElement*)a;
- MVertex *v = e->getVertex(0);
- int n = e->getNumVertices();
- x[0] = v->x(); x[1] = v->y(); x[2] = v->z();
- for(int i = 1; i < n; i++) {
- v = e->getVertex(i);
- x[0] += v->x(); x[1] += v->y(); x[2] += v->z();
- }
- double oc = 1. / (double)n;
- x[0] *= oc;
- x[1] *= oc;
- x[2] *= oc;
-}
-
-static int MElementInEle(void *a, double *x)
-{
- MElement *e = (MElement*)a;
- double uvw[3];
- e->xyz2uvw(x, uvw);
- return e->isInside(uvw[0], uvw[1], uvw[2]) ? 1 : 0;
-}
MElement *GModel::getMeshElementByCoord(SPoint3 &p)
{
if(!_octree){
Msg::Debug("Rebuilding mesh element octree");
- SBoundingBox3d bb = bounds();
- double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()};
- double size[3] = {bb.max().x() - bb.min().x(),
- bb.max().y() - bb.min().y(),
- bb.max().z() - bb.min().z()};
- const int maxElePerBucket = 100; // memory vs. speed trade-off
- _octree = Octree_Create(maxElePerBucket, min, size,
- MElementBB, MElementCentroid, MElementInEle);
- std::vector<GEntity*> entities;
- getEntities(entities);
- for(unsigned int i = 0; i < entities.size(); i++)
- for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
- Octree_Insert(entities[i]->getMeshElement(j), _octree);
- Octree_Arrange(_octree);
+ _octree = buildMElementOctree(this);
}
double P[3] = {p.x(), p.y(), p.z()};
return (MElement*)Octree_Search(P, _octree);
@@ -1183,12 +1133,17 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
discEdges.push_back((discreteEdge*) *it);
}
+ printf("coucou1\n");
+
// find boundary edges of each face and put them in a map_edges that
// associates the MEdges with the tags of the adjacent faces
std::map<MEdge, std::vector<int>, Less_Edge > map_edges;
for (std::vector<discreteFace*>::iterator it = discFaces.begin();
it != discFaces.end(); it++)
(*it)->findEdges(map_edges);
+
+ printf("coucou2\n");
+
//return if no boundary edges (torus, sphere, ...)
if (map_edges.empty()) return;
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 2f53c3098b..01c6987815 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -24,6 +24,7 @@ discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_edges)
{
+ printf("coucou 3\n");
std::set<MEdge, Less_Edge> bound_edges;
for (unsigned int iFace = 0; iFace < getNumMeshElements() ; iFace++) {
MElement *e = getMeshElement(iFace);
@@ -34,6 +35,7 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_e
else bound_edges.erase(itset);
}
}
+ printf("coucou 4\n");
// for the boundary edges, associate the tag of the current discrete face
for (std::set<MEdge, Less_Edge>::iterator itv = bound_edges.begin();
@@ -49,7 +51,8 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge> &map_e
tagFaces.push_back(tag());
itmap->second = tagFaces;
}
- }
+ }
+ printf("coucou 5\n");
}
void discreteFace::setBoundEdges(std::vector<int> tagEdges)
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 657ca15d54..1abb2c276a 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -12,6 +12,18 @@
#include "GFace.h"
#include "GModel.h"
#include "Field.h"
+#include "MElementOctree.h"
+#include "MElement.h"
+#include "MLine.h"
+#include "MTriangle.h"
+#include "MVertex.h"
+#include "Octree.h"
+#include "dofManager.h"
+#include "laplaceTerm.h"
+#include "linearSystemGMM.h"
+#include "linearSystemCSR.h"
+#include "linearSystemFull.h"
+#include "linearSystemPETSc.h"
// computes the characteristic length of the mesh at a vertex in order
// to have the geometry captured with accuracy. A parameter called
@@ -232,3 +244,280 @@ bool Extend2dMeshIn3dVolumes()
{
return CTX::instance()->mesh.lcExtendFromBoundary ? true : false;
}
+
+// ---------- backgroundMesh class -----------
+
+void backgroundMesh::set (GFace *gf){
+ if (_current) delete _current;
+ _current = new backgroundMesh(gf);
+}
+
+void backgroundMesh::unset (){
+ if (_current) delete _current;
+ _current = 0;
+}
+
+backgroundMesh::backgroundMesh (GFace *_gf)
+{
+ for (int i=0;i<_gf->triangles.size();i++){
+ MTriangle *e = _gf->triangles[i];
+ MVertex *news[3];
+ for (int j=0;j<3;j++){
+ std::map<MVertex*,MVertex*>::iterator it = _3Dto2D.find(e->getVertex(j));
+ MVertex *newv =0;
+ if (it == _3Dto2D.end()){
+ SPoint2 p;
+ bool success = reparamMeshVertexOnFace(e->getVertex(j), _gf, p);
+ newv = new MVertex (p.x(),p.y(),0.0);
+ _vertices.push_back(newv);
+ _3Dto2D[e->getVertex(j)] = newv;
+ _2Dto3D[newv] = e->getVertex(j);
+ }
+ else newv = it->second;
+ news[j] = newv;
+ }
+ _triangles.push_back(new MTriangle(news[0],news[1],news[2]));
+ }
+ _octree = buildMElementOctree (_triangles);
+ if (CTX::instance()->mesh.lcFromPoints)
+ propagate1dMesh(_gf);
+ else {
+ std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
+ for ( ; itv2 != _2Dto3D.end(); ++itv2){
+ _sizes[itv2->first] = 1.e22;
+ }
+ }
+ updateSizes(_gf);
+ _3Dto2D.clear();
+ _2Dto3D.clear();
+}
+
+backgroundMesh::~backgroundMesh (){
+ for (int i=0;i<_vertices.size();i++) delete _vertices[i];
+ for (int i=0;i<_triangles.size();i++) delete _triangles[i];
+ delete _octree;
+}
+
+void backgroundMesh::propagate1dMesh (GFace * _gf)
+{
+ std::list<GEdge*> e = _gf->edges();
+ std::list<GEdge*>::const_iterator it = e.begin();
+ std::map<MVertex*,double> sizes;
+
+ linearSystem<double> *_lsys = 0;
+#if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
+ _lsys = new linearSystemPETSc<double>;
+#elif defined(HAVE_GMM) && !defined(HAVE_TAUCS)
+ linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
+ _lsysb->setGmres(1);
+ _lsys = _lsysb;
+#elif defined(HAVE_TAUCS)
+ _lsys = new linearSystemCSRTaucs<double>;
+#else
+ _lsys = new linearSystemFull<double>;
+#endif
+
+ dofManager<double> myAssembler(_lsys);
+
+ for( ; it != e.end(); ++it ){
+ if (!(*it)->isSeam(_gf)){
+ for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
+ MVertex *v1 = (*it)->lines[i]->getVertex(0);
+ MVertex *v2 = (*it)->lines[i]->getVertex(1);
+ double d = sqrt ((v1->x()-v2->x())*(v1->x()-v2->x())+
+ (v1->y()-v2->y())*(v1->y()-v2->y())+
+ (v1->z()-v2->z())*(v1->z()-v2->z()));
+ for (int k=0;k<2;k++){
+ MVertex *v = (*it)->lines[i]->getVertex(k);
+ std::map<MVertex*,double>::iterator itv = sizes.find(v);
+ if (itv==sizes.end())
+ sizes[v] = d;
+ else
+ itv->second = 0.5* ( itv->second + d );
+ }
+ }
+ }
+ }
+ std::map<MVertex*,double>::iterator itv = sizes.begin();
+ for ( ; itv != sizes.end(); ++itv){
+ myAssembler.fixVertex(itv->first, 0, 1, itv->second);
+ }
+
+ simpleFunction<double> ONE(1.0);
+ laplaceTerm l(0, 1, &ONE);
+
+ for (unsigned int k= 0; k< _gf->triangles.size(); k++){
+ MTriangle *t = _gf->triangles[k];
+ myAssembler.numberVertex(t->getVertex(0), 0, 1);
+ myAssembler.numberVertex(t->getVertex(1), 0, 1);
+ myAssembler.numberVertex(t->getVertex(2), 0, 1);
+ }
+
+ for (unsigned int k= 0; k< _gf->triangles.size(); k++){
+ MTriangle *t = _gf->triangles[k];
+ SElement se(t);
+ l.addToMatrix(myAssembler, &se);
+ }
+ _lsys->systemSolve();
+
+ std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
+ for ( ; itv2 != _2Dto3D.end(); ++itv2){
+ MVertex *v_2D = itv2->first;
+ MVertex *v_3D = itv2->second;
+ double value = myAssembler.getDofValue(v_3D, 0, 1);
+ _sizes[v_2D] = value;
+ }
+ delete _lsys;
+}
+
+void backgroundMesh::updateSizes (GFace * _gf){
+ std::map<MVertex*,double>::iterator itv = _sizes.begin();
+ for ( ; itv != _sizes.end(); ++itv){
+ SPoint2 p;
+ MVertex *v = _2Dto3D[itv->first];
+ double lc;
+ if (v->onWhat()->dim() == 0){
+ lc = BGM_MeshSize(v->onWhat(), 0,0,v->x(),v->y(),v->z());
+ }
+ else if (v->onWhat()->dim() == 1){
+ double u;
+ v->getParameter(0,u);
+ lc = BGM_MeshSize(v->onWhat(), u,0,v->x(),v->y(),v->z());
+ }
+ else{
+ bool success = reparamMeshVertexOnFace(v, _gf, p);
+ lc = BGM_MeshSize(_gf, p.x(),p.y(),v->x(),v->y(),v->z());
+ }
+ // printf("2D -- %g %g 3D -- %g %g\n",p.x(),p.y(),v->x(),v->y());
+ itv->second = std::min(lc,itv->second);
+ itv->second = std::max(itv->second, CTX::instance()->mesh.lcMin);
+ itv->second = std::min(itv->second, CTX::instance()->mesh.lcMax);
+ }
+ // return;
+ // ---------------------
+ // now do some diffusion
+ // ---------------------
+
+
+ std::list<GEdge*> e = _gf->edges();
+ std::list<GEdge*>::const_iterator it = e.begin();
+
+ linearSystem<double> *_lsys = 0;
+#if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
+ _lsys = new linearSystemPETSc<double>;
+#elif defined(HAVE_GMM) && !defined(HAVE_TAUCS)
+ linearSystemGmm<double> *_lsysb = new linearSystemGmm<double>;
+ _lsysb->setGmres(1);
+ _lsys = _lsysb;
+#elif defined(HAVE_TAUCS)
+ _lsys = new linearSystemCSRTaucs<double>;
+#else
+ _lsys = new linearSystemFull<double>;
+#endif
+
+ dofManager<double> myAssembler(_lsys);
+
+ // M (U^{t+1} - U^{t})/DT + K U^{t+1} = 0
+
+ for( ; it != e.end(); ++it ){
+ if (!(*it)->isSeam(_gf)){
+ for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
+ MVertex *v1 = (*it)->lines[i]->getVertex(0);
+ MVertex *v2 = (*it)->lines[i]->getVertex(1);
+ MVertex *v1_2D = _3Dto2D[ v1 ];
+ MVertex *v2_2D = _3Dto2D[ v2 ];
+
+ myAssembler.fixVertex(v1, 0, 1, _sizes[v1_2D]);
+ myAssembler.fixVertex(v2, 0, 1, _sizes[v2_2D]);
+ }
+ }
+ }
+
+ double DT = 0.01;
+ simpleFunction<double> ONEOVERDT(1.0/DT);
+ simpleFunction<double> ONE(1.0);
+ helmholtzTerm<double> diffusionTerm(0,1, 1, &ONE,&ONEOVERDT);
+ helmholtzTerm<double> massTerm(0,1, 1, 0,&ONEOVERDT);
+
+ for (unsigned int k= 0; k< _gf->triangles.size(); k++){
+ MTriangle *t = _gf->triangles[k];
+ myAssembler.numberVertex(t->getVertex(0), 0, 1);
+ myAssembler.numberVertex(t->getVertex(1), 0, 1);
+ myAssembler.numberVertex(t->getVertex(2), 0, 1);
+ }
+ fullMatrix<double> mass_matrix(3,3);
+ for (unsigned int k= 0; k< _gf->getNumMeshElements(); k++){
+ MElement *e = _gf->getMeshElement(k);
+ SElement se(e);
+ diffusionTerm.addToMatrix(myAssembler, &se);
+ massTerm.elementMatrix(&se,mass_matrix);
+ for (int i=0;i<e->getNumVertices();i++){
+ double TERM = 0.0;
+ for (int j=0;j<e->getNumVertices();j++){
+ MVertex *v_2D = _3Dto2D[ e->getVertex(j) ];
+ TERM+=_sizes[v_2D] * mass_matrix(i,j);
+ }
+ myAssembler.assemble(e->getVertex(i),0,1,TERM);
+ }
+ }
+ _lsys->systemSolve();
+
+ std::map<MVertex*,MVertex*>::iterator itv2 = _2Dto3D.begin();
+ for ( ; itv2 != _2Dto3D.end(); ++itv2){
+ MVertex *v_2D = itv2->first;
+ MVertex *v_3D = itv2->second;
+ double value = myAssembler.getDofValue(v_3D, 0, 1);
+ _sizes[v_2D] = value;
+ }
+ delete _lsys;
+}
+
+double backgroundMesh::operator() (double u, double v, double w) const {
+ double uv[3] = {u, v, w};
+ double uv2[3];
+ // return 1.0;
+ MElement *e = (MElement*)Octree_Search(uv, _octree);
+ if (!e) {
+ Msg::Error("cannot find %g %g",u,v);
+ return 1.0;
+ }
+ e->xyz2uvw (uv,uv2);
+ std::map<MVertex*,double>::const_iterator itv1 = _sizes.find(e->getVertex(0));
+ std::map<MVertex*,double>::const_iterator itv2 = _sizes.find(e->getVertex(1));
+ std::map<MVertex*,double>::const_iterator itv3 = _sizes.find(e->getVertex(2));
+ // printf("%g %g -> %g\n",u,v, itv1->second * (1-uv2[0]-uv2[2]) + itv2->second * uv2[0] + itv3->second * uv2[1]);
+ return itv1->second * (1-uv2[0]-uv2[1]) + itv2->second * uv2[0] + itv3->second * uv2[1];
+}
+
+void backgroundMesh::print (const std::string &filename, GFace *gf) const {
+ FILE *f = fopen (filename.c_str(),"w");
+ fprintf(f,"View \"Background Mesh\"{\n");
+ for(int i=0;i<_triangles.size();i++){
+ MVertex *v1 = _triangles[i]->getVertex(0);
+ MVertex *v2 = _triangles[i]->getVertex(1);
+ MVertex *v3 = _triangles[i]->getVertex(2);
+ std::map<MVertex*,double>::const_iterator itv1 = _sizes.find(v1);
+ std::map<MVertex*,double>::const_iterator itv2 = _sizes.find(v2);
+ std::map<MVertex*,double>::const_iterator itv3 = _sizes.find(v3);
+ if (!gf){
+ fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
+ v1->x(),v1->y(),v1->z(),
+ v2->x(),v2->y(),v2->z(),
+ v3->x(),v3->y(),v3->z(),itv1->second,itv2->second,itv3->second);
+ }
+ else {
+ GPoint p1 = gf->point(SPoint2(v1->x(),v1->y()));
+ GPoint p2 = gf->point(SPoint2(v2->x(),v2->y()));
+ GPoint p3 = gf->point(SPoint2(v3->x(),v3->y()));
+ fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%g,%g,%g};\n",
+ p1.x(),p1.y(),p1.z(),
+ p2.x(),p2.y(),p2.z(),
+ p3.x(),p3.y(),p3.z(),itv1->second,itv2->second,itv3->second);
+ }
+
+ }
+ fprintf(f,"};\n");
+ fclose(f);
+
+}
+backgroundMesh* backgroundMesh::_current = 0;
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index 28194642cf..1cb0b2aaed 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -7,12 +7,44 @@
#define _BACKGROUND_MESH_H_
#include "STensor3.h"
+#include "simpleFunction.h"
+#include <vector>
+class Octree;
+class GFace;
+class MElement;
+class MVertex;
class GEntity;
+// this background mesh makes use
+
+class backgroundMesh : public simpleFunction<double>
+{
+ Octree *_octree;
+ std::vector<MVertex*> _vertices;
+ std::vector<MElement*> _triangles;
+ std::map<MVertex*,double> _sizes;
+ std::map<MVertex*,MVertex*> _3Dto2D;
+ std::map<MVertex*,MVertex*> _2Dto3D;
+ static backgroundMesh * _current;
+ backgroundMesh(GFace *);
+ ~backgroundMesh();
+public:
+ static void set (GFace *);
+ static void unset ();
+ static backgroundMesh * current () {return _current;}
+ void propagate1dMesh(GFace *);
+ void updateSizes(GFace *);
+ double operator () (double u, double v, double w) const;
+ void print (const std::string &filename, GFace *gf) const;
+};
+
+
double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double Z);
SMetric3 BGM_MeshMetric(GEntity *ge, double U, double V, double X, double Y, double Z);
bool Extend1dMeshIn2dSurfaces();
bool Extend2dMeshIn3dVolumes();
+
+
#endif
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index 5bc0c843b2..ee7a2f9dba 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -104,7 +104,7 @@ struct p2data{
highOrderSmoother *_s)
: gf(_gf), t1(_t1), t2(_t2), n12(_n12), s(_s)
{
- elasticityTerm el(0, 1.e3, .3333,1);
+ elasticityTerm el(0, 1.e3, .48,1);
s->moveToStraightSidedLocation(t1);
s->moveToStraightSidedLocation(t2);
m1 = new fullMatrix<double>(3 * t1->getNumVertices(),
@@ -368,8 +368,8 @@ void highOrderSmoother::optimize(GFace * gf,
// then try to swap for better configurations
smooth(gf, true);
- //int nbSwap =
- //swapHighOrderTriangles(gf,edgeVertices,faceVertices,this);
+ int nbSwap =
+ swapHighOrderTriangles(gf,edgeVertices,faceVertices,this);
// smooth(gf,true);
// smooth(gf,true);
// smooth(gf,true);
@@ -437,7 +437,7 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*> & all, GFace *gf)
lsys->setPrec(5.e-8);
#endif
dofManager<double> myAssembler(lsys);
- elasticityTerm El(0, 1.0, .333, getTag());
+ elasticityTerm El(0, 1.0, .48, getTag());
std::vector<MElement*> layer, v;
@@ -537,10 +537,10 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*> & all, GFace *gf)
}
double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
- GFace *gf,
- dofManager<double> &myAssembler,
- std::set<MVertex*> &verticesToMove,
- elasticityTerm &El)
+ GFace *gf,
+ dofManager<double> &myAssembler,
+ std::set<MVertex*> &verticesToMove,
+ elasticityTerm &El)
{
std::set<MVertex*>::iterator it;
@@ -711,8 +711,8 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
}
void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
- std::vector<MElement*>& old_elems,
- GFace* gf) {
+ std::vector<MElement*>& old_elems,
+ GFace* gf) {
printf("Smoothing a cavity...\n");
printf("Old elems : %d and %d\n", old_elems[0]->getNum(), old_elems[1]->getNum());
@@ -741,7 +741,6 @@ void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
v.insert(v.end(), gf->quadrangles.begin(),gf->quadrangles.end());
addOneLayer(v,cavity,layer);
- //addOneLayer(v,old_elems,layer);
for (unsigned int i = 0; i < layer.size(); i++){
if (find(old_elems.begin(), old_elems.end(), layer[i]) == old_elems.end()) {
@@ -750,7 +749,7 @@ void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
MVertex *vert = layer[i]->getVertex(j);
myAssembler.fixVertex(vert, 0, getTag(), 0);
myAssembler.fixVertex(vert, 1, getTag(), 0);
- printf("Fixing vertex %d\n", vert->getNum());
+ printf("Fixing vertex (internal) %d\n", vert->getNum());
}
}
}
@@ -764,11 +763,13 @@ void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
//std::map<MVertex*,SVector3> verticesToMove;
std::set<MVertex*> verticesToMove;
+ printf(" size = %d\n",_straightSidedLocation.size());
+
for (unsigned int i = 0; i < cavity.size(); i++){
for (int j = 0; j < cavity[i]->getNumVertices(); j++){
MVertex *vert = cavity[i]->getVertex(j);
its = _straightSidedLocation.find(vert);
- if (its != _straightSidedLocation.end()) {
+ if (its == _straightSidedLocation.end()) {
printf("SETTING LOCATIONS for %d\n",vert->getNum());
_straightSidedLocation[vert] =
SVector3(vert->x(), vert->y(), vert->z());
@@ -781,7 +782,7 @@ void highOrderSmoother::smooth_cavity(std::vector<MElement*>& cavity,
if (vert->onWhat()->dim() < _dim){
myAssembler.fixVertex(vert, 0, getTag(), 0);
myAssembler.fixVertex(vert, 1, getTag(), 0);
- printf("Fixing vertex %d\n", vert->getNum());
+ printf("Fixing vertex (boundary) %d\n", vert->getNum());
}
}
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index e9fd4023b8..cd518f885a 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1289,6 +1289,20 @@ void meshGFace::operator() (GFace *gf)
Msg::Debug("Type %d %d triangles generated, %d internal vertices",
gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size());
+
+ // test : recompute the background mesh using a PDE
+ /*
+ if (backgroundMesh::current()){
+ backgroundMesh::unset();
+ }
+ else{
+ backgroundMesh::set(gf);
+ char name[256];
+ sprintf(name,"bgm-%d.pos",gf->tag());
+ backgroundMesh::current()->print(name,gf);
+ (*this)(gf);
+ }
+ */
}
bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges)
@@ -1409,7 +1423,7 @@ void partitionAndRemesh(GFaceCompound *gf)
GFace *gfc = gf->model()->getFaceByTag(numf + NF + i );
meshGFace mgf;
mgf(gfc);
- //lloyd(gfc);
+ // gfc->lloyd(20,0);
for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
MTriangle *t = gfc->triangles[j];
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 93194e0ba0..85c4e19626 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -450,44 +450,6 @@ void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
}
}
-void splitEdgePassUnsorted(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
-{
- int NN1 = m.edges.size();
- int NN2 = 0;
- std::list<BDS_Edge*>::iterator it = m.edges.begin();
- while (1){
- if (NN2++ >= NN1) break;
- if (!(*it)->deleted){
- double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
- if ((*it)->numfaces() == 2 && (lone > MAXE_)){
- const double coord = 0.5;
- //const double coord = computeEdgeMiddleCoord((*it)->p1, (*it)->p2, gf,
- // m.scalingU, m.scalingV);
- BDS_Point *mid;
-
- GPoint gpp = gf->point
- (m.scalingU*(coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u),
- m.scalingV*(coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v));
- if (gpp.succeeded()){
- mid = m.add_point
- (++m.MAXPOINTNUMBER,
- coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
- coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v, gf);
- mid->lcBGM() = BGM_MeshSize
- (gf,
- (coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u)*m.scalingU,
- (coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
- mid->X,mid->Y,mid->Z);
- mid->lc() = 0.5 * ((*it)->p1->lc() + (*it)->p2->lc());
- if(!m.split_edge(*it, mid)) m.del_point(mid);
- else nb_split++;
- }
- }
- }
- ++it;
- }
-}
-
// A test for spheres
static void midpointsphere(GFace *gf, double u1, double v1, double u2,
double v2, double &u12, double &v12, double r)
@@ -560,13 +522,22 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
mid = m.add_point(++m.MAXPOINTNUMBER, gpp.x(),gpp.y(),gpp.z());
mid->u = U;
mid->v = V;
- mid->lcBGM() = BGM_MeshSize
- (gf,
- (coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
- (coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
- mid->X,mid->Y,mid->Z);
- mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
- if(!m.split_edge(e, mid)) m.del_point(mid);
+ if (backgroundMesh::current()){
+ mid->lc() = mid->lcBGM() =
+ backgroundMesh::current()->operator() (
+ (coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
+ (coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
+ 0.0);
+ }
+ else {
+ mid->lcBGM() = BGM_MeshSize
+ (gf,
+ (coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
+ (coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
+ mid->X,mid->Y,mid->Z);
+ mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
+ }
+ if(!m.split_edge(e, mid)) m.del_point(mid);
else nb_split++;
}
}
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index e3a0cbfb40..487f89ed1f 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -7,10 +7,12 @@
#include "MTriangle.h"
#include "Context.h"
#include "meshGFace.h"
+#include "BackgroundMesh.h"
void lloydAlgorithm::operator () ( GFace * gf) {
std::set<MVertex*> all;
+
// get all the points of the face ...
for (unsigned int i = 0; i < gf->getNumMeshElements(); i++){
@@ -20,6 +22,8 @@ void lloydAlgorithm::operator () ( GFace * gf) {
}
}
+ backgroundMesh::set(gf);
+ backgroundMesh::current()->print("bgm.pos",0);
// Create a triangulator
DocRecord triangulator (all.size());
@@ -39,8 +43,8 @@ void lloydAlgorithm::operator () ( GFace * gf) {
SPoint2 p;
bool success = reparamMeshVertexOnFace(*it, gf, p);
if (!success) {
- printf("loyd no succes exit \n");
- exit(1);
+ Msg::Error("Impossible to apply Lloyd to model face %d",gf->tag());
+ Msg::Error("A mesh vertex cannot be reparametrized");
return;
}
double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
@@ -75,7 +79,7 @@ void lloydAlgorithm::operator () ( GFace * gf) {
double E;
SPoint2 p(pt.where.h,pt.where.v);
if (!infiniteNorm){
- centroidOfPolygon (p,pts, cgs(i,0),cgs(i,1),E);
+ centroidOfPolygon (p,pts, cgs(i,0),cgs(i,1),E,backgroundMesh::current());
}
else {
centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E);
@@ -89,7 +93,7 @@ void lloydAlgorithm::operator () ( GFace * gf) {
if (ITER % 10 == 0){
char name[234];
sprintf(name,"LloydIter%d.pos",ITER);
- // triangulator.makePosView(name);
+ triangulator.makePosView(name);
}
for(PointNumero i = 0; i < triangulator.numPoints; i++) {
@@ -109,7 +113,7 @@ void lloydAlgorithm::operator () ( GFace * gf) {
// now create the vertices
std::vector<MVertex*> mesh_vertices;
for (int i=0; i<triangulator.numPoints;i++){
- // get the ith v2ertex
+ // get the ith vertex
PointRecord &pt = triangulator.points[i];
MVertex *v = (MVertex*)pt.data;
if (v->onWhat() == gf && !triangulator.onHull(i)){
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index ae5f1dd217..11fae3321b 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -617,16 +617,21 @@ void centroidOfOrientedBox(std::vector<SPoint2> &pts, const double &angle,
}
void centroidOfPolygon(SPoint2 &pc, std::vector<SPoint2> &pts,
- double &xc, double &yc, double &inertia)
+ double &xc, double &yc, double &inertia,
+ simpleFunction<double> *bgm)
{
double area_tot = 0;
SPoint2 center(0,0);
for (unsigned int j = 0; j < pts.size(); j++){
SPoint2 &pa = pts[j];
SPoint2 &pb = pts[(j+1)%pts.size()];
- const double area = triangle_area2d(pa,pb,pc);
- area_tot += area;
- center += ((pa+pb+pc) * (area/3.0));
+ const double area = triangle_area2d(pa,pb,pc);
+ const double lc = bgm ? (*bgm)((pa.x()+pb.x()+pc.x())/3.0,
+ (pa.y()+pb.y()+pc.y())/3.0,
+ 0.0) : 1.0;
+ const double fact = 1./(lc*lc);
+ area_tot += area*fact;
+ center += ((pa+pb+pc) * (area*fact/3.0));
}
SPoint2 x = center * (1.0/area_tot);
inertia = 0;
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index 90f04276d4..92c58b47a3 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -9,6 +9,7 @@
#include <algorithm>
#include "fullMatrix.h"
#include "SPoint2.h"
+#include "simpleFunction.h"
class binding;
class GFace;
@@ -84,8 +85,8 @@ class DocRecord{
int numTriangles; // number of triangles
Triangle *triangles; // 2D results
DocRecord(int n);
- double &x(int i){ return points[i].where.v; }
- double &y(int i){ return points[i].where.h; }
+ double &x(int i){ return points[i].where.h; }
+ double &y(int i){ return points[i].where.v; }
void* &data(int i){ return points[i].data; }
void setPoints(fullMatrix<double> *p);
~DocRecord();
@@ -108,7 +109,8 @@ void centroidOfPolygon(SPoint2 &pc,
std::vector<SPoint2> &pts,
double &xc,
double &yc,
- double &inertia);
+ double &inertia,
+ simpleFunction<double> *bgm = 0);
#endif
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 8a83381671..6de8e415a0 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -12,12 +12,6 @@
#include "dgTransformNodalValue.h"
#include "meshPartition.h"
-/*
- compute
- \int \vec{f} \cdot \grad \phi dv
-*/
-
-
// works for any number of groups
/*
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index e5bb4c9434..ce37b69cd6 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -16,6 +16,12 @@ class dgSystemOfEquations;
class dgAlgorithm {
public :
+ static void jacobianVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
+ const dgConservationLaw &claw, // the conservation law
+ const dgGroupOfElements & group,
+ fullMatrix<double> &solution,
+ fullMatrix<double> &residual);
+
static void computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group, // the group
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index e4484d04c1..248303d1e5 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -116,6 +116,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
}
// redistribution matrix
// quadrature weight x parametric gradients in quadrature points
+ // redistribition Jacobian Matrix : quadrature weight x parametric gradients x shapefunctions
_PsiDPsiDXi = fullMatrix<double> (_dimUVW*nbQP, nbNodes*nbNodes);
//reditribution of the diffusive jacobian dimUVW*dimUVW*nbIntegrationPoints x nbNodes*nbNodes
@@ -135,6 +136,11 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
(*_redistributionFluxes[2])(k,xi) = g[k][2] * weight;
(*_redistributionSource)(k,xi) = f[k] * weight;
(*_collocation)(xi,k) = f[k];
+ // for (int l=0;l<_fs.coefficients.size1();l++){
+ // (*_redistributionJacobianOfFluxes[0])(l+_fs.coefficients.size1()*k,j) = g[k][0] * f[l] * weight;
+ // (*_redistributionJacobianOfFluxes[0])(l+_fs.coefficients.size1()*k,j) = g[k][1] * f[l] * weight;
+ // (*_redistributionJacobianOfFluxes[0])(l+_fs.coefficients.size1()*k,j) = g[k][2] * f[l] * weight;
+ // }
}
for (int alpha=0; alpha<_dimUVW; alpha++) for (int beta=0; beta<_dimUVW; beta++) {
for (int i=0; i<nbNodes; i++) for (int j=0; j<nbNodes; j++) {
@@ -157,6 +163,9 @@ void dgGroupOfElements::copyPrivateDataFrom(const dgGroupOfElements *from){
dgGroupOfElements::~dgGroupOfElements(){
delete _integration;
+ // delete _redistributionJacobianOfFluxes[0];
+ // delete _redistributionJacobianOfFluxes[1];
+ // delete _redistributionJacobianOfFluxes[2];
delete _redistributionFluxes[0];
delete _redistributionFluxes[1];
delete _redistributionFluxes[2];
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index ed2ea0dc4d..9a8a25b0e1 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -67,6 +67,8 @@ class dgGroupOfElements {
// the gradient of shape functions (in parametric space)
// for both diffusive and convective fluxes
fullMatrix<double> *_redistributionFluxes[3];
+ // redistribution of the jacobian of fluxes to vertices
+ // fullMatrix<double> *_redistributionJacobianOfFluxes[3];
// redistribution for the source term
fullMatrix<double> *_redistributionSource;
// inverse mass matrix of all elements
@@ -111,6 +113,7 @@ public:
inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;}
inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
inline const fullMatrix<double> & getFluxRedistributionMatrix (int i) const {return *_redistributionFluxes[i];}
+ // inline const fullMatrix<double> & getJacobianOfFluxRedistributionMatrix (int i) const {return *_redistributionJacobianOfFluxes[i];}
inline const fullMatrix<double> & getSourceRedistributionMatrix () const {return *_redistributionSource;}
inline const fullMatrix<double> & getDPsiDXDPsiDXi() const {return _dPsiDXDPsiDXi;}
inline const fullMatrix<double> & getPsiDPsiDXi() const {return _PsiDPsiDXi;}
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index f3b52fb66f..7bda156be0 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -72,7 +72,7 @@ class helmholtzTerm : public femTerm<scalar> {
const double w = GP[i].pt[2];
const double weightDetJ = GP[i].weight * e->getJacobian(u, v, w, jac);
SPoint3 p; e->pnt(u, v, w, p);
- const scalar K = (*_k)(p.x(), p.y(), p.z());
+ const scalar K = _k ? (*_k)(p.x(), p.y(), p.z()) : 0.0;
const scalar A = _a ? (*_a)(p.x(), p.y(), p.z()) : 0.0;
inv3x3(jac, invjac) ;
e->getGradShapeFunctions(u, v, w, grads);
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 3f1f221eda..207a1060ac 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/Solver/linearSystemCSR.cpp
@@ -276,7 +276,7 @@ void sortColumns_(int NbLines,
for (int j = jptr[i]; j < jptr[i + 1] - 1; j++){
ptr[j] = j + 1;
}
- ptr[jptr[i + 1]] = 0;
+ // ptr[jptr[i + 1]] = 0;
}
delete[] count;
@@ -351,6 +351,7 @@ extern "C" {
template<>
int linearSystemCSRTaucs<double>::systemSolve()
{
+ if (!_a)return 1;
if(!sorted){
sortColumns_(_b->size(),
CSRList_Nbr(_a),
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 953dcdc192..c8989a61ab 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -24,10 +24,31 @@
//--------------------------------------------------------------
+struct compareRotatedPoints {
+ double angle;
+ const SPoint2 &left;
+ compareRotatedPoints (const SPoint2& l,const SPoint2& r) : left(l){
+ angle = atan2(r.y()-l.y(),r.x()-l.x());
+ }
+ bool operator ( ) (const SPoint2 &p1, const SPoint2 &p2) const {
+ double x1 = (p1.x()-left.x())*cos(angle) + (p1.y()-left.y())*sin(angle);
+ double x2 = (p2.x()-left.x())*cos(angle) + (p2.y()-left.y())*sin(angle);
+ if (x1<x2)return true;
+ if (x1>x2)return false;
+ double y1 =-(p1.x()-left.x())*sin(angle) + (p1.y()-left.y())*cos(angle);
+ double y2 =-(p2.x()-left.x())*sin(angle) + (p2.y()-left.y())*cos(angle);
+ if (y1<y2)return true;
+ return false;
+ }
+};
+
+
struct sort_pred {
- bool operator()(const std::pair<SPoint2,multiscaleLaplaceLevel*> &left, const std::pair<SPoint2,multiscaleLaplaceLevel*> &right) {
- return left.first.x() < right.first.x();
- }
+ compareRotatedPoints comparator;
+ sort_pred (const SPoint2 &left, const SPoint2 &right) : comparator(left,right) {}
+ bool operator()(const std::pair<SPoint2,multiscaleLaplaceLevel*> &left, const std::pair<SPoint2,multiscaleLaplaceLevel*> &right) {
+ return comparator(left.first,right.first);
+ }
};
@@ -214,6 +235,8 @@ static int intersection_segments (SPoint2 &p1, SPoint2 &p2,
}
//--------------------------------------------------------------
+
+
static void recur_compute_centers_ (double R, double a1, double a2,
multiscaleLaplaceLevel * root, int &nbElems ){
@@ -279,7 +302,7 @@ static void recur_compute_centers_ (double R, double a1, double a2,
}
//sort centers
- std::sort(centers.begin(),centers.end(), sort_pred());
+ std::sort(centers.begin(),centers.end(), sort_pred(PL,PR));
int nbrecur = 0;
for (int i=1;i<centers.size()-1;i++){
@@ -294,7 +317,7 @@ static void recur_compute_centers_ (double R, double a1, double a2,
printf("children =%d nbrecur =%d \n", root->children.size(), nbrecur);
for (int i=0;i<centers.size();i++){
printf("center i=%d (%g %g)\n", i, centers[i].first.x(), centers[i].first.x());
- if(centers[i].second) printf(" level %d recur %d elems =%d\n", centers[i].second->recur, centers[i].second->region );
+ if(centers[i].second) printf(" level %d recur %d elems =%d\n", centers[i].second->recur, centers[i].second->region, centers[i].second->elements.size() );
}
}
@@ -531,18 +554,10 @@ static void recur_cut_ (double R, double a1, double a2,
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++){
@@ -805,6 +820,7 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
root->recur = 0;
root->region = 0;
root->scale = 1.0;
+ root->_name = "Root";
int totNbElems = 0;
parametrize(*root, totNbElems);
@@ -1009,11 +1025,10 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level, int &totNbE
totNbElems += level.elements.size();
//Save multiscale meshes
- char name[245];
- sprintf(name,"multiscale_%d_%d_real.msh",level.recur, level.region);
- printLevel (name,level.elements,0,2.0);
- sprintf(name,"multiscale_%d_%d_param.msh",level.recur, level.region);
- printLevel (name,level.elements,&level.coordinates,2.0);
+ std::string name1(level._name+"real.msh");
+ std::string name2(level._name+"param.msh");
+ printLevel (name1.c_str(),level.elements,0,2.0);
+ printLevel (name2.c_str(),level.elements,&level.coordinates,2.0);
//For every small region compute a new parametrization
Msg::Info("Level (%d-%d): %d connected small regions",level.recur, level.region, regions.size());
@@ -1031,6 +1046,9 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level, int &totNbE
nextLevel->elements = regions[i];
nextLevel->recur = level.recur+1;
nextLevel->region = i;
+ std::stringstream s1 ; s1 << nextLevel->recur;
+ std::stringstream s2 ; s2 << nextLevel->region;
+ nextLevel->_name = level._name+"-"+s1.str()+"-"+s2.str();
SBoundingBox3d smallB;
for(std::set<MVertex *>::iterator itv = tooSmallv.begin(); itv !=tooSmallv.end() ; ++itv){
SPoint2 p = solution[*itv];
@@ -1123,14 +1141,17 @@ void multiscaleLaplace::cut (std::vector<MElement *> &elements)
connected_left_right(left, right);
+ printf("left.size() = %d, right.size() = %d, nbElem = %d elements.size() = %d\n",
+ left.size(),right.size(),totNbElems,elements.size());
+
+ printLevel ("cut-left.msh",left,0,2.2);
+ printLevel ("cut-right.msh",right,0,2.2);
elements.clear();
elements.insert(elements.end(),left.begin(),left.end());
elements.insert(elements.end(),right.begin(),right.end());
- //printLevel ("multiscale-all.msh",elements, 0,2.0);
- //printLevel ("multiscale-left.msh",left,0,2.0);
- //printLevel ("multiscale-right.msh",right,0,2.0);
+ printLevel ("cut-all.msh",elements, 0,2.2);
}
diff --git a/Solver/multiscaleLaplace.h b/Solver/multiscaleLaplace.h
index 7f309fadd2..d587dfd9ec 100644
--- a/Solver/multiscaleLaplace.h
+++ b/Solver/multiscaleLaplace.h
@@ -20,6 +20,7 @@ struct multiscaleLaplaceLevel {
std::vector<MElement *> elements;
std::map<MVertex*,SPoint2> coordinates;
std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > cut;
+ std::string _name;
};
class multiscaleLaplace{
diff --git a/benchmarks/boolean/JAW.geo b/benchmarks/boolean/JAW.geo
index c9d2d566b3..dad88b5f70 100644
--- a/benchmarks/boolean/JAW.geo
+++ b/benchmarks/boolean/JAW.geo
@@ -1,4 +1,4 @@
-Mesh.CharacteristicLengthFactor=0.15;
+//Mesh.CharacteristicLengthFactor=0.15;
L = 100;
@@ -16,4 +16,4 @@ 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};
+//Compound Surface(100) = {1 ... 26};
diff --git a/benchmarks/boolean/TORUS.geo b/benchmarks/boolean/TORUS.geo
index e62b092c1c..bc906c45fd 100644
--- a/benchmarks/boolean/TORUS.geo
+++ b/benchmarks/boolean/TORUS.geo
@@ -1,13 +1,13 @@
Mesh.CharacteristicLengthFactor=0.1;
X = 12;
-For I In {0:2}
+For I In {0:1}
OCCShape("Torus",{2*I*X,0,0,0,0,1,10,4},"Fuse");
OCCShape("Torus",{2*I*X,24,0,0,0,1,10,4},"Fuse");
OCCShape("Torus",{2*I*X,48,0,0,0,1,10,4},"Fuse");
EndFor
-Compound Surface(100) = {1,2,3,4,5,6,7,8,9};
+//Compound Surface(100) = {1,2,3,4,5,6,7,8,9};
diff --git a/benchmarks/step/capot.geo b/benchmarks/step/capot.geo
index e1291d8eeb..d071a1530d 100644
--- a/benchmarks/step/capot.geo
+++ b/benchmarks/step/capot.geo
@@ -1,9 +1,9 @@
-Mesh.CharacteristicLengthFactor=0.2;
+//Mesh.CharacteristicLengthFactor=0.2;
Merge "capot.brep";
Compound Line(1000) = {47,50};
Compound Line(1001) = {44,46};
-Compound Surface(100) = {1, 8, 15, 17, 16, 18, 9, 2, 3, 10, 7, 14, 11, 4, 12, 5, 6, 13} Boundary {{},{},{},{}};
+//Compound Surface(100) = {1, 8, 15, 17, 16, 18, 9, 2, 3, 10, 7, 14, 11, 4, 12, 5, 6, 13} Boundary {{},{},{},{}};
//Recombine Surface {7, 4, 5, 6, 13, 12, 3, 11, 14, 2, 15, 16, 1, 9, 10, 8, 17, 18};
--
GitLab