diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index 0bbe5ee232c39c2a1764f711391a77b8d92eb117..9793c2eca60372bbde04d97a3d43551e9e7c3cf6 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
// compute the largest inscribed ellipsoid...
#include "STensor3.h"
@@ -73,8 +78,8 @@ SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2
}
// (1-t) * m1 + t * m2
-SMetric3 interpolation (const SMetric3 &m1,
- const SMetric3 &m2,
+SMetric3 interpolation (const SMetric3 &m1,
+ const SMetric3 &m2,
const double t)
{
SMetric3 im1 = m1.invert();
@@ -85,9 +90,9 @@ SMetric3 interpolation (const SMetric3 &m1,
return im1.invert();
}
-SMetric3 interpolation (const SMetric3 &m1,
- const SMetric3 &m2,
- const SMetric3 &m3,
+SMetric3 interpolation (const SMetric3 &m1,
+ const SMetric3 &m2,
+ const SMetric3 &m3,
const double u,
const double v)
{
@@ -102,10 +107,10 @@ SMetric3 interpolation (const SMetric3 &m1,
return im1.invert();
}
-SMetric3 interpolation (const SMetric3 &m1,
- const SMetric3 &m2,
- const SMetric3 &m3,
- const SMetric3 &m4,
+SMetric3 interpolation (const SMetric3 &m1,
+ const SMetric3 &m2,
+ const SMetric3 &m3,
+ const SMetric3 &m4,
const double u,
const double v,
const double w)
diff --git a/Mesh/meshGFaceBoundaryLayers.cpp b/Mesh/meshGFaceBoundaryLayers.cpp
index 9dec85a730c019b783e756454f9efcfe85210dbf..ce42e25918a343af815fa9f0603e9ee93a198292 100644
--- a/Mesh/meshGFaceBoundaryLayers.cpp
+++ b/Mesh/meshGFaceBoundaryLayers.cpp
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include "GModel.h"
#include "GFace.h"
#include "MVertex.h"
@@ -18,7 +23,7 @@
SVector3 interiorNormal (SPoint2 p1, SPoint2 p2, SPoint2 p3){
-
+
SVector3 ez (0,0,1);
SVector3 d (p1.x()-p2.x(),p1.y()-p2.y(),0);
SVector3 h (p3.x()-0.5*(p2.x()+p1.x()),p3.y()-0.5*(p2.y()+p1.y()),0);
@@ -27,7 +32,7 @@ SVector3 interiorNormal (SPoint2 p1, SPoint2 p2, SPoint2 p3){
// printf("%g %g\n",n.x(),n.y());
if (dot(n,h) > 0)return n;
return n*(-1.);
-
+
}
double computeAngle (GFace *gf, const MEdge &e1, const MEdge &e2, SVector3 &n1, SVector3 &n2){
@@ -40,22 +45,22 @@ double computeAngle (GFace *gf, const MEdge &e1, const MEdge &e2, SVector3 &n1,
MVertex *v22 = e2.getVertex(1);
MVertex *v0,*v1,*v2;
if (v11 == v21){
- v0 = v12 ; v1 = v11 ; v2 = v22;
+ v0 = v12 ; v1 = v11 ; v2 = v22;
}
else if (v11 == v22){
- v0 = v12 ; v1 = v11 ; v2 = v21;
+ v0 = v12 ; v1 = v11 ; v2 = v21;
}
else if (v12 == v21){
- v0 = v11 ; v1 = v12 ; v2 = v22;
+ v0 = v11 ; v1 = v12 ; v2 = v22;
}
else if (v12 == v22){
- v0 = v11 ; v1 = v12 ; v2 = v21;
+ v0 = v11 ; v1 = v12 ; v2 = v21;
}
else throw;
-
+
reparamMeshEdgeOnFace(v0, v1, gf, p0, p1);
reparamMeshEdgeOnFace(v0, v2, gf, p0, p2);
-
+
SVector3 t1 (p1.x()-p0.x(),p1.y()-p0.y(),0);
SVector3 t2 (p2.x()-p1.x(),p2.y()-p1.y(),0);
t1.normalize();
@@ -63,10 +68,10 @@ double computeAngle (GFace *gf, const MEdge &e1, const MEdge &e2, SVector3 &n1,
SVector3 n = crossprod(t1,t2);
double sign = dot(t1,n2);
-
+
double a = atan2 (n.z(),cosa);
a = sign > 0 ? fabs(a) : -fabs(a);
-
+
// printf("a = %12.5e cos %12.5E sin %12.5E %g %g vs %g %g\n",a,cosa,n.z(),n1.x(),n1.y(),n2.x(),n2.y());
return a;
}
@@ -76,17 +81,17 @@ void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3])
Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(uv[0], uv[1]));
double res[2][2];
-
+
double M[2][3] = {{der.first().x(),der.first().y(),der.first().z()},
{der.second().x(),der.second().y(),der.second().z()}};
-
-
+
+
for (int i=0;i<3;i++){
for (int l=0;l<3;l++){
res[i][l] = 0;
for (int j=0;j<3;j++){
for (int k=0;k<3;k++){
- res[i][l] += M[i][j]*m(j,k)*M[l][k];
+ res[i][l] += M[i][j]*m(j,k)*M[l][k];
}
}
}
@@ -94,7 +99,7 @@ void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3])
metric[0] = res[0][0];
metric[1] = res[1][0];
metric[2] = res[1][1];
-}
+}
BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
{
@@ -108,7 +113,7 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
}
Field *bl_field = fields->get(fields->getBoundaryLayerField());
BoundaryLayerField *blf = dynamic_cast<BoundaryLayerField*> (bl_field);
-
+
if (!blf)return 0;
double _treshold = blf->fan_angle * M_PI / 180 ;
@@ -116,9 +121,9 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
BoundaryLayerColumns * _columns = new BoundaryLayerColumns;
// assume a field that i) gives us the closest point on one of the BL components
- // ii) Gives us mesh sizes in the 3 directions.
+ // ii) Gives us mesh sizes in the 3 directions.
- // build vertex to vertex connexions
+ // build vertex to vertex connexions
std::list<GEdge*> edges = gf->edges();
std::list<GEdge*> embedded_edges = gf->embeddedEdges();
edges.insert(edges.begin(), embedded_edges.begin(),embedded_edges.end());
@@ -137,10 +142,10 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
}
// printf("%d boundary points\n",_vertices.size());
-
+
// assume that the initial mesh has been created i.e. that there exist
// triangles inside the domain. Triangles are used to define
- // exterior normals
+ // exterior normals
for (int i=0;i<gf->triangles.size();i++){
SPoint2 p0,p1,p2;
MVertex *v0 = gf->triangles[i]->getVertex(0);
@@ -150,21 +155,21 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
reparamMeshEdgeOnFace(v0, v2, gf, p0, p2);
MEdge me01(v0,v1);
- SVector3 v01 = interiorNormal (p0,p1,p2);
+ SVector3 v01 = interiorNormal (p0,p1,p2);
_columns->_normals.insert(std::make_pair(me01,v01));
MEdge me02(v0,v2);
- SVector3 v02 = interiorNormal (p0,p2,p1);
+ SVector3 v02 = interiorNormal (p0,p2,p1);
_columns->_normals.insert(std::make_pair(me02,v02));
MEdge me21(v2,v1);
- SVector3 v21 = interiorNormal (p2,p1,p0);
+ SVector3 v21 = interiorNormal (p2,p1,p0);
_columns->_normals.insert(std::make_pair(me21,v21));
- }
+ }
// for all boundry points
for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
- std::vector<MVertex*> _connections;
+ std::vector<MVertex*> _connections;
std::vector<SVector3> _dirs;
double LL;
for ( std::multimap<MVertex*,MVertex*> :: iterator itm = _columns->_non_manifold_edges.lower_bound(*it);
@@ -224,7 +229,7 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
if (_connections.size() == 2){
MEdge e1 (*it,_connections[0]);
MEdge e2 (*it,_connections[1]);
- LL = 0.5 * (e1.length() + e2.length());
+ LL = 0.5 * (e1.length() + e2.length());
std::vector<SVector3> N1,N2;
for ( std::multimap<MEdge,SVector3,Less_Edge> :: iterator itm = _columns->_normals.lower_bound(e1);
itm != _columns->_normals.upper_bound(e1); ++itm) N1.push_back(itm->second);
@@ -233,16 +238,16 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
double LL;
if (N1.size() == N2.size()){
// if (N1.size() > 1)printf("%d sides\n",N1.size());
- for (int SIDE = 0; SIDE < N1.size() ; SIDE++){
+ for (int SIDE = 0; SIDE < N1.size() ; SIDE++){
// IF THE ANGLE IS GREATER THAN THRESHOLD, ADD DIRECTIONS !!
- double angle = computeAngle (gf,e1,e2,N1[SIDE],N2[SIDE]);
+ double angle = computeAngle (gf,e1,e2,N1[SIDE],N2[SIDE]);
// if (N1.size() > 1)printf("angle = %g\n",angle);
if (angle < _treshold /*&& angle > - _treshold*/){
SVector3 x = N1[SIDE]*1.01+N2[SIDE];
x.normalize();
_dirs.push_back(x);
}
- else if (angle >= _treshold){
+ else if (angle >= _treshold){
int fanSize = angle / _treshold;
// printf("ONE FAN HAS BEEN CREATED : %d %d %d %d ANGLE = %g | %g %g %g %g\n",e1.getVertex(0)->getNum(),
// e1.getVertex(1)->getNum(),e2.getVertex(0)->getNum(),e2.getVertex(1)->getNum(),
@@ -275,40 +280,40 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
ee[0] = ee[1];ee[1] = eee0;
}
_columns->addFan (*it,ee[0],ee[1],true);
-
+
for (int i=-1; i<=fanSize; i++){
double t = (double)(i+1)/ (fanSize+1);
double alpha = t * AMAX + (1.-t)* AMIN;
- // printf("%d %g %g %g %g\n",i,alpha,alpha1,alpha2,alpha2-alpha1);
+ // printf("%d %g %g %g %g\n",i,alpha,alpha1,alpha2,alpha2-alpha1);
SVector3 x (cos(alpha),sin(alpha),0);
x.normalize();
_dirs.push_back(x);
- }
+ }
}
}
}
}
-
+
// if (_dirs.size() > 1)printf("%d directions\n",_dirs.size());
- // now create the BL points
+ // now create the BL points
for (int DIR=0;DIR<_dirs.size();DIR++){
SPoint2 p;
- SVector3 n = _dirs[DIR];
-
+ SVector3 n = _dirs[DIR];
+
// < ------------------------------- > //
// N = X(p0+ e n) - X(p0) //
// = e * (dX/du n_u + dX/dv n_v) //
// < ------------------------------- > //
-
+
MVertex *current = *it;
reparamMeshVertexOnFace(current,gf,p);
-
+
int nbCol = 100;
std::vector<MVertex*> _column;
std::vector<SMetric3> _metrics;
- // printf("start with point %g %g (%g %g)\n",current->x(),current->y(),p.x(),p.y());
+ // printf("start with point %g %g (%g %g)\n",current->x(),current->y(),p.x(),p.y());
AttractorField *catt = 0;
SPoint3 _close;
double _current_distance;
@@ -317,17 +322,17 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
SMetric3 m;
double metric[3];
double l;
- (*bl_field)(current->x(),current->y(), current->z(), m, current->onWhat());
+ (*bl_field)(current->x(),current->y(), current->z(), m, current->onWhat());
if (!catt){
catt = blf->current_closest;
_close = blf->_closest_point;
_current_distance = blf -> current_distance;
}
- SPoint2 poffset (p.x() + 1.e-8 * n.x(),
+ SPoint2 poffset (p.x() + 1.e-8 * n.x(),
p.y() + 1.e-8 * n.y());
- buildMeshMetric(gf, poffset, m, metric);
+ buildMeshMetric(gf, poffset, m, metric);
const double l2 = n.x()*n.x()*metric[0] + 2*n.x()*n.y()*metric[1] + n.y()*n.y()*metric[2] ;
- l = 1./sqrt(l2);
+ l = 1./sqrt(l2);
// if (_dirs.size() > 1) printf("l = %g metric = %g %g %g dim %d tag %d \n",l,metric[0],metric[1],metric[2],current->onWhat()->dim(),current->onWhat()->tag());
// printf("%g %g\n",l,LL);
if (l >= blf->hfar){
@@ -354,7 +359,7 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
catt = blf->current_closest;
_close = blf->_closest_point;
_current_distance = blf -> current_distance;
- SPoint2 pnew (p.x() + l * n.x(),
+ SPoint2 pnew (p.x() + l * n.x(),
p.y() + l * n.y());
GPoint gp = gf->point (pnew);
MFaceVertex *_current = new MFaceVertex (gp.x(),gp.y(),gp.z(),gf,pnew.x(),pnew.y());
@@ -363,17 +368,17 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
_column.push_back(current);
// printf("pnew %g %g new point %g %g n %g %g\n",pnew.x(),pnew.y(),gp.x(),gp.y(),n.x(),n.y());
_metrics.push_back(m);
- // const double l = n[0]*m(0,0) +;
+ // const double l = n[0]*m(0,0) +;
if (_column.size() > nbCol)break; // FIXME
p = pnew;
}
_columns->addColumn(n,*it, _column, _metrics);
}
- }
- // HERE WE SHOULD FILTER THE POINTS IN ORDER NOT TO HAVE POINTS THAT ARE TOO CLOSE
+ }
+ // HERE WE SHOULD FILTER THE POINTS IN ORDER NOT TO HAVE POINTS THAT ARE TOO CLOSE
// DEBUG STUFF
-
+
FILE *f = fopen ("test.pos","w");
fprintf(f,"View \"\" {\n");
for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end() ; ++it){
@@ -385,10 +390,10 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
fprintf(f,"SP(%g,%g,%g){%d};\n",blv->x(),blv->y(),blv->z(),v->getNum());
}
}
- }
+ }
fprintf(f,"};\n");
fclose (f);
-
+
// END OF DEBUG STUFF
return _columns;
diff --git a/Mesh/meshGFaceBoundaryLayers.h b/Mesh/meshGFaceBoundaryLayers.h
index 9aae125eaa96bd7b82ba4f3d05f2f3f72f167760..385fb73c5b8a67a29eb2876f210dd0e1da4b8cb5 100644
--- a/Mesh/meshGFaceBoundaryLayers.h
+++ b/Mesh/meshGFaceBoundaryLayers.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef _MESHGFACE_BNDRYLR_
#define _MESHGFACE_BNDRYLR_
@@ -6,27 +11,29 @@
#include "MEdge.h"
#include <map>
#include <set>
+
class Field;
class GFace;
+
struct BoundaryLayerData
{
- BoundaryLayerData (const SVector3 & dir,
+ BoundaryLayerData (const SVector3 & dir,
std::vector<MVertex*> column,
std::vector<SMetric3> metrics)
: _n(dir), _column(column),_metrics(metrics)
{
-
- }
+
+ }
SVector3 _n;
std::vector<MVertex*> _column;
- std::vector<SMetric3> _metrics;
+ std::vector<SMetric3> _metrics;
};
-struct BoundaryLayerFan
+struct BoundaryLayerFan
{
MEdge _e1, _e2;
bool sense;
- BoundaryLayerFan (MEdge e1, MEdge e2 , bool s = true)
+ BoundaryLayerFan (MEdge e1, MEdge e2 , bool s = true)
: _e1(e1),_e2(e2) , sense (s)
{}
};
@@ -36,7 +43,7 @@ struct edgeColumn {
edgeColumn (const BoundaryLayerData & c1, const BoundaryLayerData & c2) : _c1(c1),_c2(c2){}
};
-class BoundaryLayerColumns
+class BoundaryLayerColumns
{
std::multimap<MVertex*,BoundaryLayerData> _data;
std::map<MVertex*,BoundaryLayerFan> _fans;
@@ -49,8 +56,8 @@ public:
iter end() {return _data.end();}
iterf beginf() {return _fans.begin();}
iterf endf() {return _fans.end();}
- BoundaryLayerColumns ()
- {
+ BoundaryLayerColumns ()
+ {
}
inline void addColumn (const SVector3 &dir, MVertex* v,
std::vector<MVertex*> _column,
@@ -60,7 +67,7 @@ public:
inline void addFan (MVertex *v, MEdge e1, MEdge e2, bool s){
_fans.insert(std::make_pair(v,BoundaryLayerFan(e1,e2,s)));
}
-
+
inline const BoundaryLayerData & getColumn ( MVertex *v, MEdge e ){
std::map<MVertex*,BoundaryLayerFan>::const_iterator it = _fans.find(v);
int N = getNbColumns (v) ;
@@ -74,7 +81,7 @@ public:
}
else Msg::Error("CANNOT HANDLE EMBEDDED LINES IN BLs");
}
-
+
inline edgeColumn getColumns ( MVertex *v1, MVertex *v2 , int side ){
Equal_Edge aaa;
MEdge e(v1,v2);
@@ -88,7 +95,7 @@ public:
// if (nbSides != 1)printf("I'm here %d sides\n",nbSides);
// Standard case, only two extruded columns from the two vertices
if (N1 == 1 && N2 == 1) return edgeColumn(getColumn(v1,0),getColumn(v2,0));
- // one fan on
+ // one fan on
if (nbSides == 1){
if (it1 != _fans.end() && it2 == _fans.end() ){
if (aaa(it1->second._e1,e))
@@ -101,8 +108,8 @@ public:
return edgeColumn(getColumn (v1,0),getColumn(v2,0));
else
return edgeColumn(getColumn (v1,0),getColumn(v2,N2-1));
- }
-
+ }
+
if (N1 == 1 || N2 == 2){
const BoundaryLayerData & c10 = getColumn(v1,0);
const BoundaryLayerData & c20 = getColumn(v2,0);
@@ -120,14 +127,14 @@ public:
Msg::Error ("Impossible Boundary Layer Configuration : one side and no fans");
// FIXME WRONG
- return N1 ? edgeColumn (getColumn (v1,0),getColumn(v1,0)) : edgeColumn (getColumn (v2,0),getColumn(v2,0));
+ return N1 ? edgeColumn (getColumn (v1,0),getColumn(v1,0)) : edgeColumn (getColumn (v2,0),getColumn(v2,0));
}
else if (nbSides == 2){
if (it1 == _fans.end() && it2 == _fans.end() ){
if (N1 != 2 || N2 != 2){
Msg::Error ("Impossible Boundary Layer Configuration");
// FIXME WRONG
- return N1 ? edgeColumn (getColumn (v1,0),getColumn(v1,0)) : edgeColumn (getColumn (v2,0),getColumn(v2,0));
+ return N1 ? edgeColumn (getColumn (v1,0),getColumn(v1,0)) : edgeColumn (getColumn (v2,0),getColumn(v2,0));
}
// printf("coucou side %d %d %d\n",side,N1,N2);
const BoundaryLayerData & c10 = getColumn(v1,0);
@@ -148,7 +155,7 @@ public:
Msg::Error ("Not yet Done in BoundaryLayerData");
}
}
-
+
inline int getNbColumns (MVertex *v) {
return _data.count (v);
}
@@ -158,14 +165,12 @@ public:
itm != _data.upper_bound(v); ++itm)
{
if (count++ == iColumn) return itm->second;
- }
+ }
}
void filterPoints();
};
-
BoundaryLayerColumns * buidAdditionalPoints2D (GFace *gf) ;
void buildMeshMetric(GFace *gf, double *uv, SMetric3 &m, double metric[3]);
-
#endif
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index 464b325127d49f3df38bbe92b37930362bd4bc34..b1fa1ad386627ca461a7afa62f260f8a617bf682 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -1,4 +1,7 @@
-
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
#include "meshMetric.h"
#include "meshGFaceOptimize.h"
@@ -48,7 +51,7 @@ meshMetric::meshMetric(GModel *gm){
}
}
}
- _octree = new MElementOctree(_elements);
+ _octree = new MElementOctree(_elements);
}
meshMetric::meshMetric(std::vector<MElement*> elements){
@@ -62,7 +65,7 @@ meshMetric::meshMetric(std::vector<MElement*> elements){
_elements.push_back(copy);
}
- _octree = new MElementOctree(_elements);
+ _octree = new MElementOctree(_elements);
}
void meshMetric::addMetric(int technique, simpleFunction<double> *fct, std::vector<double> parameters){
@@ -174,23 +177,23 @@ meshMetric::~meshMetric(){
void meshMetric::computeValues( v2t_cont adj){
- v2t_cont :: iterator it = adj.begin();
+ v2t_cont :: iterator it = adj.begin();
while (it != adj.end()) {
- std::vector<MElement*> lt = it->second;
+ std::vector<MElement*> lt = it->second;
MVertex *ver = it->first;
vals[ver]= (*_fct)(ver->x(),ver->y(),ver->z());
it++;
}
}
-// u = a + b(x-x_0) + c(y-y_0) + d(z-z_0)
+// u = a + b(x-x_0) + c(y-y_0) + d(z-z_0)
void meshMetric::computeHessian( v2t_cont adj){
int DIM = _dim + 1;
- for (int ITER=0;ITER<DIM;ITER++){
+ for (int ITER=0;ITER<DIM;ITER++){
v2t_cont :: iterator it = adj.begin();
while (it != adj.end()) {
- std::vector<MElement*> lt = it->second;
+ std::vector<MElement*> lt = it->second;
MVertex *ver = it->first;
while (lt.size() < 11) increaseStencil(ver,adj,lt); //<7
// if ( ver->onWhat()->dim() < _dim ){
@@ -209,9 +212,9 @@ void meshMetric::computeHessian( v2t_cont adj){
for(int i = 0; i < lt.size(); i++) {
MElement *e = lt[i];
int npts; IntPt *pts;
- SPoint3 p;
+ SPoint3 p;
e->getIntegrationPoints(0,&npts,&pts);
- if (ITER == 0){
+ if (ITER == 0){
SPoint3 p;
e->pnt(pts[0].pt[0],pts[0].pt[1],pts[0].pt[2],p);
b(i) = (*_fct)(p.x(),p.y(),p.z());
@@ -237,7 +240,7 @@ void meshMetric::computeHessian( v2t_cont adj){
if (ITER == 0){
double gr1 = result(1);
double gr2 = result(2);
- double gr3 = (_dim==2) ? 0.0:result(3);
+ double gr3 = (_dim==2) ? 0.0:result(3);
double norm = sqrt(gr1*gr1+gr2*gr2+gr3*gr3);
if (norm == 0.0 || _technique == meshMetric::HESSIAN) norm = 1.0;
grads[ver] = SVector3(gr1/norm,gr2/norm,gr3/norm);
@@ -289,26 +292,26 @@ void meshMetric::computeMetric(){
if (dist < _E && norm != 0.0){
double h = hmin*(hmax/hmin-1)*dist/_E + hmin;
double C = 1./(h*h) -1./(hmax*hmax);
- hlevelset(0,0) += C*gr(0)*gr(0)/norm ;
- hlevelset(1,1) += C*gr(1)*gr(1)/norm ;
- hlevelset(2,2) += C*gr(2)*gr(2)/norm ;
- hlevelset(1,0) = hlevelset(0,1) = C*gr(1)*gr(0)/norm ;
- hlevelset(2,0) = hlevelset(0,2) = C*gr(2)*gr(0)/norm ;
- hlevelset(2,1) = hlevelset(1,2) = C*gr(2)*gr(1)/norm ;
+ hlevelset(0,0) += C*gr(0)*gr(0)/norm ;
+ hlevelset(1,1) += C*gr(1)*gr(1)/norm ;
+ hlevelset(2,2) += C*gr(2)*gr(2)/norm ;
+ hlevelset(1,0) = hlevelset(0,1) = C*gr(1)*gr(0)/norm ;
+ hlevelset(2,0) = hlevelset(0,2) = C*gr(2)*gr(0)/norm ;
+ hlevelset(2,1) = hlevelset(1,2) = C*gr(2)*gr(1)/norm ;
}
// else if (dist > _E && dist < 2.*_E && norm != 0.0){
// double hmid = (hmax-hmin)/(1.*_E)*dist+(2.*hmin-1.*hmax);
// double C = 1./(hmid*hmid) -1./(hmax*hmax);
- // hlevelset(0,0) += C*gr(0)*gr(0)/norm ;
- // hlevelset(1,1) += C*gr(1)*gr(1)/norm ;
- // hlevelset(2,2) += C*gr(2)*gr(2)/norm ;
- // hlevelset(1,0) = hlevelset(0,1) = C*gr(1)*gr(0)/norm ;
- // hlevelset(2,0) = hlevelset(0,2) = C*gr(2)*gr(0)/norm ;
- // hlevelset(2,1) = hlevelset(1,2) = C*gr(2)*gr(1)/norm ;
+ // hlevelset(0,0) += C*gr(0)*gr(0)/norm ;
+ // hlevelset(1,1) += C*gr(1)*gr(1)/norm ;
+ // hlevelset(2,2) += C*gr(2)*gr(2)/norm ;
+ // hlevelset(1,0) = hlevelset(0,1) = C*gr(1)*gr(0)/norm ;
+ // hlevelset(2,0) = hlevelset(0,2) = C*gr(2)*gr(0)/norm ;
+ // hlevelset(2,1) = hlevelset(1,2) = C*gr(2)*gr(1)/norm ;
// }
H = hlevelset;
}
- //See paper Ducrot and Frey:
+ //See paper Ducrot and Frey:
//Anisotropic levelset adaptation for accurate interface capturing,
//ijnmf, 2010
else if (_technique == meshMetric::FREY ){
@@ -327,8 +330,8 @@ void meshMetric::computeMetric(){
hfrey(2,0) = hfrey(0,2) = C*gr(2)*gr(0)/(norm) + hessian(2,0)/epsGeom;
hfrey(2,1) = hfrey(1,2) = C*gr(2)*gr(1)/(norm) + hessian(2,1)/epsGeom;
// hfrey(0,0) += C*gr(0)*gr(0)/norm;
- // hfrey(1,1) += C*gr(1)*gr(1)/norm;
- // hfrey(2,2) += C*gr(2)*gr(2)/norm;
+ // hfrey(1,1) += C*gr(1)*gr(1)/norm;
+ // hfrey(2,2) += C*gr(2)*gr(2)/norm;
// hfrey(1,0) = hfrey(0,1) = gr(1)*gr(0)/(norm) ;
// hfrey(2,0) = hfrey(0,2) = gr(2)*gr(0)/(norm) ;
// hfrey(2,1) = hfrey(1,2) = gr(2)*gr(1)/(norm) ;
@@ -350,7 +353,7 @@ void meshMetric::computeMetric(){
SMetric3 metric;
if (signed_dist < _E && signed_dist > _E_moins && norm != 0.0){
- // first, compute the eigenvectors of the hessian, for a curvature-based metric
+ // first, compute the eigenvectors of the hessian, for a curvature-based metric
fullMatrix<double> eigvec_hessian(3,3);
fullVector<double> lambda_hessian(3);
hessian.eig(eigvec_hessian,lambda_hessian);
@@ -359,11 +362,11 @@ void meshMetric::computeMetric(){
ti.push_back(SVector3(eigvec_hessian(0,1),eigvec_hessian(1,1),eigvec_hessian(2,1)));
ti.push_back(SVector3(eigvec_hessian(0,2),eigvec_hessian(1,2),eigvec_hessian(2,2)));
- // compute the required characteristic length relative to the curvature and the distance to iso
+ // compute the required characteristic length relative to the curvature and the distance to iso
double maxkappa = std::max(std::max(fabs(lambda_hessian(0)),fabs(lambda_hessian(1))),fabs(lambda_hessian(2)));
// epsilon is set such that the characteristic element size in the tangent direction is h = 2pi R_min/_Np = sqrt(1/metric_maxmax) = sqrt(epsilon/kappa_max)
double un_sur_epsilon = maxkappa/((2.*3.14/_Np)*(2.*3.14/_Np));
- // metric curvature-based eigenvalues
+ // metric curvature-based eigenvalues
std::vector<double> eigenvals_curvature;
eigenvals_curvature.push_back(fabs(lambda_hessian(0))*un_sur_epsilon);
eigenvals_curvature.push_back(fabs(lambda_hessian(1))*un_sur_epsilon);
@@ -394,7 +397,7 @@ void meshMetric::computeMetric(){
for (int i=0;i<3;i++)
if (i!=grad_index) ti_index.push_back(i);
// std::cout << "gr tgr t1 t2 dots_tgr_gr_t1_t2 (" << gr(0) << "," << gr(1) << "," << gr(2) << ") (" << ti[grad_index](0) << "," << ti[grad_index](1) << "," << ti[grad_index](2) << ") (" << ti[ti_index[0]](0) << "," << ti[ti_index[0]](1) << "," << ti[ti_index[0]](2) << ") (" << ti[ti_index[1]](0) << "," << ti[ti_index[1]](1) << "," << ti[ti_index[1]](2) << ") " << fabs(dot(ti[grad_index],gr)) << " " << (dot(ti[grad_index],ti[ti_index[0]])) << " " << (dot(ti[grad_index],ti[ti_index[1]])) << std::endl;
-
+
// finally, creating the metric
std::vector<double> eigenvals;
eigenvals.push_back(std::min(std::max(eigenval_direction,metric_value_hmax),metric_value_hmin));// in gradient direction
@@ -415,13 +418,13 @@ void meshMetric::computeMetric(){
setOfMetrics[metricNumber].insert(std::make_pair(ver,metric));
setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
- it++;
+ it++;
}
if (_technique != meshMetric::EIGENDIRECTIONS && _technique!=meshMetric::EIGENDIRECTIONS_LINEARINTERP_H){
- //See paper Hachem and Coupez:
- //Finite element solution to handle complex heat and fluid flows in
+ //See paper Hachem and Coupez:
+ //Finite element solution to handle complex heat and fluid flows in
//industrial furnaces using the immersed volume technique, ijnmf, 2010
fullMatrix<double> V(3,3);
@@ -435,7 +438,7 @@ void meshMetric::computeMetric(){
// hessian.eig(Vhess,Shess);
// double h = hmin*(hmax/hmin-1)*dist/_E + hmin;
// double lam1 = Shess(0);
- // double lam2 = Shess(1);
+ // double lam2 = Shess(1);
// double lam3 = (_dim == 3)? Shess(2) : 1.;
// lambda1 = lam1;
// lambda2 = lam2/lambda1;
@@ -443,8 +446,8 @@ void meshMetric::computeMetric(){
// }
// else{
lambda1 = S(0);
- lambda2 = S(1);
- lambda3 = (_dim == 3)? S(2) : 1.;
+ lambda2 = S(1);
+ lambda3 = (_dim == 3)? S(2) : 1.;
//}
if (_technique == meshMetric::HESSIAN || (dist < _E && _technique == meshMetric::FREY)){
@@ -464,7 +467,7 @@ void meshMetric::computeMetric(){
setOfMetrics[metricNumber].insert(std::make_pair(ver,metric));
setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
- it++;
+ it++;
}
}
@@ -504,8 +507,8 @@ double meshMetric::operator() (double x, double y, double z, GEntity *ge) {
}
SPoint3 xyz(x,y,z), uvw;
double initialTol = MElement::getTolerance();
- MElement::setTolerance(1.e-4);
- MElement *e = _octree->find(x, y, z, _dim);
+ MElement::setTolerance(1.e-4);
+ MElement *e = _octree->find(x, y, z, _dim);
MElement::setTolerance(initialTol);
if (e){
e->xyz2uvw(xyz,uvw);
@@ -515,7 +518,7 @@ double meshMetric::operator() (double x, double y, double z, GEntity *ge) {
}
double value = e->interpolate(val,uvw[0],uvw[1],uvw[2]);
delete [] val;
- return value;
+ return value;
}
return 1.e22;
}
@@ -529,8 +532,8 @@ void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEnti
metr = SMetric3(1.e-22);
SPoint3 xyz(x,y,z), uvw;
double initialTol = MElement::getTolerance();
- MElement::setTolerance(1.e-4);
- MElement *e = _octree->find(x, y, z, _dim);
+ MElement::setTolerance(1.e-4);
+ MElement *e = _octree->find(x, y, z, _dim);
MElement::setTolerance(initialTol);
if (e){
@@ -560,13 +563,13 @@ void meshMetric::printMetric(const char* n){
FILE *f = fopen (n,"w");
fprintf(f,"View \"\"{\n");
- //std::map<MVertex*,SMetric3 >::const_iterator it = _hessian.begin();
+ //std::map<MVertex*,SMetric3 >::const_iterator it = _hessian.begin();
std::map<MVertex*,SMetric3>::const_iterator it= _nodalMetrics.begin();
//for (; it != _nodalMetrics.end(); ++it){
for (; it != _hessian.end(); ++it){
MVertex *v = it->first;
SMetric3 h = it->second;
- double lapl = h(0,0)+h(1,1)+h(2,2);
+ double lapl = h(0,0)+h(1,1)+h(2,2);
//fprintf(f, "SP(%g,%g,%g){%g};\n", it->first->x(),it->first->y(),it->first->z(),lapl);
fprintf(f,"TP(%g,%g,%g){%g,%g,%g,%g,%g,%g,%g,%g,%g};\n",
it->first->x(),it->first->y(),it->first->z(),
@@ -575,16 +578,16 @@ void meshMetric::printMetric(const char* n){
h(2,0),h(2,1),h(2,2)
);
- }
+ }
fprintf(f,"};\n");
- fclose (f);
+ fclose (f);
}
double meshMetric::getLaplacian (MVertex *v) {
MVertex *vNew = _vertexMap[v->getNum()];
std::map<MVertex*, SMetric3 >::const_iterator it = _hessian.find(vNew);
SMetric3 h = it->second;
- return h(0,0)+h(1,1)+h(2,2);
+ return h(0,0)+h(1,1)+h(2,2);
}
/*void meshMetric::curvatureContributionToMetric (){
@@ -597,14 +600,14 @@ double meshMetric::getLaplacian (MVertex *v) {
SMetric3 curvMetric = max_edge_curvature_metric((GVertex*)v->onWhat(),l);
it->second = intersection(it->second,curvMetric);
its->second = std::min(l,its->second);
- }
+ }
else if (v->onWhat()->dim() == 1){
double u,l;
v->getParameter(0,u);
SMetric3 curvMetric = max_edge_curvature_metric((GEdge*)v->onWhat(),u,l);
it->second = intersection(it->second,curvMetric);
its->second = std::min(l,its->second);
- }
+ }
}
}*/
diff --git a/Mesh/meshMetric.h b/Mesh/meshMetric.h
index 238338df20296714385b9b6a3772f93057852dbb..262f418716edb9cab590b98dc575a30fee414c81 100644
--- a/Mesh/meshMetric.h
+++ b/Mesh/meshMetric.h
@@ -1,10 +1,17 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef _MESH_METRIC_H_
#define _MESH_METRIC_H_
+
#include <map>
#include <algorithm>
#include "STensor3.h"
#include "Field.h"
#include "meshGFaceOptimize.h"
+
template <class scalar> class simpleFunction;
class MVertex;
class gLevelset;
@@ -38,10 +45,10 @@ class meshMetric: public Field {
nodalMetricTensor _nodalMetrics,_hessian;
nodalField _nodalSizes, _detMetric;
- std::map<int,nodalMetricTensor> setOfMetrics;
- std::map<int,nodalField> setOfSizes;
- std::map<int,nodalField> setOfDetMetric;
-
+ std::map<int,nodalMetricTensor> setOfMetrics;
+ std::map<int,nodalField> setOfSizes;
+ std::map<int,nodalField> setOfDetMetric;
+
public:
meshMetric(std::vector<MElement*> elements);
meshMetric(GModel *gm);
@@ -49,9 +56,9 @@ class meshMetric: public Field {
~meshMetric();
// compute a new metric and add it to the set of metrics
- // parameters[1] = lcmin (default : in global gmsh options CTX::instance()->mesh.lcMin)
- // parameters[2] = lcmax (default : in global gmsh options CTX::instance()->mesh.lcMax)
- // Available algorithms ("techniques"):
+ // parameters[1] = lcmin (default : in global gmsh options CTX::instance()->mesh.lcMin)
+ // parameters[2] = lcmax (default : in global gmsh options CTX::instance()->mesh.lcMax)
+ // Available algorithms ("techniques"):
// 1: fct is a LS, metric based on Coupez technique
// parameters[0] = thickness of the interface (mandatory)
// 2: metric based on the hessian of fct
@@ -73,7 +80,7 @@ class meshMetric: public Field {
void computeMetric() ;
void computeValues( v2t_cont adj);
void computeHessian( v2t_cont adj);
-
+
double getLaplacian (MVertex *v);
virtual bool isotropic () const {return false;}
virtual const char *getName()
@@ -88,9 +95,10 @@ class meshMetric: public Field {
// get metric at point(x,y,z) (previously computes intersection of metrics if not done yet)
virtual double operator() (double x, double y, double z, GEntity *ge=0) ;
virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
-
+
void printMetric(const char* n);
// export pos files of fct, fct gradients (fct is the lattest fct passed to meshMetric !!) and resulting metric (intersection of all computed metrics)
void exportInfo(const char *fileendname);
};
+
#endif
diff --git a/Numeric/BergotBasis.cpp b/Numeric/BergotBasis.cpp
index 636c6e500bce4d0c3666f85730385b10e18fe31c..bc8ce13da41f92770cebeafd83f486dae3c5fcfc 100644
--- a/Numeric/BergotBasis.cpp
+++ b/Numeric/BergotBasis.cpp
@@ -1,12 +1,18 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include "BergotBasis.h"
-BergotBasis::BergotBasis(int p): order(p) {
- // allocate function information and fill
+BergotBasis::BergotBasis(int p): order(p)
+{
+ // allocate function information and fill
iOrder = new int[size()];
jOrder = new int[size()];
kOrder = new int[size()];
-
+
legendre = new LegendrePolynomials(order);
int index = 0;
@@ -14,11 +20,11 @@ BergotBasis::BergotBasis(int p): order(p) {
for (int j=0;j<=order;j++) {
int mIJ = std::max(i,j);
for (int k=0;k<= (int) (order - mIJ);k++,index++) {
-
+
iOrder[index] = i;
jOrder[index] = j;
kOrder[index] = k;
-
+
if (jacobi.find(mIJ) == jacobi.end()) {
jacobi[mIJ] = new JacobiPolynomials(2*mIJ+2,0,order);
}
@@ -27,27 +33,26 @@ BergotBasis::BergotBasis(int p): order(p) {
}
}
-BergotBasis::~BergotBasis() {
-
+BergotBasis::~BergotBasis()
+{
delete [] iOrder;
delete [] jOrder;
delete [] kOrder;
-
+
delete legendre;
std::map<int,JacobiPolynomials*>::iterator jIter = jacobi.begin();
for (;jIter!=jacobi.end();++jIter) delete jIter->second;
}
-
-int BergotBasis::size() const {
+int BergotBasis::size() const
+{
return (2*order*order*order + 9*order*order + 13*order + 6)/6;
}
-
-void BergotBasis::f(double u, double v, double w, double* val) const {
-
+void BergotBasis::f(double u, double v, double w, double* val) const
+{
double uhat = (w == 1.) ? 1 : u/(1.-w);
-
+
std::vector<double> uFcts(order+1);
legendre->f(uhat,&(uFcts[0]));
@@ -65,8 +70,8 @@ void BergotBasis::f(double u, double v, double w, double* val) const {
double* wf = wFcts[a];
jIter->second->f(what,wf);
}
-
- // recombine to find shape function values
+
+ // recombine to find shape function values
int index = 0;
for (int i=0;i<=order;i++) {
@@ -79,7 +84,7 @@ void BergotBasis::f(double u, double v, double w, double* val) const {
}
}
}
-
+
jIter=jacobi.begin();
for (;jIter!=jacobi.end();jIter++) {
@@ -88,14 +93,15 @@ void BergotBasis::f(double u, double v, double w, double* val) const {
}
}
-void BergotBasis::df(double u, double v, double w, double grads[][3]) const {
+void BergotBasis::df(double u, double v, double w, double grads[][3]) const
+{
std::vector<double> uFcts(order+1);
legendre->f(u,&(uFcts[0]));
- std::vector<double> uGrads(order+1);
+ std::vector<double> uGrads(order+1);
legendre->df(u,&(uGrads[0]));
- std::vector<double> vFcts(order+1);
+ std::vector<double> vFcts(order+1);
legendre->f(v,&(vFcts[0]));
std::vector<double> vGrads(order+1);
@@ -105,17 +111,17 @@ void BergotBasis::df(double u, double v, double w, double grads[][3]) const {
std::map<int,double* > wFcts;
std::map<int,double* > wGrads;
std::map<int,JacobiPolynomials*>::const_iterator jIter=jacobi.begin();
-
+
for (;jIter!=jacobi.end();jIter++) {
int a = jIter->first;
- wFcts[a] = new double[order+1];
+ wFcts[a] = new double[order+1];
double* wf = wFcts[a];
jIter->second->f(w,wf);
wGrads[a] = new double[order+1];
double* wg = wGrads[a];
jIter->second->df(w,wg);
}
-
+
// now recombine to find the shape function gradients
int index = 0;
@@ -126,7 +132,7 @@ void BergotBasis::df(double u, double v, double w, double grads[][3]) const {
double* wf = (wFcts .find(mIJ))->second;
double* wg = (wGrads.find(mIJ))->second;
-
+
double wPowM2 = pow_int(1.-w,((int) mIJ) - 2);
double wPowM1 = wPowM2*(1.-w);
double wPowM0 = wPowM1*(1.-w);
@@ -142,8 +148,8 @@ void BergotBasis::df(double u, double v, double w, double grads[][3]) const {
grads[index][2] += v * uFcts[i] * vGrads[j] * wf[k] * wPowM2;
grads[index][2] -= mIJ * uFcts[i] * vFcts[j] * wf[k] * wPowM1;
grads[index][2] += 2 * uFcts[i] * vFcts[j] * wg[k] * wPowM0;
-
-
+
+
}
}
}
diff --git a/Numeric/BergotBasis.h b/Numeric/BergotBasis.h
index a1bc9be529d44882b127a0d5c68e8643ed5481b0..cecf7b1bbe7e6f0c7a2bee7eaacf82ac8bfbfb4e 100644
--- a/Numeric/BergotBasis.h
+++ b/Numeric/BergotBasis.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef BERGOTBASIS_H
#define BERGOTBASIS_H
@@ -17,13 +22,13 @@ class BergotBasis : public polynomialBasis {
private:
int order; //!< maximal order of surrounding functional spaces (on triangle / quad)
- int *iOrder; //!< order of \f$\hat \xi \f$ polynomial
- int *jOrder; //!< order of \f$\hat \eta \f$ polynomial
- int *kOrder; //!< order of \f$\hat \zeta \f$ polynomial
+ int *iOrder; //!< order of \f$\hat \xi \f$ polynomial
+ int *jOrder; //!< order of \f$\hat \eta \f$ polynomial
+ int *kOrder; //!< order of \f$\hat \zeta \f$ polynomial
- //! list of Legendre polynomials up to order p
+ //! list of Legendre polynomials up to order p
LegendrePolynomials* legendre;
-
+
//! list of Jacobi polynomials up to order p in function of index i (\f$ \alpha = 2*i + 2\f$)
std::map<int,JacobiPolynomials*> jacobi;
diff --git a/Numeric/GaussIntegration.cpp b/Numeric/GaussIntegration.cpp
index 848ea626247973722530b66871077bfc2d413252..61531226ced6ca1f952f8f34a254d32d833bad63 100644
--- a/Numeric/GaussIntegration.cpp
+++ b/Numeric/GaussIntegration.cpp
@@ -1,8 +1,13 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include "GaussIntegration.h"
#include "GmshDefines.h"
-static void pts2fullMatrix(int npts, IntPt *pts, fullMatrix<double> &pMatrix,
- fullMatrix<double> &wMatrix)
+static void pts2fullMatrix(int npts, IntPt *pts, fullMatrix<double> &pMatrix,
+ fullMatrix<double> &wMatrix)
{
pMatrix.resize(npts,3);
wMatrix.resize(npts,1);
@@ -20,7 +25,7 @@ void gaussIntegration::getTriangle(int order, fullMatrix<double> &pts,
pts2fullMatrix(getNGQTPts(order),getGQTPts(order),pts,weights);
}
-void gaussIntegration::getLine(int order, fullMatrix<double> &pts,
+void gaussIntegration::getLine(int order, fullMatrix<double> &pts,
fullMatrix<double> &weights)
{
pts2fullMatrix(getNGQLPts(order),getGQLPts(order),pts,weights);
diff --git a/Numeric/GaussJacobi1D.cpp b/Numeric/GaussJacobi1D.cpp
index 56ce75f0a06613c48b5a7ee967d91e0460a1e2b9..5bad8a1e831a45369e9daa6f5f46cb60c73dc424 100644
--- a/Numeric/GaussJacobi1D.cpp
+++ b/Numeric/GaussJacobi1D.cpp
@@ -1,8 +1,13 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
// code to generate this file (using alglib) :
#if 0
//alglib
#include "ap.h"
-#include "gq.h"
+#include "gq.h"
int main(int argc, char **argv)
{
int maxA = 4;
diff --git a/Numeric/GaussJacobi1D.h b/Numeric/GaussJacobi1D.h
index 769e8fd6edf3ec8af22357caa85e77770dc03942..84a639999a416ff7673073ec4d9cd71d9f6d405a 100644
--- a/Numeric/GaussJacobi1D.h
+++ b/Numeric/GaussJacobi1D.h
@@ -1,4 +1,11 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef _GAUSS_JACOBI_1D_H_
#define _GAUSS_JACOBI_1D_H_
-void getGaussJacobiQuadrature(int a, int b, int n, double **pt, double **wt)
+
+void getGaussJacobiQuadrature(int a, int b, int n, double **pt, double **wt);
+
#endif
diff --git a/Numeric/jacobiPolynomials.cpp b/Numeric/jacobiPolynomials.cpp
index 911a1d0ed38080a87eb0b454c7ec75bf8ed6d623..7a03c9c9b1f7d37ab1f1ac8bc8b05aa0312d5edf 100644
--- a/Numeric/jacobiPolynomials.cpp
+++ b/Numeric/jacobiPolynomials.cpp
@@ -1,43 +1,47 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include "jacobiPolynomials.h"
-inline double Pochhammer(double x,int n) {
+inline double Pochhammer(double x,int n)
+{
double result(1.);
for (int i=0;i<n;i++) result *= (x+i);
return result;
}
-
JacobiPolynomials::JacobiPolynomials(double a, double b, int o):
alpha(a),beta(b),n(o),alphaPlusBeta(a+b),a2MinusB2(a*a-b*b) {}
-JacobiPolynomials::~JacobiPolynomials() {;}
-void JacobiPolynomials::f(double u, double *val) const {
+JacobiPolynomials::~JacobiPolynomials() {;}
+void JacobiPolynomials::f(double u, double *val) const
+{
val[0] = 1.;
if (n>=1) val[1] = 0.5*(2.*(alpha+1.) + (alphaPlusBeta + 2.)*(u-1.));
-
+
for (int i=1;i<n;i++) {
-
+
double ii = (double) i;
double twoI = 2.*ii;
-
+
double a1i = 2.*(ii+1.)*(ii+alphaPlusBeta+1.)*(twoI+alphaPlusBeta);
double a2i = (twoI+alphaPlusBeta+1.)*(a2MinusB2);
double a3i = Pochhammer(twoI+alphaPlusBeta,3);
double a4i = 2.*(ii+alpha)*(ii+beta)*(twoI+alphaPlusBeta+2.);
-
+
val[i+1] = ((a2i + a3i * u)* val[i] - a4i * val[i-1])/a1i;
-
- }
+ }
}
-
-void JacobiPolynomials::df(double u, double *val) const {
-
+void JacobiPolynomials::df(double u, double *val) const
+{
std::vector<double> tmp(n+1);
f(u,&(tmp[0]));
-
+
val[0] = 0;
if (n>=1) val[1] = 0.5*(alphaPlusBeta + 2.);
@@ -47,7 +51,7 @@ void JacobiPolynomials::df(double u, double *val) const {
double g2 = aa*(1.-u*u);
double g1 = ii*(alpha - beta - aa*u);
double g0 = 2.*(ii+alpha)*(ii+beta);
-
+
val[i] = (g1 * tmp[i] + g0 * tmp[i-1])/g2;
}
}
diff --git a/Numeric/jacobiPolynomials.h b/Numeric/jacobiPolynomials.h
index 31f88d644aae7d348f19ef330b75929da2545a3c..7b88a950c1be2b35727ca225166e8ef7f3c5c1f5 100644
--- a/Numeric/jacobiPolynomials.h
+++ b/Numeric/jacobiPolynomials.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef JACOBIPOLYNOMIALS_H
#define JACOBIPOLYNOMIALS_H
diff --git a/Numeric/legendrePolynomials.cpp b/Numeric/legendrePolynomials.cpp
index 397c59e0f504c16ac68fcd9d2178dd20a01090dd..3968e53c4eb48dee71f19fbfac5235cf273030c6 100644
--- a/Numeric/legendrePolynomials.cpp
+++ b/Numeric/legendrePolynomials.cpp
@@ -1,34 +1,41 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include "legendrePolynomials.h"
LegendrePolynomials::LegendrePolynomials(int o): n(o) {}
LegendrePolynomials::~LegendrePolynomials() {;}
-void LegendrePolynomials::f(double u, double *val) const {
+void LegendrePolynomials::f(double u, double *val) const
+{
val[0] = 1;
-
+
for (int i=0;i<n;i++) {
-
+
double a1i = i+1;
double a3i = 2.*i+1;
double a4i = i;
-
+
val[i+1] = a3i * u * val[i];
if (i>0) val[i+1] -= a4i * val[i-1];
val[i+1] /= a1i;
}
}
-void LegendrePolynomials::df(double u, double *val) const {
+void LegendrePolynomials::df(double u, double *val) const
+{
std::vector<double> tmp(n+1);
f(u,&(tmp[0]));
- val[0] = 0;
+ val[0] = 0;
double g2 = (1.-u*u);
-
+
for (int i=1;i<=n;i++) {
double g1 = - u*i;
double g0 = (double) i;
-
+
val[i] = (g1 * tmp[i] + g0 * tmp[i-1])/g2;
}
}
diff --git a/Numeric/legendrePolynomials.h b/Numeric/legendrePolynomials.h
index ab67139a6feec347a7187246edddc350d5b3b111..ce7dd998b99f72adec4af6ca31e93aae01f6caa6 100644
--- a/Numeric/legendrePolynomials.h
+++ b/Numeric/legendrePolynomials.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef LEGENDREPOLYNOMIALS_H
#define LEGENDREPOLYNOMIALS_H
diff --git a/Qt/GLWidget.cpp b/Qt/GLWidget.cpp
index 333e5476844190e5441b5068565139ec9d868e1d..82c880ab66525f8c53f7ac41126b6c13767da120 100644
--- a/Qt/GLWidget.cpp
+++ b/Qt/GLWidget.cpp
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include <QtGui>
#include <QtOpenGL>
#include <math.h>
diff --git a/Qt/GLWidget.h b/Qt/GLWidget.h
index 18e90dd025c4523712cafea658015bc3e2bd5b3a..0f5dd542e261a2f3dff6f6c2873563b48bf645fb 100644
--- a/Qt/GLWidget.h
+++ b/Qt/GLWidget.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef _GLWIDGET_H_
#define _GLWIDGET_H_
diff --git a/Qt/graphicWindow.cpp b/Qt/graphicWindow.cpp
index e75716872ee4f9ac26cecb8ea2b842785fc8a6d6..de1057bc5664f6dc2950a5c92eabb7d28fd06e18 100644
--- a/Qt/graphicWindow.cpp
+++ b/Qt/graphicWindow.cpp
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include <QtGui>
#include "GLWidget.h"
#include "graphicWindow.h"
@@ -5,7 +10,7 @@
graphicWindow::graphicWindow()
{
_glWidget = new GLWidget();
-
+
_xSlider = new QSlider(Qt::Vertical);
_xSlider->setRange(0, 360 * 16);
_xSlider->setSingleStep(16);
@@ -13,10 +18,10 @@ graphicWindow::graphicWindow()
_xSlider->setTickInterval(15 * 16);
_xSlider->setTickPosition(QSlider::TicksRight);
_xSlider->setValue(15 * 16);
-
+
connect(_xSlider, SIGNAL(valueChanged(int)), _glWidget, SLOT(setXRotation(int)));
connect(_glWidget, SIGNAL(xRotationChanged(int)), _xSlider, SLOT(setValue(int)));
-
+
QHBoxLayout *mainLayout = new QHBoxLayout;
mainLayout->addWidget(_glWidget);
mainLayout->addWidget(_xSlider);
diff --git a/Qt/graphicWindow.h b/Qt/graphicWindow.h
index 74b08cfd84e83d0b6a210bbafc7ae5be9adf8454..3e05cd7de9a0b2adbc51aba9159ce1c636c1bc94 100644
--- a/Qt/graphicWindow.h
+++ b/Qt/graphicWindow.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef _GRAPHIC_WINDOW_H_
#define _GRAPHIC_WINDOW_H_
diff --git a/Solver/SElement.cpp b/Solver/SElement.cpp
index 746c7d44c3c2ba8360da046a1cafa796784f9245..fc02f393e5e47b577e4e7a90833cb3b2c370ccf2 100644
--- a/Solver/SElement.cpp
+++ b/Solver/SElement.cpp
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include "SElement.h"
// FIXME: this will change in the future (the base SElement should no
diff --git a/Solver/STensor33.cpp b/Solver/STensor33.cpp
index 96582fe0ae878a93f54296fd1ef93b9e3b8a5cc2..1335b7329bbd616367028fcc361a75d74c75b090 100644
--- a/Solver/STensor33.cpp
+++ b/Solver/STensor33.cpp
@@ -1,20 +1,29 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Eric Bechet
+//
+
#include "STensor33.h"
void STensor33::print (const char *s) const
{
char format[512];
const char l[256]="%12.5E %12.5E %12.5E \n";
- sprintf (format, " tensor3 %s : \n %s %s %s \n %s %s %s \n %s %s %s \n",s, l,l,l, l,l,l, l,l,l);
+ sprintf (format, " tensor3 %s : \n %s %s %s \n %s %s %s \n %s %s %s \n",
+ s, l,l,l, l,l,l, l,l,l);
printf(format,s,_val[ 0],_val[ 1],_val[ 2],
_val[ 3],_val[ 4],_val[ 5],
_val[ 6],_val[ 7],_val[ 8],
-
+
_val[ 9],_val[10],_val[11],
_val[12],_val[13],_val[14],
_val[15],_val[16],_val[17],
-
+
_val[18],_val[19],_val[20],
_val[21],_val[22],_val[23],
_val[24],_val[25],_val[26]);
}
-
diff --git a/Solver/STensor33.h b/Solver/STensor33.h
index f14ee8ac9071f53def48a92c79fc12ba68c2816f..0e32fc86abef15d43170206e5f16c996078fe4ff 100644
--- a/Solver/STensor33.h
+++ b/Solver/STensor33.h
@@ -1,3 +1,12 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Eric Bechet
+//
+
#ifndef _STENSOR33_H_
#define _STENSOR33_H_
@@ -5,7 +14,6 @@
#include "fullMatrix.h"
#include "Numeric.h"
-
// concrete class for general 3rd-order tensor in three-dimensional space
class STensor33 {
@@ -184,7 +192,4 @@ inline double operator*(const STensor33 &m, const STensor33 &t)
return val;
}
-
-
#endif
-
diff --git a/Solver/STensor43.cpp b/Solver/STensor43.cpp
index 90904f64c9954a6040fa88290fdbce7e8b91b1bc..dc33abe653adb38303a6a7902de63507bd4fe79c 100644
--- a/Solver/STensor43.cpp
+++ b/Solver/STensor43.cpp
@@ -1,3 +1,11 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Eric Bechet
+//
#include "STensor43.h"
@@ -9,15 +17,12 @@ void STensor43::print (const char *s) const
printf(format,s,_val[ 0],_val[ 1],_val[ 2], _val[ 3],_val[ 4],_val[ 5], _val[ 6],_val[ 7],_val[ 8],
_val[ 9],_val[10],_val[11], _val[12],_val[13],_val[14], _val[15],_val[16],_val[17],
_val[18],_val[19],_val[20], _val[21],_val[22],_val[23], _val[24],_val[25],_val[26],
-
+
_val[27],_val[28],_val[29], _val[30],_val[31],_val[32], _val[33],_val[34],_val[35],
_val[36],_val[37],_val[38], _val[39],_val[40],_val[41], _val[42],_val[43],_val[44],
_val[45],_val[46],_val[47], _val[48],_val[49],_val[50], _val[51],_val[52],_val[53],
-
+
_val[54],_val[55],_val[56], _val[57],_val[58],_val[59], _val[60],_val[61],_val[62],
_val[63],_val[64],_val[65], _val[66],_val[67],_val[68], _val[69],_val[70],_val[71],
_val[72],_val[73],_val[74], _val[75],_val[76],_val[77], _val[78],_val[79],_val[80]);
}
-
-
-
diff --git a/Solver/STensor43.h b/Solver/STensor43.h
index c764fde277f0aeba8b4fb61c0c364ec7e42a9fd6..30563863025355147ebe9aa54eb3eefc597371db 100644
--- a/Solver/STensor43.h
+++ b/Solver/STensor43.h
@@ -1,3 +1,11 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Eric Bechet
+//
#ifndef _STENSOR43_H_
#define _STENSOR43_H_
@@ -6,7 +14,6 @@
#include "fullMatrix.h"
#include "Numeric.h"
-
// concrete class for general 3rd-order tensor in three-dimensional space
class STensor43 {
@@ -268,5 +275,5 @@ inline double operator*( const STensor43 &m , const STensor43 &t)
val+=m(i,j,k,l)*t(l,k,j,i);
return val;
}
-#endif
+#endif
diff --git a/Solver/dofManager.cpp b/Solver/dofManager.cpp
index a2a6c14253bfd4f3a088bcce07a1f29c41f6116a..aefd1b566eef5d593dd652c414701eda5f2300e0 100644
--- a/Solver/dofManager.cpp
+++ b/Solver/dofManager.cpp
@@ -1,8 +1,15 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include <dofManager.h>
#include "GmshConfig.h"
+
#ifdef HAVE_MPI
#include "mpi.h"
#endif
+
template<>
void dofManager<double>::scatterSolution()
{
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 68c30abedfc9b98895e79734d27dcd176b485641..72f5ce82e973b002fd2dcfcfa1101a12654994e1 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef _FUNCTION_SPACE_H_
#define _FUNCTION_SPACE_H_
diff --git a/Solver/groupOfElements.cpp b/Solver/groupOfElements.cpp
index 5da5342ab5ee6f980c352c6f68d4c8c5c9510f5d..4d3a68d709d48530d312cea0fb6823a6684be156 100644
--- a/Solver/groupOfElements.cpp
+++ b/Solver/groupOfElements.cpp
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include "groupOfElements.h"
#include "GModel.h"
#include "GEntity.h"
diff --git a/Solver/groupOfElements.h b/Solver/groupOfElements.h
index 42289900ec44237c7ff050dbc225fdb17a6d5199..4ade2d7072dd842d3beec4495d3def0f0044f13c 100644
--- a/Solver/groupOfElements.h
+++ b/Solver/groupOfElements.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef _GROUPOFELEMENTS_H_
#define _GROUPOFELEMENTS_H_
diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp
index 80094e0b0d5a13d47a30c6fa9c91bf064c6160af..527cb74eb75f27b4147b456592d787fe4904c039 100644
--- a/Solver/linearSystemPETSc.cpp
+++ b/Solver/linearSystemPETSc.cpp
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include "GmshConfig.h"
#if defined(HAVE_PETSC)
#include "petsc.h"
@@ -13,9 +18,8 @@ template class linearSystemPETSc<double>;
template class linearSystemPETSc<std::complex<double> >;
#endif
-
-
-void linearSystemPETScBlockDouble::_kspCreate() {
+void linearSystemPETScBlockDouble::_kspCreate()
+{
KSPCreate(_sequential ? PETSC_COMM_SELF : PETSC_COMM_WORLD, &_ksp);
if (this->_parameters.count("petscPrefix"))
KSPAppendOptionsPrefix(_ksp, this->_parameters["petscPrefix"].c_str());
@@ -23,7 +27,8 @@ void linearSystemPETScBlockDouble::_kspCreate() {
_kspAllocated = true;
}
-void linearSystemPETScBlockDouble::addToMatrix(int row, int col, const fullMatrix<double> &val)
+void linearSystemPETScBlockDouble::addToMatrix(int row, int col,
+ const fullMatrix<double> &val)
{
if (!_entriesPreAllocated)
preAllocateEntries();
@@ -58,7 +63,8 @@ void linearSystemPETScBlockDouble::addToMatrix(int row, int col, const fullMatri
#endif
}
-void linearSystemPETScBlockDouble::addToRightHandSide(int row, const fullMatrix<double> &val)
+void linearSystemPETScBlockDouble::addToRightHandSide(int row,
+ const fullMatrix<double> &val)
{
for (int ii = 0; ii < _blockSize; ii++) {
PetscInt i = row * _blockSize + ii;
@@ -67,12 +73,14 @@ void linearSystemPETScBlockDouble::addToRightHandSide(int row, const fullMatrix<
}
}
-void linearSystemPETScBlockDouble::getFromMatrix(int row, int col, fullMatrix<double> &val ) const
+void linearSystemPETScBlockDouble::getFromMatrix(int row, int col,
+ fullMatrix<double> &val ) const
{
Msg::Error("getFromMatrix not implemented for PETSc");
}
-void linearSystemPETScBlockDouble::getFromRightHandSide(int row, fullMatrix<double> &val) const
+void linearSystemPETScBlockDouble::getFromRightHandSide(int row,
+ fullMatrix<double> &val) const
{
for (int i = 0; i < _blockSize; i++) {
int ii = row*_blockSize +i;
@@ -95,7 +103,6 @@ void linearSystemPETScBlockDouble::addToSolution(int row, const fullMatrix<doubl
}
}
-
void linearSystemPETScBlockDouble::getFromSolution(int row, fullMatrix<double> &val) const
{
for (int i = 0; i < _blockSize; i++) {
@@ -123,7 +130,8 @@ void linearSystemPETScBlockDouble::allocate(int nbRows)
MatSetSizes(_a,nbRows * _blockSize, nbRows * _blockSize, PETSC_DETERMINE, PETSC_DETERMINE);
if (Msg::GetCommSize() > 1 && !_sequential) {
MatSetType(_a, MATMPIBAIJ);
- } else {
+ }
+ else {
MatSetType(_a, MATSEQBAIJ);
}
if (_parameters.count("petscPrefix"))
@@ -186,7 +194,9 @@ int linearSystemPETScBlockDouble::systemSolve()
KSPSolve(_ksp, _b, _x);
return 1;
}
-void linearSystemPETScBlockDouble::insertInSparsityPattern (int i, int j) {
+
+void linearSystemPETScBlockDouble::insertInSparsityPattern (int i, int j)
+{
i -= _localRowStart;
if (i<0 || i>= _localSize) return;
_sparsity.insertEntry (i,j);
@@ -256,7 +266,6 @@ void linearSystemPETScBlockDouble::zeroSolution()
}
}
-
linearSystemPETScBlockDouble::linearSystemPETScBlockDouble(bool sequential)
{
_entriesPreAllocated = false;
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index af30b617160abce860e13d392b24c3963decd6cc..b13a884e1ae79f11c131892a9332a36b9cade51a 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#include "Context.h"
#include "multiscaleLaplace.h"
#include "GmshConfig.h"
@@ -32,12 +37,12 @@ struct compareRotatedPoints {
}
bool operator ( ) (const SPoint2 &p1, const SPoint2 &p2) const {
//sort from left (x=-1) to right (sin=0), cos(=1)
- 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);
+ 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);
+ 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;
}
@@ -52,7 +57,7 @@ struct sort_pred {
}
};
-static void sort_centers_dist(std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers,
+static void sort_centers_dist(std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > ¢ers,
const SPoint2 &leftP){
std::vector<std::pair<SPoint2,multiscaleLaplaceLevel*> > myCenters(centers);
@@ -87,7 +92,7 @@ static void sort_centers_dist(std::vector<std::pair<SPoint2,multiscaleLaplaceLev
//--------------------------------------------------------------
static int intersection_segments_b (SPoint2 &p1, SPoint2 &p2,
- SPoint2 &q1, SPoint2 &q2,
+ SPoint2 &q1, SPoint2 &q2,
double x[2]){
double P1[2] = {p1.x(),p1.y()};
double P2[2] = {p2.x(),p2.y()};
@@ -111,7 +116,7 @@ static void recur_connect (MVertex *v,
if (touched.find(v) != touched.end())return;
touched.insert(v);
- for (std::multimap <MVertex*,MEdge>::iterator it = v2e.lower_bound(v);
+ for (std::multimap <MVertex*,MEdge>::iterator it = v2e.lower_bound(v);
it != v2e.upper_bound(v) ; ++it){
group.insert(it->second);
for (int i=0;i<it->second.getNumVertices();++i){
@@ -121,10 +126,10 @@ static void recur_connect (MVertex *v,
}
//--------------------------------------------------------------
-static void connected_bounds (std::vector<MElement*> &elements,
+static void connected_bounds (std::vector<MElement*> &elements,
std::vector<std::vector<MEdge> > &boundaries){
- std::vector<MEdge> bEdges;
+ std::vector<MEdge> bEdges;
for(unsigned int i = 0; i < elements.size(); i++){
for(int j = 0; j < elements[i]->getNumEdges(); j++){
MEdge me = elements[i]->getEdge(j);
@@ -132,9 +137,9 @@ static void connected_bounds (std::vector<MElement*> &elements,
bEdges.push_back(me);
else
bEdges.erase(std::find(bEdges.begin(), bEdges.end(),me));
- }
- }
-
+ }
+ }
+
std::multimap<MVertex*,MEdge> v2e;
for (unsigned int i = 0; i < bEdges.size(); ++i){
for (int j=0;j<bEdges[i].getNumVertices();j++){
@@ -171,24 +176,24 @@ static double getSizeBB(std::vector<MEdge> &me){
bb+=pt1;
bb+=pt2;
}
-
+
//double H = norm(SVector3(bb.max(), bb.min()));
//printf("H=%g \n", H);
obbox = SOrientedBoundingBox::buildOBB(vertices);
double H = obbox.getMaxSize();
-
- return H;
+
+ return H;
}
//--------------------------------------------------------------
-static void ordering_dirichlet(std::vector<MElement*> &elements,
+static void ordering_dirichlet(std::vector<MElement*> &elements,
std::vector<std::pair<MVertex*,double> > &dirichletNodes){
//finding all boundaries
std::vector<std::vector<MEdge> > boundaries;
connected_bounds(elements,boundaries);
-
+
//largest boundary is dirichlet boundary
std::vector<MEdge> dirichletEdges;
double maxSize = 0.0;
@@ -196,7 +201,7 @@ static void ordering_dirichlet(std::vector<MElement*> &elements,
std::vector<MEdge> iBound = boundaries[i];
double size = getSizeBB(iBound);
if (size > maxSize) {
- dirichletEdges = iBound;
+ dirichletEdges = iBound;
maxSize = size;
}
}
@@ -207,14 +212,14 @@ static void ordering_dirichlet(std::vector<MElement*> &elements,
double tot_length = 0.0;
for(unsigned int i = 0; i < dirichletEdges.size(); i++ ){
MVertex *v0 = dirichletEdges[i].getVertex(0);
- MVertex *v1 = dirichletEdges[i].getVertex(1);
- double len = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ MVertex *v1 = dirichletEdges[i].getVertex(1);
+ double len = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
(v0->y() - v1->y()) * (v0->y() - v1->y()) +
- (v0->z() - v1->z()) * (v0->z() - v1->z())) ;
+ (v0->z() - v1->z()) * (v0->z() - v1->z())) ;
tot_length += len;
temp.push_back(dirichletEdges[i]);
}
-
+
dirichletNodes.push_back(std::make_pair(dirichletEdges[0].getVertex(0),0.0));
MVertex *current_v = dirichletEdges[0].getVertex(1);
temp.erase(temp.begin());
@@ -228,7 +233,7 @@ static void ordering_dirichlet(std::vector<MElement*> &elements,
found = true;
current_v = v1;
temp.erase(itl);
- double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ 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);
@@ -239,7 +244,7 @@ static void ordering_dirichlet(std::vector<MElement*> &elements,
found = true;
current_v = v0;
temp.erase(itl);
- double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ 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);
@@ -248,13 +253,13 @@ static void ordering_dirichlet(std::vector<MElement*> &elements,
}
}
if(!found) return ;
- }
+ }
return;
}
//--------------------------------------------------------------
static int intersection_segments (SPoint2 &p1, SPoint2 &p2,
- SPoint2 &q1, SPoint2 &q2,
+ SPoint2 &q1, SPoint2 &q2,
double x[2]){
double A[2][2];
A[0][0] = p2.x()-p1.x();
@@ -265,13 +270,13 @@ static int intersection_segments (SPoint2 &p1, SPoint2 &p2,
sys2x2(A,b,x);
return (x[0] >= 0.0 && x[0] <= 1. &&
- x[1] >= 0.0 && x[1] <= 1.);
-
+ x[1] >= 0.0 && x[1] <= 1.);
+
}
//--------------------------------------------------------------
static void recur_compute_centers_ (double R, double a1, double a2,
multiscaleLaplaceLevel * root, int &nbElems ){
-
+
nbElems += root->elements.size();
root->radius = R;
@@ -282,12 +287,12 @@ static void recur_compute_centers_ (double R, double a1, double a2,
double fact = 2.5;
SPoint2 PL (fact*R*cos(a1),fact*R*sin(a1));
SPoint2 PR (fact*R*cos(a2),fact*R*sin(a2));
-
+
std::vector<SPoint2> centersChild;
centersChild.clear();
for (unsigned int i=0;i<root->children.size();i++){
- multiscaleLaplaceLevel* m = root->children[i];
- centers.push_back(std::make_pair(m->center,m));
+ multiscaleLaplaceLevel* m = root->children[i];
+ centers.push_back(std::make_pair(m->center,m));
m->radius = 0.0;
for (std::map<MVertex*,SPoint2>::iterator it = m->coordinates.begin();
it != m->coordinates.end() ; ++it){
@@ -297,8 +302,8 @@ static void recur_compute_centers_ (double R, double a1, double a2,
}
centersChild.push_back(m->center);
}
-
- //add the center of real holes ...
+
+ //add the center of real holes ...
std::vector<std::vector<MEdge> > boundaries;
connected_bounds(root->elements, boundaries);
int added = 0;
@@ -328,26 +333,26 @@ static void recur_compute_centers_ (double R, double a1, double a2,
(c.y() - p.y())*(c.y() - p.y()));
if (dist < 0.5*rad) newCenter = false;//0.6
}
-
+
if (std::abs(rad/root->radius) < 0.65 && std::abs(rad) < 0.95 && newCenter){//0.6
added++;
- centers.push_back(std::make_pair(c,zero));
+ centers.push_back(std::make_pair(c,zero));
}
}
}
if (added != toadd && root->children.size()> 0) {
- printf("!!!!!!!! ARG added = %d != %d (bounds=%d, child=%d)\n", added,
- (int)(boundaries.size() - 1 - root->children.size()), (int)boundaries.size(),
+ printf("!!!!!!!! ARG added = %d != %d (bounds=%d, child=%d)\n", added,
+ (int)(boundaries.size() - 1 - root->children.size()), (int)boundaries.size(),
(int)root->children.size());
}
//sort centers from left to right
std::sort(centers.begin(),centers.end(), sort_pred(PL,PR));
-
+
//sort from distances
//sort_centers_dist(centers, PL);
- centers.insert(centers.begin(), std::make_pair(PL,zero));
+ 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++){
@@ -355,8 +360,8 @@ static void recur_compute_centers_ (double R, double a1, double a2,
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);
}
}
@@ -389,9 +394,9 @@ static void recur_cut_edges_ (multiscaleLaplaceLevel *root,
std::map<MVertex *, SPoint2>::iterator it1 = root->coordinates.find(it->getVertex(1));
if (it0 != root->coordinates.end() && it1 != root->coordinates.end()){
SPoint2 q1 = root->coordinates[it->getVertex(0)];
- SPoint2 q2 = root->coordinates[it->getVertex(1)];
+ SPoint2 q2 = root->coordinates[it->getVertex(1)];
double x[2];
- int inters = intersection_segments (p1,p2,q1,q2,x);
+ 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(),
@@ -424,7 +429,7 @@ static void recur_cut_elements_ (multiscaleLaplaceLevel * root,
for (unsigned int i = 0; i < root->elements.size(); i++){
MVertex *c[3] = {0,0,0};
for (int j=0;j<3;j++){
- MEdge ed = root->elements[i]->getEdge(j);
+ MEdge ed = root->elements[i]->getEdge(j);
std::map<MEdge,MVertex*,Less_Edge> :: iterator it = cutEdges.find(ed);
if (it != cutEdges.end()){
c[j] = it->second;
@@ -515,7 +520,7 @@ 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<MEdge,Less_Edge> &theCut,
std::set<MElement*> &leftSet){
if (leftSet.find(e) != leftSet.end())return;
@@ -530,7 +535,7 @@ static void recur_leftCut_ (MElement *e,
}
}
}
-
+
}
//--------------------------------------------------------------
@@ -579,8 +584,8 @@ static void keepConnected(std::vector<MElement*> &goodSize, std::vector<MElement
connectedRegions (goodSize,regGoodSize);
if (regGoodSize.size() > 0){
int index=0;
- int maxSize= regGoodSize[0].size();
- for (unsigned int i=1;i< regGoodSize.size() ; i++){
+ int maxSize= regGoodSize[0].size();
+ for (unsigned int i=1;i< regGoodSize.size() ; i++){
int size = regGoodSize[i].size();
if(size > maxSize){
index = i;
@@ -588,7 +593,7 @@ static void keepConnected(std::vector<MElement*> &goodSize, std::vector<MElement
}
}
goodSize.clear();
- for (unsigned int i=0;i< regGoodSize.size() ; i++){
+ 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());
}
@@ -596,8 +601,8 @@ static void keepConnected(std::vector<MElement*> &goodSize, std::vector<MElement
}
//--------------------------------------------------------------
static void recur_cut_ (double R, double a1, double a2,
- multiscaleLaplaceLevel * root,
- std::vector<MElement *> &left,
+ multiscaleLaplaceLevel * root,
+ std::vector<MElement *> &left,
std::vector<MElement *> &right){
SPoint2 PL (R*cos(a1),R*sin(a1));
@@ -608,7 +613,7 @@ static void recur_cut_ (double R, double a1, double a2,
(PL.y()-PR.y())*(PL.y()-PR.y()));
SPoint2 farLeft (0.5*(PL.x()+PR.x()) - (PR.y()-PL.y())/d ,
0.5*(PL.y()+PR.y()) + (PR.x()-PL.x())/d );
-
+
for (unsigned int i = 0; i < root->elements.size(); i++){
SPoint2 pp (0,0);
for (int j=0; j<root->elements[i]->getNumVertices(); j++){
@@ -618,20 +623,20 @@ static void recur_cut_ (double R, double a1, double a2,
int nbIntersect = 0;
for (unsigned int j = 0; j < centers.size() - 1; j++){
double x[2];
- nbIntersect += intersection_segments (centers[j].first,centers[j+1].first,pp,farLeft,x);
+ nbIntersect += intersection_segments (centers[j].first,centers[j+1].first,pp,farLeft,x);
}
if (nbIntersect %2 != 0)
left.push_back(root->elements[i]);
else
right.push_back(root->elements[i]);
}
-
+
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() );
+ 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);
}
@@ -639,7 +644,7 @@ static void recur_cut_ (double R, double a1, double a2,
}
//--------------------------------------------------------------
-static void connected_left_right (std::vector<MElement *> &left,
+static void connected_left_right (std::vector<MElement *> &left,
std::vector<MElement *> &right ){
//connected left
@@ -668,7 +673,7 @@ static void printCut(std::map<MEdge,MVertex*,Less_Edge> &cutEdges, std::set<MEd
}
fprintf(f1,"};\n");
fclose(f1);
-
+
printf("Writing edges.pos \n");
std::set<MEdge,Less_Edge>::iterator itc = theCut.begin();
FILE *f2 = fopen("edges.pos","w");
@@ -687,7 +692,7 @@ static void printLevel(const char* fn,
double version,
double *dx = 0)
{
- if(!CTX::instance()->mesh.saveAll) return;
+ if(!CTX::instance()->mesh.saveAll) return;
std::set<MVertex*> vs;
for (unsigned int i = 0; i < elements.size(); i++)
@@ -696,9 +701,9 @@ static void printLevel(const char* fn,
bool binary = false;
FILE *fp = fopen (fn, "w");
- fprintf(fp, "$MeshFormat\n");
+ fprintf(fp, "$MeshFormat\n");
fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
- fprintf(fp, "$EndMeshFormat\n");
+ fprintf(fp, "$EndMeshFormat\n");
fprintf(fp, "$Nodes\n%d\n", (int)vs.size());
std::set<MVertex*> :: iterator it = vs.begin();
@@ -713,13 +718,13 @@ static void printLevel(const char* fn,
else fprintf(fp, "%d %g %g %g\n", (*it)->getIndex(),(*it)->x(), (*it)->y(), (*it)->z());
}
fprintf(fp, "$EndNodes\n");
-
+
fprintf(fp, "$Elements\n%d\n", (int)elements.size());
for (unsigned int i = 0; i < elements.size(); i++){
elements[i]->writeMSH(fp, version);
}
fprintf(fp, "$EndElements\n");
-
+
fclose(fp);
}
//--------------------------------------------------------------
@@ -733,27 +738,27 @@ static double localSize(MElement *e, std::map<MVertex*,SPoint2> &solution){
for(int j = 0; j<e->getNumVertices(); ++j){
SPoint2 p = solution[e->getVertex(j)];
local += SPoint3(p.x(),p.y(),0.0);
- }
+ }
return local.max().distance(local.min());
-// MVertex* v0 = e->getVertex(0);
-// MVertex* v1 = e->getVertex(1);
-// MVertex* v2 = e->getVertex(2);
-// double p0[3] = {v0->x(), v0->y(), v0->z()};
+// MVertex* v0 = e->getVertex(0);
+// MVertex* v1 = e->getVertex(1);
+// MVertex* v2 = e->getVertex(2);
+// double p0[3] = {v0->x(), v0->y(), v0->z()};
// double p1[3] = {v1->x(), v1->y(), v1->z()};
// double p2[3] = {v2->x(), v2->y(), v2->z()};
// double a_3D = fabs(triangle_area(p0, p1, p2));
// SPoint2 s1 = solution[v0];
// SPoint2 s2 = solution[v1];
// SPoint2 s3 = solution[v2];
-// double q0[3] = {s1.x(), s1.y(), 0.0};
+// double q0[3] = {s1.x(), s1.y(), 0.0};
// double q1[3] = {s2.x(), s2.y(), 0.0};
// double q2[3] = {s3.x(), s3.y(), 0.0};
-// double a_2D = fabs(triangle_area(q0, q1, q2));
-
+// double a_2D = fabs(triangle_area(q0, q1, q2));
+
// return a_2D; //a_2D / a_3D;
-
-
+
+
}
@@ -795,7 +800,7 @@ static void one2OneMap(std::vector<MElement *> &elements, std::map<MVertex*,SPoi
// 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;
@@ -807,20 +812,20 @@ static void one2OneMap(std::vector<MElement *> &elements, std::map<MVertex*,SPoi
// }
// }
// 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;
-
+
// }
// }
@@ -828,9 +833,9 @@ static void one2OneMap(std::vector<MElement *> &elements, std::map<MVertex*,SPoi
}
//--------------------------------------------------------------
-static bool checkOrientation(std::vector<MElement *> &elements,
+static bool checkOrientation(std::vector<MElement *> &elements,
std::map<MVertex*,SPoint2> &solution, int iter) {
-
+
double a_old = 0, a_new;
bool oriented = true;
int count = 0;
@@ -843,7 +848,7 @@ static bool checkOrientation(std::vector<MElement *> &elements,
double p1[2] = {v2[0],v2[1]};
double p2[2] = {v3[0],v3[1]};
a_new = robustPredicates::orient2d(p0, p1, p2);
- if(count == 0) a_old=a_new;
+ if(count == 0) a_old=a_new;
if(a_new*a_old < 0.){
oriented = false;
break;
@@ -852,8 +857,8 @@ static bool checkOrientation(std::vector<MElement *> &elements,
a_old = a_new;
}
count++;
- }
-
+ }
+
int iterMax = 1;
if(!oriented && iter < iterMax){
if (iter == 0) Msg::Debug("Parametrization is NOT 1 to 1 : applying cavity checks.");
@@ -870,14 +875,14 @@ static bool checkOrientation(std::vector<MElement *> &elements,
}
//--------------------------------------------------------------
-multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
- std::map<MVertex*, SPoint3> &allCoordinates)
+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;
+ //if (!CTX::instance()->mesh.smoothInternalEdges)return;
- //Find the boundary loops
+ //Find the boundary loops
//The loop with the largest equivalent radius is the Dirichlet boundary
std::vector<std::pair<MVertex*,double> > boundaryNodes;
ordering_dirichlet(elements,boundaryNodes);
@@ -898,26 +903,26 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
root->_name = "Root";
parametrize(*root);
-
+
//fill the coordinates
- std::vector<double> iScale;
+ std::vector<double> iScale;
std::vector<SPoint2> iCenter;
fillCoordinates(*root, allCoordinates, iScale, iCenter);
//Compute centers for the cut
int nbElems = 0;
recur_compute_centers_ (1.0, M_PI, 0.0, root, nbElems);
-
+
//Split the mesh in left and right
//or Cut the mesh in left and right
- splitElems(elements);
+ splitElems(elements);
//cutElems(elements);
}
-void multiscaleLaplace::fillCoordinates (multiscaleLaplaceLevel & level,
- std::map<MVertex*, SPoint3> &allCoordinates,
- std::vector<double> &iScale,
+void multiscaleLaplace::fillCoordinates (multiscaleLaplaceLevel & level,
+ std::map<MVertex*, SPoint3> &allCoordinates,
+ std::vector<double> &iScale,
std::vector<SPoint2> &iCenter){
iScale.push_back(level.scale);
@@ -935,13 +940,13 @@ void multiscaleLaplace::fillCoordinates (multiscaleLaplaceLevel & level,
}
}
-
+
for (unsigned int i=0;i<level.children.size();i++){
- multiscaleLaplaceLevel* m = level.children[i];
+ multiscaleLaplaceLevel* m = level.children[i];
fillCoordinates(*m, allCoordinates, iScale, iCenter);
}
-
+
}
void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
@@ -953,8 +958,8 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
for(int j = 0; j<e->getNumVertices(); ++j){
allNodes.insert(e->getVertex(j));
}
- }
-
+ }
+
//Parametrize level
std::map<MVertex*,SPoint2> solution;
parametrize_method(level, allNodes, solution, HARMONIC);
@@ -986,11 +991,11 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
//Only keep the connected elements vectors goodSize (the rest goes into tooSmall)
keepConnected(goodSize, tooSmall);
- //Add the not too small regions to the level.elements
+ //Add the not too small regions to the level.elements
std::vector<std::vector<MElement*> > regions_, regions ;
regions.clear(); regions_.clear();
connectedRegions (tooSmall,regions_);
- for (unsigned int i=0;i< regions_.size() ; i++){
+ for (unsigned int i=0;i< regions_.size() ; i++){
bool really_small_elements = false;
for (unsigned int k=0; k<regions_[i].size() ; k++){
MElement *e = regions_[i][k];
@@ -1003,14 +1008,14 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
}
else
goodSize.insert(goodSize.begin(), regions_[i].begin(), regions_[i].end() );
- }
+ }
//check for convex small regions patches
for (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){
+ 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));
@@ -1041,7 +1046,7 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
// for (int i=0;i< regions.size() ; i++){
// SPoint2 c = cents[i];
// double rad = rads[i];
- // for (std::vector<MElement*>::iterator it = regions[i].begin(); it != regions[i].end(); ++it){
+ // for (std::vector<MElement*>::iterator it = regions[i].begin(); it != regions[i].end(); ++it){
// MElement *e = *it;
// SPoint2 p0 = solution[e->getVertex(0)];
// SPoint2 p1 = solution[e->getVertex(1)];
@@ -1058,7 +1063,7 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
// }
// keepConnected(regions[i], goodSize);
//}
-
+
level.elements.clear();
level.elements = goodSize;
@@ -1083,16 +1088,16 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
//For every small region compute a new parametrization
Msg::Info("Level (%d-%d): %d connected small regions",level.recur, level.region, regions.size());
- for (unsigned int i = 0; i < regions.size(); i++){
+ for (unsigned int i = 0; i < regions.size(); i++){
std::set<MVertex*> tooSmallv;
tooSmallv.clear();
for (unsigned int k=0; k<regions[i].size() ; k++){
MElement *e = regions[i][k];
for(int j = 0; j<e->getNumVertices(); ++j){
tooSmallv.insert(e->getVertex(j));
- }
+ }
}
-
+
multiscaleLaplaceLevel *nextLevel = new multiscaleLaplaceLevel;
nextLevel->elements = regions[i];
nextLevel->recur = level.recur+1;
@@ -1108,7 +1113,7 @@ 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){
MVertex *v = *itv;
if (goodSizev.find(v) != goodSizev.end()){
@@ -1125,9 +1130,9 @@ void multiscaleLaplace::parametrize(multiscaleLaplaceLevel & level){
}
-void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
+void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
std::set<MVertex*> &allNodes,
- std::map<MVertex*,SPoint2> &solution,
+ std::map<MVertex*,SPoint2> &solution,
typeOfMapping tom)
{
@@ -1178,15 +1183,15 @@ void multiscaleLaplace::parametrize_method (multiscaleLaplaceLevel & level,
// solve
if (myAssembler.sizeOfR() != 0) _lsys->systemSolve();
- // get the values
+ // get the values
int count = 0;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value;
- myAssembler.getDofValue(v, 0, 1, value);
+ myAssembler.getDofValue(v, 0, 1, value);
if (step == 0)solution[v] = SPoint2(value,0);
else solution[v] = SPoint2(solution[v][0],value);
- }
+ }
_lsys->clear();
}
@@ -1202,10 +1207,10 @@ void multiscaleLaplace::cutElems(std::vector<MElement *> &elements)
std::set<MVertex*> cutVertices;
elements.clear();
- recur_cut_edges_ (root, cutEdges,cutVertices);
+ recur_cut_edges_ (root, cutEdges,cutVertices);
recur_cut_elements_ (root,cutEdges,cutVertices,theCut, 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++){
@@ -1216,13 +1221,13 @@ void multiscaleLaplace::cutElems(std::vector<MElement *> &elements)
leftS.clear();
std::vector<MElement*> left,right;
recur_leftCut_ (elements[0], e2e, theCut, leftS);
-
+
for (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() );
@@ -1233,9 +1238,9 @@ void multiscaleLaplace::cutElems(std::vector<MElement *> &elements)
elements.insert(elements.end(),left.begin(),left.end());
elements.insert(elements.end(),right.begin(),right.end());
- printLevel ("Rootcut-left.msh",left,0,2.2);
- printLevel ("Rootcut-right.msh",right,0,2.2);
- printLevel ("Rootcut-all.msh",elements, 0,2.2);
+ printLevel ("Rootcut-left.msh",left,0,2.2);
+ printLevel ("Rootcut-right.msh",right,0,2.2);
+ printLevel ("Rootcut-all.msh",elements, 0,2.2);
//exit(1);
}
void multiscaleLaplace::splitElems(std::vector<MElement *> &elements)
@@ -1246,9 +1251,9 @@ void multiscaleLaplace::splitElems(std::vector<MElement *> &elements)
recur_cut_ (1.0, M_PI, 0.0, root,left,right);
connected_left_right(left, right);
- printLevel ("Rootsplit-left.msh",left,0,2.2);
- printLevel ("Rootsplit-right.msh",right,0,2.2);
- printLevel ("Rootsplit-all.msh",elements, 0,2.2);
+ printLevel ("Rootsplit-left.msh",left,0,2.2);
+ printLevel ("Rootsplit-right.msh",right,0,2.2);
+ 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);
diff --git a/Solver/multiscaleLaplace.h b/Solver/multiscaleLaplace.h
index 62f57cc8b974d1845d08803bb5fa0a566f25a03f..b4be14fa30179c7c93a4529f0861642878bce317 100644
--- a/Solver/multiscaleLaplace.h
+++ b/Solver/multiscaleLaplace.h
@@ -1,3 +1,8 @@
+// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
#ifndef _MULTISCALE_LAPLACE_H_
#define _MULTISCALE_LAPLACE_H_
@@ -25,23 +30,23 @@ struct multiscaleLaplaceLevel {
class multiscaleLaplace{
public:
- multiscaleLaplace (std::vector<MElement *> &elements,
- std::map<MVertex*, SPoint3> &allCoordinates);
- void cutElems (std::vector<MElement *> &elements);
- void splitElems (std::vector<MElement *> &elements);
+ multiscaleLaplace (std::vector<MElement *> &elements,
+ std::map<MVertex*, SPoint3> &allCoordinates);
+ void cutElems (std::vector<MElement *> &elements);
+ void splitElems (std::vector<MElement *> &elements);
typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3} typeOfMapping;
multiscaleLaplaceLevel* root;
void fillCoordinates (multiscaleLaplaceLevel & level,
std::map<MVertex*, SPoint3> &allCoordinates,
- std::vector<double> &iScale,
+ std::vector<double> &iScale,
std::vector<SPoint2> &iCenter);
- void parametrize (multiscaleLaplaceLevel &);
- void parametrize_method (multiscaleLaplaceLevel & level,
+ void parametrize (multiscaleLaplaceLevel &);
+ void parametrize_method (multiscaleLaplaceLevel & level,
std::set<MVertex*> &allNodes,
- std::map<MVertex*,SPoint2> &solution,
+ std::map<MVertex*,SPoint2> &solution,
typeOfMapping tom);
-
+
};
#endif
diff --git a/contrib/Netgen/nglib_gmsh.cpp b/contrib/Netgen/nglib_gmsh.cpp
index cb1527981e094ee84dddd15b5cc053efd277e0b0..4af5637139e281b3068ba8690c2bb2b9a3df69e1 100644
--- a/contrib/Netgen/nglib_gmsh.cpp
+++ b/contrib/Netgen/nglib_gmsh.cpp
@@ -1,5 +1,5 @@
// Interface to the Netgen meshing kernel for Gmsh. This file replaces
-// the original nglib.cpp file from the Netgen distribution
+// the original nglib.cpp file from the Netgen distribution.
#include <GmshMessage.h>
#include <linalg.hpp>
@@ -24,11 +24,11 @@ namespace nglib
int index;
char txt[1024];
public:
- mystreambuf() : index(0) {}
+ mystreambuf() : index(0) {}
int sync()
- {
+ {
txt[index] = '\0';
- if(!index || (index == 1 && (txt[0] == '.' || txt[0] == '+' ||
+ if(!index || (index == 1 && (txt[0] == '.' || txt[0] == '+' ||
txt[0] == ' ' || txt[0] == '*'))){
// ignore these messages
}
@@ -38,8 +38,8 @@ namespace nglib
else
Msg::Info(txt);
}
- index = 0;
- return 0;
+ index = 0;
+ return 0;
}
int overflow(int ch)
{
@@ -53,7 +53,7 @@ namespace nglib
index++;
}
}
- return 0;
+ return 0;
}
};
@@ -85,11 +85,11 @@ namespace nglib
// Create a new netgen mesh object
Ng_Mesh * Ng_NewMesh ()
{
- Mesh * mesh = new Mesh;
+ Mesh * mesh = new Mesh;
mesh->AddFaceDescriptor (FaceDescriptor (1, 1, 0, 1));
return (Ng_Mesh*)(void*)mesh;
}
-
+
// Delete an existing netgen mesh object
void Ng_DeleteMesh (Ng_Mesh * mesh)
{
@@ -97,10 +97,10 @@ namespace nglib
{
// Delete the Mesh structures
((Mesh*)mesh)->DeleteMesh();
-
+
// Now delete the Mesh class itself
delete (Mesh*)mesh;
-
+
// Set the Ng_Mesh pointer to NULL
mesh = NULL;
}
@@ -112,7 +112,7 @@ namespace nglib
Mesh * m = (Mesh*)mesh;
m->AddPoint (Point3d (x[0], x[1], x[2]));
}
-
+
// Manually add a surface element of a given type to an existing mesh object
void Ng_AddSurfaceElement (Ng_Mesh * mesh, Ng_Surface_Element_Type et,
int * pi)
@@ -125,7 +125,7 @@ namespace nglib
el.PNum(3) = pi[2];
m->AddSurfaceElement (el);
}
-
+
// Manually add a volume element of a given type to an existing mesh object
void Ng_AddVolumeElement (Ng_Mesh * mesh, Ng_Volume_Element_Type et,
int * pi)
@@ -139,19 +139,19 @@ namespace nglib
el.PNum(4) = pi[3];
m->AddVolumeElement (el);
}
-
+
// Obtain the number of points in the mesh
int Ng_GetNP (Ng_Mesh * mesh)
{
return ((Mesh*)mesh) -> GetNP();
}
-
+
// Obtain the number of volume elements in the mesh
int Ng_GetNE (Ng_Mesh * mesh)
{
return ((Mesh*)mesh) -> GetNE();
}
-
+
// Return point coordinates of a given point index in the mesh
void Ng_GetPoint (Ng_Mesh * mesh, int num, double * x)
{
@@ -160,7 +160,7 @@ namespace nglib
x[1] = p.Y();
x[2] = p.Z();
}
-
+
// Return the volume element at a given index "pi"
Ng_Volume_Element_Type
Ng_GetVolumeElement (Ng_Mesh * mesh, int num, int * pi)
@@ -180,24 +180,24 @@ namespace nglib
}
return et;
}
-
+
// Generates volume mesh from an existing surface mesh
Ng_Result Ng_GenerateVolumeMesh (Ng_Mesh * mesh, Ng_Meshing_Parameters * mp)
{
Mesh * m = (Mesh*)mesh;
-
+
// Philippose - 30/08/2009
- // Do not locally re-define "mparam" here... "mparam" is a global
- // object
+ // Do not locally re-define "mparam" here... "mparam" is a global
+ // object
//MeshingParameters mparam;
mp->Transfer_Parameters();
-
+
m->CalcLocalH(mparam.grading);
-
+
MeshVolume (mparam, *m);
RemoveIllegalElements (*m);
OptimizeVolume (mparam, *m);
-
+
return NG_OK;
}
@@ -205,11 +205,11 @@ namespace nglib
Ng_Result Ng_GenerateVolumeMesh (Ng_Mesh * mesh, double maxh)
{
Mesh *m = (Mesh*)mesh;
-
+
MeshingParameters mparam;
mparam.uselocalh = 1;
mparam.maxh = maxh;
-
+
try{
m->CalcLocalH(mparam.grading);
MeshVolume(mparam, *m);
@@ -226,11 +226,11 @@ namespace nglib
Ng_Result Ng_OptimizeVolumeMesh(Ng_Mesh *mesh, double maxh)
{
Mesh *m = (Mesh*)mesh;
-
+
MeshingParameters mparam;
mparam.uselocalh = 1;
mparam.maxh = maxh;
-
+
try{
m->CalcLocalH(mparam.grading);
//MeshVolume(mparam, *m);
@@ -242,39 +242,39 @@ namespace nglib
}
return NG_OK;
}
-
+
// ------------------ Begin - Meshing Parameters related functions ------------------
// Constructor for the local nglib meshing parameters class
Ng_Meshing_Parameters :: Ng_Meshing_Parameters()
{
uselocalh = 1;
-
+
maxh = 1000;
minh = 0.0;
-
+
fineness = 0.5;
grading = 0.3;
-
+
elementsperedge = 2.0;
elementspercurve = 2.0;
-
+
closeedgeenable = 0;
closeedgefact = 2.0;
-
+
second_order = 0;
quad_dominated = 0;
-
+
meshsize_filename = 0;
-
+
optsurfmeshenable = 1;
optvolmeshenable = 1;
-
+
optsteps_2d = 3;
optsteps_3d = 3;
-
+
invert_tets = 0;
invert_trigs = 0;
-
+
check_overlap = 1;
check_overlapping_boundary = 1;
}
@@ -283,68 +283,68 @@ namespace nglib
void Ng_Meshing_Parameters :: Reset_Parameters()
{
uselocalh = 1;
-
+
maxh = 1000;
minh = 0;
-
+
fineness = 0.5;
grading = 0.3;
-
+
elementsperedge = 2.0;
elementspercurve = 2.0;
-
+
closeedgeenable = 0;
closeedgefact = 2.0;
-
+
second_order = 0;
quad_dominated = 0;
-
+
meshsize_filename = 0;
-
+
optsurfmeshenable = 1;
optvolmeshenable = 1;
-
+
optsteps_2d = 3;
optsteps_3d = 3;
-
+
invert_tets = 0;
invert_trigs = 0;
-
+
check_overlap = 1;
check_overlapping_boundary = 1;
}
- //
+ //
void Ng_Meshing_Parameters :: Transfer_Parameters()
{
mparam.uselocalh = uselocalh;
-
+
mparam.maxh = maxh;
mparam.minh = minh;
-
+
mparam.grading = grading;
mparam.curvaturesafety = elementspercurve;
mparam.segmentsperedge = elementsperedge;
-
+
mparam.secondorder = second_order;
mparam.quad = quad_dominated;
-
+
mparam.meshsizefilename = meshsize_filename;
-
+
mparam.optsteps2d = optsteps_2d;
mparam.optsteps3d = optsteps_3d;
-
+
mparam.inverttets = invert_tets;
mparam.inverttrigs = invert_trigs;
-
+
mparam.checkoverlap = check_overlap;
mparam.checkoverlappingboundary = check_overlapping_boundary;
}
-
+
} // End of namespace nglib
// compatibility functions:
-namespace netgen
+namespace netgen
{
char geomfilename[255];
@@ -391,7 +391,7 @@ namespace netgen
void Render()
{
- ;
+ ;
}
} // End of namespace netgen
diff --git a/utils/misc/find_missing_copyright.sh b/utils/misc/find_missing_copyright.sh
new file mode 100755
index 0000000000000000000000000000000000000000..023e9a57889b41b7711de6b136c84fea04d0a161
--- /dev/null
+++ b/utils/misc/find_missing_copyright.sh
@@ -0,0 +1,3 @@
+#!/bin/sh
+
+find . -type f ! -exec grep -q 'Copyright' {} \; -print