diff --git a/CMakeLists.txt b/CMakeLists.txt
index c00ecbc275f40978aaa6d71180be477e582207d2..a9004953242f0c20234f16f1f0fbdac943fac7b1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -635,7 +635,7 @@ if(HAVE_MESH)
endif(ENABLE_TETGEN AND EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/contrib/Tetgen1.5/tetgen.h)
if(HAVE_TETGEN)
message("WARNING: By including Tetgen you have to comply with Tetgen's "
- "special licensing requirements stated in contrib/Tetgen/LICENSE.")
+ "special licensing requirements stated in contrib/Tetgen*/LICENSE.")
endif(HAVE_TETGEN)
endif(HAVE_MESH)
diff --git a/Mesh/simple3D.cpp b/Mesh/simple3D.cpp
index c9baa841c7d3bca0005ecc99013ea72edba2bbe4..4fd499293cdbc1d5385813f295ed979b5b564ab6 100755
--- a/Mesh/simple3D.cpp
+++ b/Mesh/simple3D.cpp
@@ -14,9 +14,9 @@
#include <algorithm>
#include "directions3D.h"
-#if defined(HAVE_RTREE)
+#if defined(HAVE_RTREE)
#include "rtree.h"
-#endif
+#endif
#define k1 0.7 //k1*h is the minimal distance between two nodes
#define k2 0.5 //k2*h is the minimal distance to the boundary
@@ -38,14 +38,14 @@ class Metric{
void set_m32(double);
void set_m13(double);
void set_m23(double);
- void set_m33(double);
+ void set_m33(double);
double get_m11();
double get_m21();
double get_m31();
- double get_m12();
+ double get_m12();
double get_m22();
double get_m32();
- double get_m13();
+ double get_m13();
double get_m23();
double get_m33();
};
@@ -93,7 +93,7 @@ double infinity_distance(SPoint3 p1,SPoint3 p2,Metric m){
double x1,y1,z1;
double x2,y2,z2;
double a,b,c,d,e,f,g,h,i;
-
+
a = m.get_m11();
b = m.get_m21();
c = m.get_m31();
@@ -103,16 +103,16 @@ double infinity_distance(SPoint3 p1,SPoint3 p2,Metric m){
g = m.get_m13();
h = m.get_m23();
i = m.get_m33();
-
+
x1 = a*p1.x() + b*p1.y() + c*p1.z();
y1 = d*p1.x() + e*p1.y() + f*p1.z();
z1 = g*p1.x() + h*p1.y() + i*p1.z();
-
+
x2 = a*p2.x() + b*p2.y() + c*p2.z();
y2 = d*p2.x() + e*p2.y() + f*p2.z();
z2 = g*p2.x() + h*p2.y() + i*p2.z();
-
- distance = std::max(std::max(fabs(x2-x1),fabs(y2-y1)),fabs(z2-z1));
+
+ distance = std::max(std::max(fabs(x2-x1),fabs(y2-y1)),fabs(z2-z1));
return distance;
}
@@ -122,13 +122,13 @@ bool rtree_callback(Node* neighbour,void* w){
Metric m;
Node *spawn,*parent;
Wrapper* wrapper;
-
+
wrapper = static_cast<Wrapper*>(w);
spawn = wrapper->get_spawn();
parent = wrapper->get_parent();
h = spawn->get_size();
m = spawn->get_metric();
-
+
if(neighbour!=parent){
distance = infinity_distance(spawn->get_point(),neighbour->get_point(),m);
if(distance<k1*h){
@@ -136,7 +136,7 @@ bool rtree_callback(Node* neighbour,void* w){
return false;
}
}
-
+
return true;
}
@@ -310,9 +310,9 @@ void Filler::treat_model(){
GRegion* gr;
GModel* model = GModel::current();
GModel::riter it;
-
- Frame_field::init_model();
-
+
+ Frame_field::init_model();
+
for(it=model->firstRegion();it!=model->lastRegion();it++)
{
gr = *it;
@@ -320,13 +320,14 @@ void Filler::treat_model(){
treat_region(gr);
}
}
-
+
Frame_field::clear();
}
void Filler::treat_region(GRegion* gr){
-#if defined(HAVE_RTREE)
- int i,j;
+#if defined(HAVE_RTREE)
+ unsigned int i;
+ int j;
int count;
bool ok;
double x,y,z;
@@ -346,7 +347,7 @@ void Filler::treat_region(GRegion* gr){
RTree<Node*,double,3,double> rtree;
octree = new MElementOctree(gr->model());
-
+
for(i=0;i<gr->getNumMeshElements();i++){
element = gr->getMeshElement(i);
for(j=0;j<element->getNumVertices();j++){
@@ -354,7 +355,7 @@ void Filler::treat_region(GRegion* gr){
old_vertices.insert(vertex);
}
}
-
+
for(it=old_vertices.begin();it!=old_vertices.end();it++){
if((*it)->onWhat()->dim()<3){
boundary_vertices.push_back(*it);
@@ -407,7 +408,7 @@ void Filler::treat_region(GRegion* gr){
wrapper.set_too_close(0);
wrapper.set_spawn(spawn);
wrapper.set_parent(parent);
- rtree.Search(spawn->min,spawn->max,rtree_callback,&wrapper);
+ rtree.Search(spawn->min,spawn->max,rtree_callback,&wrapper);
if(!wrapper.get_too_close()){
fifo.push(spawn);
rtree.Insert(spawn->min,spawn->max,spawn);
@@ -422,14 +423,14 @@ void Filler::treat_region(GRegion* gr){
printf("%d\n",count);
count++;
}
-
+
deleter(gr);
std::vector<GRegion*> regions;
regions.push_back(gr);
meshGRegion mesher(regions); //?
mesher(gr); //?
MeshDelaunayVolume(regions);
-
+
delete octree;
for(i=0;i<garbage.size();i++) delete garbage[i];
for(i=0;i<new_vertices.size();i++) delete new_vertices[i];
@@ -439,22 +440,22 @@ void Filler::treat_region(GRegion* gr){
Metric Filler::get_metric(double x,double y,double z){
Metric m;
- Matrix m2;
-
- m2 = Frame_field::search(x,y,z);
-
+ Matrix m2;
+
+ m2 = Frame_field::search(x,y,z);
+
m.set_m11(m2.get_m11());
m.set_m21(m2.get_m21());
m.set_m31(m2.get_m31());
-
+
m.set_m12(m2.get_m12());
m.set_m22(m2.get_m22());
m.set_m32(m2.get_m32());
-
+
m.set_m13(m2.get_m13());
m.set_m23(m2.get_m23());
m.set_m33(m2.get_m33());
-
+
return m;
}
@@ -467,28 +468,28 @@ Metric Filler::get_clf_metric(double x,double y,double z,GEntity* ge){
Metric m;
Field* field;
FieldManager* manager;
-
+
m = Metric();
manager = ge->model()->getFields();
if(manager->getBackgroundField()>0){
field = manager->get(manager->getBackgroundField());
if(field){
(*field)(x,y,z,v1,v2,v3,ge);
-
- m.set_m11(v1.x());
+
+ m.set_m11(v1.x());
m.set_m21(v1.y());
m.set_m31(v1.z());
-
+
m.set_m12(v2.x());
m.set_m22(v2.y());
m.set_m32(v2.z());
-
+
m.set_m13(v3.x());
m.set_m23(v3.y());
m.set_m33(v3.z());
}
}
-
+
return m;
}
@@ -496,8 +497,8 @@ double Filler::get_clf_size(double x,double y,double z,GEntity* ge){
double h;
Field* field;
FieldManager* manager;
-
- h = 0.25;
+
+ h = 0.25;
manager = ge->model()->getFields();
if(manager->getBackgroundField()>0){
field = manager->get(manager->getBackgroundField());
@@ -505,7 +506,7 @@ double Filler::get_clf_size(double x,double y,double z,GEntity* ge){
h = (*field)(x,y,z,ge);
}
}
-
+
return h;
}
@@ -521,20 +522,20 @@ bool Filler::far_from_boundary(MElementOctree* octree,Node* node){
double h;
SPoint3 point;
MElement *e1,*e2,*e3,*e4,*e5,*e6;
-
+
point = node->get_point();
x = point.x();
y = point.y();
z = point.z();
h = node->get_size();
-
+
e1 = (MElement*)octree->find(x+k2*h,y,z,3,true);
e2 = (MElement*)octree->find(x-k2*h,y,z,3,true);
e3 = (MElement*)octree->find(x,y+k2*h,z,3,true);
e4 = (MElement*)octree->find(x,y-k2*h,z,3,true);
e5 = (MElement*)octree->find(x,y,z+k2*h,3,true);
e6 = (MElement*)octree->find(x,y,z-k2*h,3,true);
-
+
if(e1!=NULL && e2!=NULL && e3!=NULL && e4!=NULL && e5!=NULL && e6!=NULL) return 1;
else return 0;
}
@@ -544,14 +545,14 @@ void Filler::compute_parameters(Node* node,GEntity* ge){
double h;
Metric m;
SPoint3 point;
-
+
point = node->get_point();
x = point.x();
y = point.y();
z = point.z();
m = get_metric(x,y,z);
h = get_size(x,y,z);
-
+
node->set_size(h);
node->set_metric(m);
node->min[0] = x - sqrt3*h;
@@ -574,44 +575,44 @@ void Filler::offsprings(GEntity* ge,MElementOctree* octree,Node* node,Node* n1,N
double h1,h2,h3,h4,h5,h6;
Metric m;
SPoint3 point;
-
+
point = node->get_point();
x = point.x();
y = point.y();
z = point.z();
h = node->get_size();
m = node->get_metric();
-
+
h1 = improvement(ge,octree,point,h,SVector3(m.get_m11(),m.get_m21(),m.get_m31()));
x1 = x + h1*m.get_m11();
y1 = y + h1*m.get_m21();
z1 = z + h1*m.get_m31();
-
+
h2 = improvement(ge,octree,point,h,SVector3(-m.get_m11(),-m.get_m21(),-m.get_m31()));
x2 = x - h2*m.get_m11();
y2 = y - h2*m.get_m21();
z2 = z - h2*m.get_m31();
-
+
h3 = improvement(ge,octree,point,h,SVector3(m.get_m12(),m.get_m22(),m.get_m32()));
x3 = x + h3*m.get_m12();
y3 = y + h3*m.get_m22();
z3 = z + h3*m.get_m32();
-
+
h4 = improvement(ge,octree,point,h,SVector3(-m.get_m12(),-m.get_m22(),-m.get_m32()));
x4 = x - h4*m.get_m12();
y4 = y - h4*m.get_m22();
z4 = z - h4*m.get_m32();
-
+
h5 = improvement(ge,octree,point,h,SVector3(m.get_m13(),m.get_m23(),m.get_m33()));
x5 = x + h5*m.get_m13();
y5 = y + h5*m.get_m23();
z5 = z + h5*m.get_m33();
-
+
h6 = improvement(ge,octree,point,h,SVector3(-m.get_m13(),-m.get_m23(),-m.get_m33()));
x6 = x - h6*m.get_m13();
y6 = y - h6*m.get_m23();
z6 = z - h6*m.get_m33();
-
+
*n1 = Node(SPoint3(x1,y1,z1));
*n2 = Node(SPoint3(x2,y2,z2));
*n3 = Node(SPoint3(x3,y3,z3));
@@ -625,7 +626,7 @@ double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,doub
double average;
double h_farther;
double coeffA,coeffB;
-
+
x = point.x() + h_nearer*direction.x();
y = point.y() + h_nearer*direction.y();
z = point.z() + h_nearer*direction.z();
@@ -633,7 +634,7 @@ double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,doub
h_farther = get_size(x,y,z);
}
else h_farther = h_nearer;
-
+
coeffA = 1.0;
coeffB = 0.2;
if(h_farther>h_nearer){
@@ -648,16 +649,16 @@ double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,doub
int Filler::get_nbr_new_vertices(){
return new_vertices.size();
}
-
+
MVertex* Filler::get_new_vertex(int i){
return new_vertices[i];
}
void Filler::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file){
- file << "SL ("
+ file << "SL ("
<< p1.x() << ", " << p1.y() << ", " << p1.z() << ", "
- << p2.x() << ", " << p2.y() << ", " << p2.z() << ")"
- << "{10, 20};\n";
+ << p2.x() << ", " << p2.y() << ", " << p2.z() << ")"
+ << "{10, 20};\n";
}
void Filler::print_node(Node* node,std::ofstream& file){
@@ -671,34 +672,34 @@ void Filler::print_node(Node* node,std::ofstream& file){
double h;
Metric m;
SPoint3 point;
-
+
point = node->get_point();
x = point.x();
y = point.y();
z = point.z();
h = node->get_size();
m = node->get_metric();
-
+
x1 = x + h*m.get_m11();
y1 = y + h*m.get_m21();
z1 = z + h*m.get_m31();
-
+
x2 = x - h*m.get_m11();
y2 = y - h*m.get_m21();
z2 = z - h*m.get_m31();
-
+
x3 = x + h*m.get_m12();
y3 = y + h*m.get_m22();
z3 = z + h*m.get_m32();
-
+
x4 = x - h*m.get_m12();
y4 = y - h*m.get_m22();
z4 = z - h*m.get_m32();
-
+
x5 = x + h*m.get_m13();
y5 = y + h*m.get_m23();
z5 = z + h*m.get_m33();
-
+
x6 = x - h*m.get_m13();
y6 = y - h*m.get_m23();
z6 = z - h*m.get_m33();
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index 079ce7a1686b4005fc2d2a2419ce781f28f3cf6f..40342e57d243f226d4251a57b42fc6970272740a 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -43,7 +43,7 @@ extern "C"
}
}
-GMSH_DistancePlugin::GMSH_DistancePlugin()
+GMSH_DistancePlugin::GMSH_DistancePlugin()
{
_fileName = "text";
_minScale = 0.0;
@@ -91,28 +91,29 @@ StringXString *GMSH_DistancePlugin::getOptionStr(int iopt)
}
void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
- std::map<MVertex*,double > _distance_map )
+ std::map<MVertex*, double > _distance_map)
{
_fileName = DistanceOptions_String[0].def;
_minScale = (double) DistanceOptions_Number[4].def;
_maxScale = (double) DistanceOptions_Number[5].def;
-
+
double minDist=1.e4;
double maxDist=0.0;
- for (std::map<MVertex*,double >::iterator itv=_distance_map.begin(); itv!=_distance_map.end(); ++itv){
+ for (std::map<MVertex*,double >::iterator itv=_distance_map.begin();
+ itv != _distance_map.end(); ++itv){
double dist = itv->second;
if (dist>maxDist) maxDist = dist;
if (dist<minDist) minDist = dist;
itv->second = dist;
}
-
+
Msg::Info("Writing %s", _fileName.c_str());
FILE *fName = fopen(_fileName.c_str(),"w");
fprintf(fName, "View \"distance \"{\n");
-
+
for (unsigned int ii=0; ii<_entities.size(); ii++) {
if (_entities[ii]->dim() == _maxDim) {
- for (unsigned int i=0; i<_entities[ii]->getNumMeshElements(); i++) {
+ for (unsigned int i=0; i<_entities[ii]->getNumMeshElements(); i++) {
MElement *e = _entities[ii]->getMeshElement(i);
int numNodes = e->getNumVertices();
if (e->getNumChildren())
@@ -120,7 +121,7 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
std::vector<double> x(numNodes), y(numNodes), z(numNodes);
std::vector<double> *out = _data->incrementList(1, e->getType(), numNodes);
std::vector<MVertex*> nods;
-
+
if (!e->getNumChildren())
for(int i = 0; i < numNodes; i++)
nods.push_back(e->getVertex(i));
@@ -128,17 +129,18 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
for(int i = 0; i < e->getNumChildren(); i++)
for(int j = 0; j < e->getChild(i)->getNumVertices(); j++)
nods.push_back(e->getChild(i)->getVertex(j));
-
- for (int nod=0; nod<numNodes; nod++) out->push_back((nods[nod])->x());
- for (int nod=0; nod<numNodes; nod++) out->push_back((nods[nod])->y());
- for (int nod=0; nod<numNodes; nod++) out->push_back((nods[nod])->z());
-
+
+ for (int nod = 0; nod < numNodes; nod++) out->push_back((nods[nod])->x());
+ for (int nod = 0; nod < numNodes; nod++) out->push_back((nods[nod])->y());
+ for (int nod = 0; nod < numNodes; nod++) out->push_back((nods[nod])->z());
+
if (_maxDim == 2)
switch (numNodes) {
case 2: fprintf(fName,"SL("); break;
case 3: fprintf(fName,"ST("); break;
case 4: fprintf(fName,"SQ("); break;
- default: Msg::Error("Error in Plugin 'Distance' (numNodes=%d).",numNodes); break;
+ default: Msg::Error("Error in Plugin 'Distance' (numNodes=%d)",
+ numNodes); break;
}
else if (_maxDim == 3)
switch (numNodes) {
@@ -146,9 +148,10 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
case 8: fprintf(fName,"SH("); break;
case 6: fprintf(fName,"SI("); break;
case 5: fprintf(fName,"SY("); break;
- default: Msg::Error("Error in Plugin 'Distance' (numNodes=%d).",numNodes); break;
+ default: Msg::Error("Error in Plugin 'Distance' (numNodes=%d)",
+ numNodes); break;
}
-
+
std::vector<double> dist;
for (int j=0; j<numNodes; j++) {
MVertex *v = nods[j];
@@ -159,13 +162,14 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
std::map<MVertex*, double>::iterator it = _distance_map.find(v);
dist.push_back(it->second);
}
-
+
fprintf(fName,"){");
for (unsigned int i=0; i<dist.size(); i++) {
if (_minScale>0 && _maxScale>0)
- dist[i] = _minScale+((dist[i]-minDist)/(maxDist-minDist))*(_maxScale-_minScale);
+ dist[i] = _minScale + ((dist[i] - minDist) / (maxDist - minDist)) *
+ (_maxScale - _minScale);
else if (_minScale>0 && _maxScale<0)
- dist[i] = _minScale+dist[i];
+ dist[i] = _minScale + dist[i];
out->push_back(dist[i]);
if (i)
fprintf(fName, ",%.16g", dist[i]);
@@ -173,12 +177,12 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
fprintf(fName, "%.16g", dist[i]);
}
fprintf(fName,"};\n");
-
+
}
}
}
fprintf(fName,"};\n");
-
+
fclose(fName);
}
@@ -189,7 +193,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
int id_face = (int) DistanceOptions_Number[2].def;
double type = (double) DistanceOptions_Number[3].def;
int ortho = (int) DistanceOptions_Number[6].def;
-
+
PView *view = new PView();
_data = getDataList(view);
#if defined(HAVE_SOLVER)
@@ -203,17 +207,17 @@ PView *GMSH_DistancePlugin::execute(PView *v)
#endif
dofManager<double> * dofView = new dofManager<double>(lsys);
#endif
-
+
std::vector<GEntity*> _entities;
GModel::current()->getEntities(_entities);
if (!_entities.size() || !_entities[_entities.size()-1]->getMeshElement(0)) {
Msg::Error("This plugin needs a mesh !");
return view;
}
-
+
GEntity* ge = _entities[_entities.size()-1];
int integrationPointTetra[2] = {0,0};
-
+
int numnodes = 0;
for (unsigned int i = 0; i < _entities.size()-1; i++)
numnodes += _entities[i]->mesh_vertices.size();
@@ -224,13 +228,13 @@ PView *GMSH_DistancePlugin::execute(PView *v)
std::vector<SPoint3> pts;
std::vector<double> distances;
std::vector<MVertex* > pt2Vertex;
- pts.clear();
+ pts.clear();
distances.clear();
pt2Vertex.clear();
pts.reserve(totNumNodes);
distances.reserve(totNumNodes);
pt2Vertex.reserve(totNumNodes);
-
+
std::map<MVertex*,double> _distanceE_map;
std::map<MVertex*,int> _isInYarn_map;
std::vector<int> index;
@@ -241,11 +245,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
std::vector<int> isInYarn2;
std::vector<SPoint3> closePts;
std::vector<SPoint3> closePts2;
-
+
for (int i=0; i<totNumNodes; i++) {
distances.push_back(1.e22);
}
-
+
int k = 0;
for (unsigned int i=0; i<_entities.size(); i++){
GEntity* ge = _entities[i];
@@ -262,29 +266,29 @@ PView *GMSH_DistancePlugin::execute(PView *v)
k++;
}
}
-
+
// Compute geometrical distance to mesh boundaries
//------------------------------------------------------
if (type < 0.0 ) {
-
+
bool existEntity = false;
-
+
for (unsigned int i=0; i<_entities.size(); i++) {
GEntity* g2 = _entities[i];
- int tag = g2->tag();
int gDim = g2->dim();
std::vector<int> phys = g2->getPhysicalEntities();
bool computeForEntity = false;
- for(int k = 0; k<phys.size(); k++) {
+ for(unsigned int k = 0; k<phys.size(); k++) {
int tagp = phys[k];
- if (id_pt==0 && id_line==0 && id_face==0 && gDim==_maxDim-1 )
+ if (id_pt == 0 && id_line == 0 && id_face == 0 && gDim == _maxDim - 1)
computeForEntity = true;
- else if ( (tagp==id_pt && gDim==0) || (tagp==id_line && gDim==1) || (tagp==id_face && gDim==2) )
+ else if ((tagp == id_pt && gDim == 0) || (tagp == id_line && gDim == 1) ||
+ (tagp == id_face && gDim == 2))
computeForEntity = true;
}
if (computeForEntity) {
existEntity = true;
- for (unsigned int k=0; k<g2->getNumMeshElements(); k++) {
+ for (unsigned int k = 0; k < g2->getNumMeshElements(); k++) {
std::vector<double> iDistances;
std::vector<SPoint3> iClosePts;
std::vector<double> iDistancesE;
@@ -294,10 +298,12 @@ PView *GMSH_DistancePlugin::execute(PView *v)
MVertex *v2 = e->getVertex(1);
SPoint3 p1(v1->x(), v1->y(), v1->z());
SPoint3 p2(v2->x(), v2->y(), v2->z());
- if ((e->getNumVertices() == 2 && order==1) || (e->getNumVertices() == 3 && order==2)) {
- signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
+ if ((e->getNumVertices() == 2 && order == 1) ||
+ (e->getNumVertices() == 3 && order == 2)) {
+ signedDistancesPointsLine(iDistances, iClosePts, pts, p1, p2);
}
- else if ((e->getNumVertices() == 3 && order == 1) || (e->getNumVertices() == 6 && order==2)) {
+ else if ((e->getNumVertices() == 3 && order == 1) ||
+ (e->getNumVertices() == 6 && order == 2)) {
MVertex *v3 = e->getVertex(2);
SPoint3 p3 (v3->x(),v3->y(),v3->z());
signedDistancesPointsTriangle(iDistances, iClosePts, pts, p1, p2, p3);
@@ -307,9 +313,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
distances[kk] = std::abs(iDistances[kk]);
MVertex *v = pt2Vertex[kk];
_distance_map[v] = distances[kk];
-/* TO DO (by AM)
+ /* TO DO (by AM)
_closePts_map[v] = iClosePts[kk];
-*/
+ */
}
}
}
@@ -321,38 +327,38 @@ PView *GMSH_DistancePlugin::execute(PView *v)
if (id_face != 0) Msg::Error("The Physical Surface does not exist !");
return view;
}
-
+
printView(_entities, _distance_map);
-
-/* TO DO (by AM)
- printView(_entities, _closePts_map);
-*/
+
+ /* TO DO (by AM)
+ printView(_entities, _closePts_map);
+ */
}
-
+
// Compute PDE for distance function
//-----------------------------------
else if (type > 0.0) {
-
+
#if defined(HAVE_SOLVER)
-
+
bool existEntity = false;
SBoundingBox3d bbox;
for(unsigned int i = 0; i < _entities.size(); i++){
GEntity* ge = _entities[i];
- int tag = ge->tag();
int gDim = ge->dim();
bool fixForEntity = false;
std::vector<int> phys = ge->getPhysicalEntities();
- for(int k=0; k<phys.size(); k++) {
+ for(unsigned int k = 0; k < phys.size(); k++) {
int tagp = phys[k];
- if (id_pt==0 && id_line==0 && id_face==0 && gDim==_maxDim-1 )
+ if (id_pt == 0 && id_line == 0 && id_face == 0 && gDim == _maxDim - 1)
fixForEntity = true;
- else if ( (tagp==id_pt && gDim==0) || (tagp==id_line && gDim==1) || (tagp==id_face && gDim==2) )
+ else if ((tagp == id_pt && gDim == 0) || (tagp == id_line && gDim == 1) ||
+ (tagp == id_face && gDim == 2) )
fixForEntity = true;
}
if (fixForEntity) {
existEntity = true;
- for (unsigned int i=0; i<ge->getNumMeshElements(); ++i) {
+ for (unsigned int i = 0; i < ge->getNumMeshElements(); ++i) {
MElement *t = ge->getMeshElement(i);
for (int k=0; k<t->getNumVertices(); k++) {
MVertex *v = t->getVertex(k);
@@ -362,46 +368,47 @@ PView *GMSH_DistancePlugin::execute(PView *v)
}
}
}
-
+
if (!existEntity){
if (id_pt != 0) Msg::Error("The Physical Point does not exist !");
if (id_line != 0) Msg::Error("The Physical Line does not exist !");
if (id_face != 0) Msg::Error("The Physical Surface does not exist !");
return view;
}
-
+
std::vector<MElement *> allElems;
- for(unsigned int ii=0; ii<_entities.size(); ii++){
+ for(unsigned int ii = 0; ii < _entities.size(); ii++){
if(_entities[ii]->dim() == _maxDim) {
GEntity *ge = _entities[ii];
- for(unsigned int i=0; i<ge->getNumMeshElements(); ++i) {
+ for(unsigned int i = 0; i < ge->getNumMeshElements(); ++i) {
MElement *t = ge->getMeshElement(i);
allElems.push_back(t);
- for (int k=0; k<t->getNumVertices(); k++)
+ for (int k = 0; k < t->getNumVertices(); k++)
dofView->numberVertex(t->getVertex(k), 0, 1);
}
}
}
-
- double L = norm(SVector3(bbox.max(), bbox.min()));
+
+ double L = norm(SVector3(bbox.max(), bbox.min()));
double mu = type*L;
-
+
simpleFunction<double> DIFF(mu*mu), ONE(1.0);
distanceTerm distance(GModel::current(), 1, &DIFF, &ONE);
-
- for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
+
+ for (std::vector<MElement* >::iterator it = allElems.begin();
+ it != allElems.end(); it++){
SElement se((*it));
distance.addToMatrix(*dofView, &se);
}
groupOfElements gr(allElems);
distance.addToRightHandSide(*dofView, gr);
-
+
Msg::Info("Distance Computation: Assembly done");
lsys->systemSolve();
Msg::Info("Distance Computation: System solved");
-
- for (std::map<MVertex*,double >::iterator itv = _distance_map.begin();
- itv != _distance_map.end() ; ++itv) {
+
+ for (std::map<MVertex*,double >::iterator itv = _distance_map.begin();
+ itv != _distance_map.end() ; ++itv) {
MVertex *v = itv->first;
double value;
dofView->getDofValue(v, 0, 1, value);
@@ -409,20 +416,19 @@ PView *GMSH_DistancePlugin::execute(PView *v)
double dist = -mu * log(1. - value);
itv->second = dist;
}
-
+
printView(_entities, _distance_map);
-
+
#endif
}
-
+
_data->setName("distance");
_data->Time.push_back(0);
_data->setFileName(_fileName.c_str());
_data->finalize();
-
- //compute also orthogonal vector to distance field
- // A Uortho = -C DIST
+ // compute also orthogonal vector to distance field
+ // A Uortho = -C DIST
//------------------------------------------------
if (ortho > 0) {
#if defined(HAVE_SOLVER)
@@ -437,9 +443,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
#endif
dofManager<double> myAssembler(lsys2);
simpleFunction<double> ONE(1.0);
-
+
double dMax = 1.0; //EMI TO CHANGE
-
+
std::vector<MElement *> allElems;
for(unsigned int ii = 0; ii < _entities.size(); ii++){
if (_entities[ii]->dim() == _maxDim) {
@@ -447,56 +453,59 @@ PView *GMSH_DistancePlugin::execute(PView *v)
for (unsigned int i=0; i<ge->getNumMeshElements(); ++i) {
MElement *t = ge->getMeshElement(i);
double vMean = 0.0;
- for (int k=0; k<t->getNumVertices(); k++) {
+ for (int k = 0; k < t->getNumVertices(); k++) {
std::map<MVertex*, double>::iterator it = _distance_map.find(t->getVertex(k));
vMean += it->second;
}
- vMean/=t->getNumVertices();
+ vMean /= t->getNumVertices();
if (vMean < dMax)
allElems.push_back(ge->getMeshElement(i));
}
}
- }
-
+ }
+
int mid = (int)floor(allElems.size() / 2.);
MElement *e = allElems[mid];
MVertex *vFIX = e->getVertex(0);
myAssembler.fixVertex(vFIX, 0, 1, 0.0);
-
- for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
+
+ for (std::vector<MElement* >::iterator it = allElems.begin();
+ it != allElems.end(); it++){
MElement *t = *it;
for(int k = 0; k < t->getNumVertices(); k++)
myAssembler.numberVertex(t->getVertex(k), 0, 1);
}
-
+
orthogonalTerm *ortho;
ortho = new orthogonalTerm(GModel::current(), 1, &ONE, &_distance_map);
// if (type < 0)
// ortho = new orthogonalTerm(GModel::current(), 1, &ONE, view);
// else
// ortho = new orthogonalTerm(GModel::current(), 1, &ONE, dofView);
-
- for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
+
+ for (std::vector<MElement* >::iterator it = allElems.begin();
+ it != allElems.end(); it++){
SElement se((*it));
ortho->addToMatrix(myAssembler, &se);
}
groupOfElements gr(allElems);
ortho->addToRightHandSide(myAssembler, gr);
-
+
Msg::Info("Orthogonal Computation: Assembly done");
lsys2->systemSolve();
Msg::Info("Orthogonal Computation: System solved");
-
+
PView *view2 = new PView();
PViewDataList *data2 = getDataList(view2);
data2->setName("ortogonal field");
-
+
Msg::Info("Writing orthogonal.pos");
FILE * f5 = fopen("orthogonal.pos","w");
fprintf(f5,"View \"orthogonal\"{\n");
- for (std::vector<MElement* >::iterator it = allElems.begin(); it != allElems.end(); it++){
+ for (std::vector<MElement* >::iterator it = allElems.begin();
+ it != allElems.end(); it++){
MElement *e = *it;
-
+
int numNodes = e->getNumVertices();
if (e->getType() == TYPE_POLYG)
numNodes = e->getNumChildren() * e->getChild(0)->getNumVertices();
@@ -504,7 +513,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
std::vector<double> *out2 = data2->incrementList(1, e->getType(), numNodes);
std::vector<MVertex*> nods;
std::vector<double> orth;
-
+
if(!e->getNumChildren())
for(int i=0; i<numNodes; i++)
nods.push_back(e->getVertex(i));
@@ -512,11 +521,11 @@ PView *GMSH_DistancePlugin::execute(PView *v)
for(int i = 0; i < e->getNumChildren(); i++)
for(int j = 0; j < e->getChild(i)->getNumVertices(); j++)
nods.push_back(e->getChild(i)->getVertex(j));
-
- for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->x());
- for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->y());
- for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->z());
-
+
+ for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->x());
+ for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->y());
+ for(int nod = 0; nod < numNodes; nod++) out2->push_back((nods[nod])->z());
+
if (_maxDim == 2)
switch (numNodes) {
case 2: fprintf(f5,"SL("); break;
@@ -532,7 +541,7 @@ PView *GMSH_DistancePlugin::execute(PView *v)
case 5: fprintf(f5,"SY("); break;
default: Msg::Fatal("Error in Plugin 'Distance' (numNodes=%g).",numNodes); break;
}
-
+
for (int j=0; j<numNodes; j++) {
MVertex *v = nods[j];
if (j)
@@ -550,9 +559,9 @@ PView *GMSH_DistancePlugin::execute(PView *v)
fprintf(f5,",%g", orth[i]);
else
fprintf(f5,"%g", orth[i]);
- }
+ }
fprintf(f5,"};\n");
-
+
}
fprintf(f5,"};\n");
fclose(f5);
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index 59b0e56a54a3517a6cb7b4dee8a8ab4d79616192..3935fd52356d70585ca9a220704471966cf3b87d 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -169,7 +169,7 @@ bool PViewDataGModel::finalize(bool computeMinMax, const std::string &interpolat
// for which we know the interpolation: it's constant)
int types[] = {TYPE_PNT, TYPE_LIN, TYPE_TRI, TYPE_QUA, TYPE_TET, TYPE_HEX,
TYPE_PRI, TYPE_PYR, TYPE_POLYG, TYPE_POLYH};
- for(int i = 0; i < sizeof(types) / sizeof(types[0]); i++){
+ for(unsigned int i = 0; i < sizeof(types) / sizeof(types[0]); i++){
if(!haveInterpolationMatrices(types[i])){
MElement *e = _getOneElementOfGivenType(model, types[i]);
if(e){
diff --git a/Post/PViewDataList.cpp b/Post/PViewDataList.cpp
index 5516811652de4050992e7414e04556bf3c9dd86b..35575f7369510447334ce69169796256d4aaacfe 100644
--- a/Post/PViewDataList.cpp
+++ b/Post/PViewDataList.cpp
@@ -19,7 +19,7 @@ PViewDataList::PViewDataList(bool isAdapted)
NbSS(0), NbVS(0), NbTS(0), NbSH(0), NbVH(0), NbTH(0),
NbSI(0), NbVI(0), NbTI(0), NbSY(0), NbVY(0), NbTY(0),
NbSD(0), NbVD(0), NbTD(0), NbT2(0), NbT3(0),
- _lastElement(-1), _lastDimension(-1), _lastNumNodes(-1),
+ _lastElement(-1), _lastDimension(-1), _lastNumNodes(-1),
_lastNumComponents(-1), _lastNumValues(-1), _lastNumEdges(-1),
_lastType(-1), _lastXYZ(0), _lastVal(0), _isAdapted(isAdapted)
{
@@ -101,18 +101,18 @@ bool PViewDataList::finalize(bool computeMinMax, const std::string &interpolatio
}
int PViewDataList::getNumScalars(int step)
-{
+{
return NbSP + NbSL + NbST + NbSQ + NbSS + NbSH + NbSI + NbSY + NbSG + NbSD;
}
int PViewDataList::getNumVectors(int step)
{
- return NbVP + NbVL + NbVT + NbVQ + NbVS + NbVH + NbVI + NbVY + NbVG + NbVD;
+ return NbVP + NbVL + NbVT + NbVQ + NbVS + NbVH + NbVI + NbVY + NbVG + NbVD;
}
int PViewDataList::getNumTensors(int step)
{
- return NbTP + NbTL + NbTT + NbTQ + NbTS + NbTH + NbTI + NbTY + NbTG + NbTD;
+ return NbTP + NbTL + NbTT + NbTQ + NbTS + NbTH + NbTI + NbTY + NbTG + NbTD;
}
int PViewDataList::getNumElements(int step, int ent)
@@ -215,7 +215,7 @@ void PViewDataList::_stat(std::vector<double> &list, int nbcomp, int nbelm,
}
int nb = list.size() / nbelm;
- for(unsigned int ele = 0; ele < nbelm; ele ++){
+ for(int ele = 0; ele < nbelm; ele ++){
int i = ele * nb;
if(type == TYPE_POLYG || type == TYPE_POLYH){
int t = (type == TYPE_POLYG) ? 0 : 1;
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index aaaf4ec24ee406f73d5a85ac82acd6a6519f54b1..83466b6c117b36882157552e7aa5424a5ecd7be8 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -2,6 +2,10 @@
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Emilie Marchandise
+//
#include "Context.h"
#include "multiscaleLaplace.h"
@@ -28,8 +32,6 @@
#include "MTriangle.h"
#include "robustPredicates.h"
-//--------------------------------------------------------------
-
struct compareRotatedPoints {
double angle;
const SPoint2 &left;
@@ -49,72 +51,23 @@ struct compareRotatedPoints {
}
};
-
struct sort_pred {
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) {
+ 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);
}
};
-static void sort_centers_dist(std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers,
- const SPoint2 &leftP){
-
- std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > myCenters(centers);
- std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > newCenters;
- SPoint2 lPoint = leftP;
-
- //printf("size centers =%d \n", myCenters.size());
- while (!myCenters.empty()){
- //printf("size centers loop =%d \n", myCenters.size());
- double dist = 1.e6;
- int numClosest = 0;
- for (int i= 0; i < myCenters.size(); i++){
- SPoint2 c = myCenters[i].first;
- double distc = sqrt((c.x() - lPoint.x())*(c.x() - lPoint.x())+
- (c.y() - lPoint.y())*(c.y() - lPoint.y()));
- //printf("distc =%g \n", distc);
- if (distc < dist){
- dist = distc;
- numClosest = i;
- }
- }
- //printf("numClosest =%d \n", numClosest);
- lPoint = myCenters[numClosest].first;
- newCenters.push_back(myCenters[numClosest]);
- myCenters.erase(myCenters.begin()+numClosest);
- }
-
- centers = newCenters;
-
-};
-
-
-//--------------------------------------------------------------
-static int intersection_segments_b (SPoint2 &p1, SPoint2 &p2,
- SPoint2 &q1, SPoint2 &q2,
- double x[2]){
- double P1[2] = {p1.x(),p1.y()};
- double P2[2] = {p2.x(),p2.y()};
- double Q1[2] = {q1.x(),q1.y()};
- double Q2[2] = {q2.x(),q2.y()};
-
- double PQ1 = robustPredicates::orient2d(P1,P2,Q1);
- double PQ2 = robustPredicates::orient2d(P1,P2,Q2);
-
- double QP1 = robustPredicates::orient2d(Q1,Q2,P1);
- 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;
+ 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);
@@ -126,10 +79,10 @@ static void recur_connect (MVertex *v,
}
}
-//--------------------------------------------------------------
-static void connected_bounds (std::vector<MElement*> &elements,
- std::vector<std::vector<MEdge> > &boundaries){
+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++){
@@ -160,9 +113,9 @@ static void connected_bounds (std::vector<MElement*> &elements,
return;
}
-//--------------------------------------------------------------
-static double getSizeBB(std::vector<MEdge> &me){
+static double getSizeBB(std::vector<MEdge> &me)
+{
SBoundingBox3d bb ;
SOrientedBoundingBox obbox ;
@@ -185,12 +138,11 @@ static double getSizeBB(std::vector<MEdge> &me){
double H = obbox.getMaxSize();
return H;
-
}
-//--------------------------------------------------------------
-static void ordering_dirichlet(std::vector<MElement*> &elements,
- std::vector<std::pair<MVertex*,double> > &dirichletNodes){
+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);
@@ -237,7 +189,8 @@ static void ordering_dirichlet(std::vector<MElement*> &elements,
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);
+ double iLength = dirichletNodes[dirichletNodes.size()-1].second +
+ (length / tot_length);
dirichletNodes.push_back(std::make_pair(v0,iLength));
break;
}
@@ -248,7 +201,8 @@ static void ordering_dirichlet(std::vector<MElement*> &elements,
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);
+ double iLength = dirichletNodes[dirichletNodes.size()-1].second +
+ (length / tot_length);
dirichletNodes.push_back(std::make_pair(v1,iLength));
break;
}
@@ -258,10 +212,11 @@ static void ordering_dirichlet(std::vector<MElement*> &elements,
return;
}
-//--------------------------------------------------------------
+
static int intersection_segments (SPoint2 &p1, SPoint2 &p2,
SPoint2 &q1, SPoint2 &q2,
- double x[2]){
+ double x[2])
+{
double A[2][2];
A[0][0] = p2.x()-p1.x();
A[0][1] = q1.x()-q2.x();
@@ -274,10 +229,10 @@ 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, int &nbElems ){
+static void recur_compute_centers_ (double R, double a1, double a2,
+ multiscaleLaplaceLevel * root, int &nbElems)
+{
nbElems += root->elements.size();
root->radius = R;
@@ -328,7 +283,8 @@ static void recur_compute_centers_ (double R, double a1, double a2,
//check if the center has not been added
bool newCenter = true;
- for (std::vector<SPoint2>::iterator it2 = centersChild.begin(); it2 != centersChild.end(); it2++){
+ for (std::vector<SPoint2>::iterator it2 = centersChild.begin();
+ it2 != centersChild.end(); it2++){
SPoint2 p = *it2;
double dist = sqrt ((c.x() - p.x())*(c.x() - p.x())+
(c.y() - p.y())*(c.y() - p.y()));
@@ -356,41 +312,42 @@ static void recur_compute_centers_ (double R, double a1, double a2,
centers.insert(centers.begin(), std::make_pair(PL,zero));
centers.push_back(std::make_pair(PR,zero));
- for (unsigned int i=1;i<centers.size()-1;i++){
- multiscaleLaplaceLevel* m1 = centers[i-1].second;
+ for (unsigned int i = 1; i < centers.size() - 1; i++){
multiscaleLaplaceLevel* m2 = centers[i].second;
- multiscaleLaplaceLevel* m3 = centers[i+1].second;
if (m2){
- a1 = myatan2 (centers[i-1].first.y()- m2->center.y() , centers[i-1].first.x()-m2->center.x());
- a2 = myatan2 (centers[i+1].first.y()- m2->center.y() , centers[i+1].first.x()-m2->center.x());
+ a1 = myatan2 (centers[i-1].first.y()- m2->center.y(),
+ centers[i-1].first.x()-m2->center.x());
+ a2 = myatan2 (centers[i+1].first.y()- m2->center.y(),
+ centers[i+1].first.x()-m2->center.x());
recur_compute_centers_ (m2->radius, a1, a2, m2, nbElems);
}
}
}
-//--------------------------------------------------------------
+
static void recur_cut_edges_ (multiscaleLaplaceLevel *root,
std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
- std::set<MVertex*> &cutVertices){
+ std::set<MVertex*> &cutVertices)
+{
const double EPS = 0.001;
std::multimap<MEdge,MElement*,Less_Edge> e2e;
std::set<MEdge,Less_Edge> allEdges;
- for (int i=0;i<root->elements.size();++i){
- for (int j=0;j<root->elements[i]->getNumEdges();j++){
+ for (unsigned int i = 0; i < root->elements.size(); ++i){
+ for (int j = 0; j < root->elements[i]->getNumEdges(); j++){
e2e.insert(std::make_pair(root->elements[i]->getEdge(j),root->elements[i]));
allEdges.insert(root->elements[i]->getEdge(j));
}
}
std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers = root->cut;
- for (unsigned int i=0;i< centers.size()-1;i++){
+ for (unsigned int i = 0; i < centers.size() - 1; i++){
SPoint2 p1 = centers[i].first;
SPoint2 p2 = centers[i+1].first;
//printf("*************** line p1p2 (%g %g) -- (%g %g) \n",p1.x(),p1.y(),p2.x(),p2.y());
for (std::set <MEdge,Less_Edge>::iterator it = allEdges.begin();
it != allEdges.end() ; ++it){
- if( cutEdges.find(*it) == cutEdges.end()){//e2e.count(*it) == 2 &&
+ if(cutEdges.find(*it) == cutEdges.end()){//e2e.count(*it) == 2 &&
std::map<MVertex *, SPoint2>::iterator it0 = root->coordinates.find(it->getVertex(0));
std::map<MVertex *, SPoint2>::iterator it1 = root->coordinates.find(it->getVertex(1));
if (it0 != root->coordinates.end() && it1 != root->coordinates.end()){
@@ -399,9 +356,10 @@ static void recur_cut_edges_ (multiscaleLaplaceLevel *root,
double x[2];
int inters = intersection_segments (p1,p2,q1,q2,x);
if (inters && x[1] > EPS && x[1] < 1.-EPS){
- MVertex *newv = new MVertex ((1.-x[1])*it->getVertex(0)->x() + x[1]*it->getVertex(1)->x(),
- (1.-x[1])*it->getVertex(0)->y() + x[1]*it->getVertex(1)->y(),
- (1.-x[1])*it->getVertex(0)->z() + x[1]*it->getVertex(1)->z());
+ MVertex *newv = new MVertex
+ ((1.-x[1])*it->getVertex(0)->x() + x[1]*it->getVertex(1)->x(),
+ (1.-x[1])*it->getVertex(0)->y() + x[1]*it->getVertex(1)->y(),
+ (1.-x[1])*it->getVertex(0)->z() + x[1]*it->getVertex(1)->z());
cutEdges[*it] = newv;
root->coordinates[newv] = q1*(1.-x[1]) + q2*x[1] ;
}
@@ -418,13 +376,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<MElement*> &_all)
+{
std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers = root->cut;
std::vector<MElement*> newElements;
for (unsigned int i = 0; i < root->elements.size(); i++){
@@ -438,25 +396,30 @@ static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
}
if (c[0] && c[1]){
newElements.push_back(new MTriangle (c[0],root->elements[i]->getVertex(1),c[1]));
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],root->elements[i]->getVertex(2)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],
+ root->elements[i]->getVertex(2)));
newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[0],c[1]));
theCut.insert(MEdge(c[0],c[1]));
}
else if (c[0] && c[2]){
newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],c[2]));
- newElements.push_back(new MTriangle (c[0],root->elements[i]->getVertex(1),root->elements[i]->getVertex(2)));
+ newElements.push_back(new MTriangle (c[0],root->elements[i]->getVertex(1),
+ root->elements[i]->getVertex(2)));
newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[2],c[0]));
theCut.insert(MEdge(c[0],c[2]));
}
else if (c[1] && c[2]){
newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[2],c[1]));
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),root->elements[i]->getVertex(1),c[2]));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),
+ root->elements[i]->getVertex(1),c[2]));
newElements.push_back(new MTriangle (c[2],root->elements[i]->getVertex(1),c[1]));
theCut.insert(MEdge(c[1],c[2]));
}
else if (c[0]){
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],root->elements[i]->getVertex(2)));
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[0],root->elements[i]->getVertex(1)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[0],
+ root->elements[i]->getVertex(2)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(2),c[0],
+ root->elements[i]->getVertex(1)));
if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end()){
theCut.insert(MEdge(c[0],root->elements[i]->getVertex(0)));
}
@@ -468,8 +431,10 @@ static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
}
}
else if (c[1]){
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),c[1],root->elements[i]->getVertex(0)));
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[1],root->elements[i]->getVertex(2)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),c[1],
+ root->elements[i]->getVertex(0)));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),c[1],
+ root->elements[i]->getVertex(2)));
if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end()){
theCut.insert(MEdge(c[1],root->elements[i]->getVertex(0)));
}
@@ -481,8 +446,10 @@ static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
}
}
else if (c[2]){
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),root->elements[i]->getVertex(1), c[2]));
- newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),root->elements[i]->getVertex(2), c[2]));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(0),
+ root->elements[i]->getVertex(1), c[2]));
+ newElements.push_back(new MTriangle (root->elements[i]->getVertex(1),
+ root->elements[i]->getVertex(2), c[2]));
if (cutVertices.find (root->elements[i]->getVertex(0)) != cutVertices.end()){
theCut.insert(MEdge(c[2],root->elements[i]->getVertex(0)));
}
@@ -518,12 +485,11 @@ static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
}
}
-//--------------------------------------------------------------
static void recur_leftCut_ (MElement *e,
std::multimap<MEdge,MElement*,Less_Edge> &e2e,
std::set<MEdge,Less_Edge> &theCut,
- std::set<MElement*> &leftSet){
-
+ std::set<MElement*> &leftSet)
+{
if (leftSet.find(e) != leftSet.end())return;
leftSet.insert(e);
//printf("insert in left %d \n", e->getNum());
@@ -538,14 +504,14 @@ static void recur_leftCut_ (MElement *e,
}
}
-//--------------------------------------------------------------
// starting form a list of elements, returns
// lists of lists that are all simply connected
static void recur_connect (const MEdge &e,
std::multimap<MEdge,MElement*,Less_Edge> &e2e,
std::set<MElement*> &group,
- std::set<MEdge,Less_Edge> &touched){
+ std::set<MEdge,Less_Edge> &touched)
+{
if (touched.find(e) != touched.end())return;
touched.insert(e);
for (std::multimap <MEdge,MElement*,Less_Edge>::iterator it = e2e.lower_bound(e);
@@ -557,7 +523,6 @@ static void recur_connect (const MEdge &e,
}
}
-//--------------------------------------------------------------
static void connectedRegions (const std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > ®ions)
{
@@ -578,9 +543,9 @@ static void connectedRegions (const std::vector<MElement*> &elements,
e2e.erase(*it);
}
}
-//--------------------------------------------------------------
-static void keepConnected(std::vector<MElement*> &goodSize, std::vector<MElement*> &tooSmall){
+static void keepConnected(std::vector<MElement*> &goodSize, std::vector<MElement*> &tooSmall)
+{
std::vector<std::vector<MElement*> > regGoodSize;
connectedRegions (goodSize,regGoodSize);
if (regGoodSize.size() > 0){
@@ -594,18 +559,22 @@ static void keepConnected(std::vector<MElement*> &goodSize, std::vector<MElement
}
}
goodSize.clear();
- for (unsigned int i=0;i< regGoodSize.size() ; i++){
- if (i == index) goodSize.insert(goodSize.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
- else tooSmall.insert(tooSmall.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
+ for (unsigned int i = 0; i < regGoodSize.size(); i++){
+ if ((int)i == index)
+ goodSize.insert(goodSize.begin(), regGoodSize[i].begin(),
+ regGoodSize[i].end());
+ else
+ tooSmall.insert(tooSmall.begin(), regGoodSize[i].begin(),
+ regGoodSize[i].end());
}
}
}
-//--------------------------------------------------------------
+
static void recur_cut_ (double R, double a1, double a2,
multiscaleLaplaceLevel * root,
std::vector<MElement *> &left,
- std::vector<MElement *> &right){
-
+ 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;
@@ -633,21 +602,20 @@ static void recur_cut_ (double R, double a1, double a2,
}
for (unsigned 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){
- a1 = myatan2 (centers[i-1].first.y() - m2->center.y() , centers[i-1].first.x() - m2->center.x() );
- a2 = myatan2 (centers[i+1].first.y() - m2->center.y() , centers[i+1].first.x() - m2->center.x() );
+ a1 = myatan2 (centers[i-1].first.y() - m2->center.y(),
+ centers[i-1].first.x() - m2->center.x());
+ a2 = myatan2 (centers[i+1].first.y() - m2->center.y(),
+ centers[i+1].first.x() - m2->center.x());
recur_cut_ (m2->radius, a1, a2, m2, left, right);
}
}
}
-//--------------------------------------------------------------
static void connected_left_right (std::vector<MElement *> &left,
- std::vector<MElement *> &right ){
-
+ std::vector<MElement *> &right)
+{
//connected left
keepConnected(left, right);
@@ -656,11 +624,11 @@ static void connected_left_right (std::vector<MElement *> &left,
left[i]->setPartition(1);
for (unsigned int i= 0; i< right.size(); i++)
right[i]->setPartition(2);
-
}
-//--------------------------------------------------------------
-static void printCut(std::map<MEdge,MVertex*,Less_Edge> &cutEdges, std::set<MEdge,Less_Edge> &theCut, std::set<MVertex*> cutVertices){
+static void printCut(std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
+ std::set<MEdge,Less_Edge> &theCut, std::set<MVertex*> cutVertices)
+{
printf("Writing points.pos \n");
std::map<MEdge,MVertex*,Less_Edge>::iterator ite = cutEdges.begin();
FILE *f1 = fopen("points.pos","w");
@@ -680,13 +648,14 @@ static void printCut(std::map<MEdge,MVertex*,Less_Edge> &cutEdges, std::set<MEd
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(),
+ 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);
}
-//--------------------------------------------------------------
+
static void printLevel(const char* fn,
std::vector<MElement *> &elements,
std::map<MVertex*,SPoint2> *coordinates,
@@ -713,7 +682,8 @@ static void printLevel(const char* fn,
(*it)->setIndex(index++);
SPoint2 p = (coordinates) ? (*coordinates)[*it] : SPoint2(0,0);
if (coordinates) {
- if (dx)fprintf(fp, "%d %22.15E %22.15E 0\n", (*it)->getIndex(), dx[2]*(p.x()-dx[0]), dx[2]*(p.y()-dx[1]));
+ if (dx)fprintf(fp, "%d %22.15E %22.15E 0\n", (*it)->getIndex(),
+ dx[2]*(p.x()-dx[0]), dx[2]*(p.y()-dx[1]));
else fprintf(fp, "%d %22.15E %22.15E 0\n", (*it)->getIndex(), p.x(), p.y());
}
else fprintf(fp, "%d %g %g %g\n", (*it)->getIndex(),(*it)->x(), (*it)->y(), (*it)->z());
@@ -728,13 +698,9 @@ static void printLevel(const char* fn,
fclose(fp);
}
-//--------------------------------------------------------------
-
-
-
-
-static double localSize(MElement *e, std::map<MVertex*,SPoint2> &solution){
+static double localSize(MElement *e, std::map<MVertex*,SPoint2> &solution)
+{
SBoundingBox3d local;
for(int j = 0; j<e->getNumVertices(); ++j){
SPoint2 p = solution[e->getVertex(j)];
@@ -758,18 +724,14 @@ static double localSize(MElement *e, std::map<MVertex*,SPoint2> &solution){
// double a_2D = fabs(triangle_area(q0, q1, q2));
// return a_2D; //a_2D / a_3D;
-
-
}
-
-//--------------------------------------------------------------
static void printLevel_onlysmall(const char* fn,
std::vector<MElement *> &elements,
std::map<MVertex*,SPoint2> *coordinates,
double version,
- double tolerance){
-
+ double tolerance)
+{
std::vector<MElement *> small;
double dx[3] = {0,0,0};
int COUNT = 0;
@@ -791,55 +753,9 @@ static void printLevel_onlysmall(const char* fn,
printLevel(fn,small,coordinates,version,dx);
}
-//-------------------------------------------------------------
-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;
-
-}
-
-//--------------------------------------------------------------
-
multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
std::map<MVertex*, SPoint3> &allCoordinates)
{
-
//To go through this execute gmsh with the option -optimize_hom
//if (!CTX::instance()->mesh.smoothInternalEdges)return;
@@ -878,14 +794,13 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
//or Cut the mesh in left and right
splitElems(elements);
//cutElems(elements);
-
}
void multiscaleLaplace::fillCoordinates (multiscaleLaplaceLevel & level,
std::map<MVertex*, SPoint3> &allCoordinates,
std::vector<double> &iScale,
- std::vector<SPoint2> &iCenter){
-
+ std::vector<SPoint2> &iCenter)
+{
iScale.push_back(level.scale);
iCenter.push_back(level.center);
@@ -901,17 +816,14 @@ void multiscaleLaplace::fillCoordinates (multiscaleLaplaceLevel & level,
}
}
-
for (unsigned int i=0;i<level.children.size();i++){
multiscaleLaplaceLevel* m = level.children[i];
fillCoordinates(*m, allCoordinates, iScale, iCenter);
- }
-
-
+ }
}
-void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
-
+void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level)
+{
//Compute all nodes for the level
std::set<MVertex*> allNodes;
for(unsigned int i = 0; i < level.elements.size(); ++i){
@@ -940,7 +852,7 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
MElement *e = level.elements[i];
std::vector<SPoint2> localCoord;
double local_size = localSize(e,solution);
- if (local_size < 1.e-6*global_size)
+ if (local_size < 1.e-6*global_size)
tooSmall.push_back(e);
else goodSize.push_back(e);
}
@@ -968,18 +880,18 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
}
//check for convex small regions patches
- for (int i=0;i< regions.size() ; i++){
+ for (unsigned int i = 0; i < regions.size(); i++){
std::vector<MElement*> &elemR = regions[i];
v2t_cont adj;
buildVertexToElement (elemR,adj);
for (std::vector<MElement*>::iterator it = elemR.begin(); it != elemR.end(); ++it){
int nbNeigh = 0;
MElement *e = *it;
- v2t_cont :: iterator it0 = adj.find(e->getVertex(0));
+ v2t_cont::iterator it0 = adj.find(e->getVertex(0));
if(it0 != adj.end()) nbNeigh += it0->second.size();
- v2t_cont :: iterator it1 = adj.find(e->getVertex(1));
+ v2t_cont::iterator it1 = adj.find(e->getVertex(1));
if(it1 != adj.end()) nbNeigh += it1->second.size();
- v2t_cont :: iterator it2 = adj.find(e->getVertex(2));
+ v2t_cont::iterator it2 = adj.find(e->getVertex(2));
if(it2 != adj.end()) nbNeigh += it2->second.size();
std::vector<MElement*>::iterator itp;
if (nbNeigh < 12) {
@@ -992,8 +904,8 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
keepConnected(elemR, goodSize);
}
tooSmall.clear();
- for (int i=0;i< regions.size() ; i++)
- tooSmall.insert(tooSmall.begin(), regions[i].begin(), regions[i].end());
+ for (unsigned int i = 0; i < regions.size(); i++)
+ tooSmall.insert(tooSmall.begin(), regions[i].begin(), regions[i].end());
keepConnected(goodSize, tooSmall);
regions.clear();
@@ -1022,7 +934,8 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
printLevel_onlysmall (name3.c_str(),level.elements,&level.coordinates,2.2,1.e-15);
//For every small region compute a new parametrization
- Msg::Info("Level (%d-%d): %d connected small regions",level.recur, level.region, regions.size());
+ Msg::Info("Level (%d-%d): %d connected small regions",
+ level.recur, level.region, regions.size());
for (unsigned int i = 0; i < regions.size(); i++){
std::set<MVertex*> tooSmallv;
tooSmallv.clear();
@@ -1041,7 +954,8 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
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){
+ for(std::set<MVertex *>::iterator itv = tooSmallv.begin();
+ itv !=tooSmallv.end(); ++itv){
SPoint2 p = solution[*itv];
nextLevel->center += p;
smallB += SPoint3(p.x(),p.y(),0.0);
@@ -1049,27 +963,28 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
nextLevel->center *= (1./(double)tooSmallv.size());
nextLevel->scale = smallB.max().distance(smallB.min());
- for(std::set<MVertex *>::iterator itv = tooSmallv.begin(); itv !=tooSmallv.end() ; ++itv){
+ for(std::set<MVertex *>::iterator itv = tooSmallv.begin();
+ itv !=tooSmallv.end() ; ++itv){
MVertex *v = *itv;
if (goodSizev.find(v) != goodSizev.end()){
- nextLevel->coordinates[v] = (solution[v]-nextLevel->center)*(1./nextLevel->scale);
+ nextLevel->coordinates[v] = (solution[v]-nextLevel->center) *
+ (1./nextLevel->scale);
}
}
// recursively continue if tooSmall is not empty
if (!tooSmallv.empty()){
- Msg::Info("Level (%d-%d) Multiscale Laplace (reg[%d] = %d too small)",level.recur,level.region, i, tooSmallv.size());
+ Msg::Info("Level (%d-%d) Multiscale Laplace (reg[%d] = %d too small)",
+ level.recur,level.region, i, tooSmallv.size());
level.children.push_back(nextLevel);
parametrize (*nextLevel);
}
}
-
}
void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
std::set<MVertex*> &allNodes,
std::map<MVertex*,SPoint2> &solution)
{
-
linearSystem<double> *_lsys;
#if defined(HAVE_PETSC)
_lsys = new linearSystemPETSc<double>;
@@ -1094,7 +1009,8 @@ void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
}
// do the numbering
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ for(std::set<MVertex *>::iterator itv = allNodes.begin();
+ itv != allNodes.end(); ++itv){
MVertex *v = *itv;
myAssembler.numberVertex(v, 0, 1);
}
@@ -1113,8 +1029,8 @@ void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
if (myAssembler.sizeOfR() != 0) _lsys->systemSolve();
// get the values
- int count = 0;
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ for(std::set<MVertex *>::iterator itv = allNodes.begin();
+ itv != allNodes.end(); ++itv){
MVertex *v = *itv;
double value;
myAssembler.getDofValue(v, 0, 1, value);
@@ -1130,7 +1046,6 @@ void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
void multiscaleLaplace::cutElems(std::vector<MElement *> &elements)
{
-
std::map<MEdge,MVertex*,Less_Edge> cutEdges;
std::set<MEdge,Less_Edge> theCut;
std::set<MVertex*> cutVertices;
@@ -1141,8 +1056,8 @@ void multiscaleLaplace::cutElems(std::vector<MElement *> &elements)
printCut(cutEdges, theCut, cutVertices);
std::multimap<MEdge,MElement*,Less_Edge> e2e;
- for (int i=0;i<elements.size();++i){
- for (int j=0;j<elements[i]->getNumEdges();j++){
+ for (unsigned 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]));
}
}
@@ -1151,15 +1066,16 @@ void multiscaleLaplace::cutElems(std::vector<MElement *> &elements)
std::vector<MElement*> left,right;
recur_leftCut_ (elements[0], e2e, theCut, leftS);
- for (int i=0;i< elements.size();i++){
+ for (unsigned int i = 0; i < elements.size(); i++){
MElement *e = elements[i];
if (leftS.find(e) != leftS.end()) left.push_back(e);
else right.push_back(e);
}
connected_left_right(left, right);
- if (left.size()== 0 || right.size() == 0) {
- printf("KO size left=%d, right=%d not good (zero elems)\n", (int) left.size(), (int) right.size() );
+ if (left.size() == 0 || right.size() == 0) {
+ printf("KO size left=%d, right=%d not good (zero elems)\n",
+ (int) left.size(), (int) right.size() );
exit(1);
}
@@ -1172,11 +1088,10 @@ void multiscaleLaplace::cutElems(std::vector<MElement *> &elements)
printLevel ("Rootcut-all.msh",elements, 0,2.2);
//exit(1);
}
+
void multiscaleLaplace::splitElems(std::vector<MElement *> &elements)
{
-
std::vector<MElement*> left,right;
- int totNbElems = 0;
recur_cut_ (1.0, M_PI, 0.0, root,left,right);
connected_left_right(left, right);
@@ -1185,7 +1100,7 @@ void multiscaleLaplace::splitElems(std::vector<MElement *> &elements)
printLevel ("Rootsplit-all.msh",elements, 0,2.2);
printLevel ("Rootsplit-left-param.msh",left,&root->coordinates,2.2);
- //printLevel_onlysmall ("Rootsplit-left-param10.msh",left,&root->coordinates,2.2,1.e-10);
+ // printLevel_onlysmall ("Rootsplit-left-param10.msh",left,&root->coordinates,2.2,1.e-10);
// printLevel_onlysmall ("Rootsplit-left-param12.msh",left,&root->coordinates,2.2,1.e-12);
// printLevel_onlysmall ("Rootsplit-left-param15.msh",left,&root->coordinates,2.2,1.e-15);
@@ -1198,13 +1113,12 @@ void multiscaleLaplace::splitElems(std::vector<MElement *> &elements)
// printLevel_onlysmall ("Rootsplit-all-param15.msh",elements,&root->coordinates,2.2,1.e-15);
if ( elements.size() != left.size()+right.size()) {
- Msg::Error("Cutting laplace wrong nb elements (%d) != left + right (%d)", elements.size(), left.size()+right.size());
+ Msg::Error("Cutting laplace wrong nb elements (%d) != left + right (%d)",
+ elements.size(), left.size()+right.size());
exit(1);
}
elements.clear();
elements.insert(elements.end(),left.begin(),left.end());
elements.insert(elements.end(),right.begin(),right.end());
-
}
-
diff --git a/Solver/orthogonalTerm.h b/Solver/orthogonalTerm.h
index fb8dedb2e899ccfd548a59ead830f7883806068e..c82f730b9286538a1b59d35e968d41ec0b50d9a4 100644
--- a/Solver/orthogonalTerm.h
+++ b/Solver/orthogonalTerm.h
@@ -16,11 +16,14 @@ class orthogonalTerm : public helmholtzTerm<double> {
bool withDof;
std::map<MVertex*, double > *_distance_map;
public:
- orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k, std::map<MVertex*, double > *distance_map)
- : helmholtzTerm<double>(gm, iField, iField, k, 0), _distance_map(distance_map), withDof(false) {}
- orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k, PView *pview)
- : helmholtzTerm<double>(gm, iField, iField, k, 0), _pview(pview), withDof(false) {}
- orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k, dofManager<double> *dofView)
+ orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k,
+ std::map<MVertex*, double > *distance_map)
+ : helmholtzTerm<double>(gm, iField, iField, k, 0), withDof(false),
+ _distance_map(distance_map) {}
+ orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k, PView *pview)
+ : helmholtzTerm<double>(gm, iField, iField, k, 0), _pview(pview), withDof(false) {}
+ orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k,
+ dofManager<double> *dofView)
: helmholtzTerm<double>(gm, iField, iField, k, 0), _dofView(dofView), withDof(true) {}
void elementVector(SElement *se, fullVector<double> &m) const
{
@@ -38,20 +41,20 @@ class orthogonalTerm : public helmholtzTerm<double> {
e->getIntegrationPoints(integrationOrder, &npts, &GP);
fullMatrix<double> mat(nbSF, nbSF);
mat.setAll(0.);
-
+
for(int i = 0; i < npts; i++){
const double u = GP[i].pt[0];
const double v = GP[i].pt[1];
const double w = GP[i].pt[2];
const double weight = GP[i].weight;
- const double detJ = e->getJacobian(u, v, w, jac);
+ const double detJ = e->getJacobian(u, v, w, jac);
SPoint3 p; e->pnt(u, v, w, p);
- inv3x3(jac, invjac);
+ inv3x3(jac, invjac);
e->getGradShapeFunctions(u, v, w, grads);
for(int j = 0; j < nbSF; j++){
- Grads[j] = SVector3(invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
+ Grads[j] = SVector3(invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
invjac[0][2] * grads[j][2],
- invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
+ invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
invjac[1][2] * grads[j][2],
invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
invjac[2][2] * grads[j][2]);
@@ -59,12 +62,12 @@ class orthogonalTerm : public helmholtzTerm<double> {
SVector3 N (jac[2][0], jac[2][1], jac[2][2]);
for(int j = 0; j < nbSF; j++)
for(int k = 0; k <= j; k++)
- mat(j, k) += dot(crossprod(Grads[j], Grads[k]), N) * weight * detJ;
+ mat(j, k) += dot(crossprod(Grads[j], Grads[k]), N) * weight * detJ;
}
for(int j = 0; j < nbSF; j++)
for(int k = 0; k < j; k++)
mat(k, j) = -1.* mat(j, k);
-
+
//2) compute vector m(i) = mat(i,j)*val(j)
fullVector<double> val(nbSF);
val.scale(0.);
@@ -89,7 +92,7 @@ class orthogonalTerm : public helmholtzTerm<double> {
/* /\* } *\/ */
/* } */
- m.scale(0.);
+ m.scale(0.);
for(int i = 0; i < nbSF; i++)
for(int j = 0; j < nbSF; j++)
m(i) += -mat(i, j) * val(j);