diff --git a/Geo/Chain.cpp b/Geo/Chain.cpp
index 2751451020452559d180c23014d18f06593e2c87..9aab03d5f447be30cb4f9705c40283788db11257 100644
--- a/Geo/Chain.cpp
+++ b/Geo/Chain.cpp
@@ -218,7 +218,7 @@ ElemChain ElemChain::getBoundaryElemChain(int i) const
bool ElemChain::inEntity(GEntity* e) const
{
if(_vertexCache[e].empty()) {
- for(int i = 0; i < e->getNumMeshElements(); i++)
+ for(unsigned int i = 0; i < e->getNumMeshElements(); i++)
for(int j = 0; j < e->getMeshElement(i)->getNumVertices(); j++)
_vertexCache[e].insert(e->getMeshElement(i)->getVertex(j));
}
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index e2a4b40e6bc3821dfc2d45a6e649e825b0d077a0..9547e3c0360cded10bbdbc8ff290ee54329a7298 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -285,7 +285,6 @@ double GEdge::length(const double &u0, const double &u1, const int nbQuadPoints)
use a golden section algorithm that minimizes
min_t \|x(t)-y\|
-
*/
const double GOLDEN = (1. + sqrt(5.)) / 2.;
@@ -356,7 +355,8 @@ GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const
return point(t);
}
-bool GEdge::computeDistanceFromMeshToGeometry (double &d2, double &dmax) {
+bool GEdge::computeDistanceFromMeshToGeometry (double &d2, double &dmax)
+{
d2 = 0.0; dmax = 0.0;
if (geomType() == Line) return true;
if (!lines.size())return false;
@@ -364,7 +364,7 @@ bool GEdge::computeDistanceFromMeshToGeometry (double &d2, double &dmax) {
int npts;
lines[0]->getIntegrationPoints(2*lines[0]->getPolynomialOrder(), &npts, &pts);
- for (int i=0;i<lines.size();i++){
+ for (unsigned int i = 0; i < lines.size(); i++){
MLine *l = lines[i];
double t[256];
@@ -386,11 +386,11 @@ bool GEdge::computeDistanceFromMeshToGeometry (double &d2, double &dmax) {
double tinit = l->interpolate(t,pts[j].pt[0],0,0);
GPoint pc = closestPoint(p, tinit);
if (!pc.succeeded())continue;
- double dsq =
- (pc.x()-p.x())*(pc.x()-p.x()) +
- (pc.y()-p.y())*(pc.y()-p.y()) +
+ double dsq =
+ (pc.x()-p.x())*(pc.x()-p.x()) +
+ (pc.y()-p.y())*(pc.y()-p.y()) +
(pc.z()-p.z())*(pc.z()-p.z());
- d2 += pts[i].weight * fabs(l->getJacobianDeterminant(pts[j].pt[0],0,0)) * dsq;
+ d2 += pts[i].weight * fabs(l->getJacobianDeterminant(pts[j].pt[0],0,0)) * dsq;
dmax = std::max(dmax,sqrt(dsq));
}
}
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 2f8e145200f2cd74173db8addc510c25e83e2f82..8ed0305eb7d825585f814ec5e4c2c07514f97be8 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -530,10 +530,12 @@ int GModel::mesh(int dimension)
#endif
}
-int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<double>* > f, std::vector<std::vector<double> > parameters, int niter, bool meshAll)
+int GModel::adaptMesh(std::vector<int> technique,
+ std::vector<simpleFunction<double>* > f,
+ std::vector<std::vector<double> > parameters,
+ int niter, bool meshAll)
{
#if defined(HAVE_MESH)
-
if (getNumMeshElements() == 0) mesh(getDim());
int nbElemsOld = getNumMeshElements();
int nbElems;
@@ -549,7 +551,7 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
fields->reset();
meshMetric *metric = new meshMetric(this);
- for (int imetric=0;imetric<technique.size();imetric++){
+ for (unsigned int imetric = 0; imetric < technique.size(); imetric++){
metric->addMetric(technique[imetric], f[imetric], parameters[imetric]);
}
fields->setBackgroundField(metric);
@@ -601,11 +603,11 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
}
}
- if (elements.size() == 0)return -1;
+ if (elements.size() == 0) return -1;
fields->reset();
meshMetric *metric = new meshMetric(this);
- for (int imetric=0;imetric<technique.size();imetric++){
+ for (unsigned int imetric = 0; imetric < technique.size(); imetric++){
metric->addMetric(technique[imetric], f[imetric], parameters[imetric]);
}
fields->setBackgroundField(metric);
@@ -1338,6 +1340,7 @@ static void recurConnectMElementsByMFace(const MFace &f,
printf("group pf %d elements found\n",(int)group.size());
}
+/*
static void recurConnectMElementsByMFaceOld(const MFace &f,
std::multimap<MFace, MElement*, Less_Face> &e2f,
std::set<MElement*> &group,
@@ -1354,6 +1357,7 @@ static void recurConnectMElementsByMFaceOld(const MFace &f,
}
}
}
+*/
static int connectedVolumes(std::vector<MElement*> &elements,
std::vector<std::vector<MElement*> > ®s)
@@ -2084,7 +2088,7 @@ GEdge* GModel::addCompoundEdge(std::vector<GEdge*> edges, int num){
}
else{
Curve *c = Create_Curve(num, MSH_SEGM_COMPOUND, 1, NULL, NULL, -1, -1, 0., 1.);
- for(int i= 0; i< edges.size(); i++)
+ for(unsigned int i= 0; i < edges.size(); i++)
c->compound.push_back(edges[i]->tag());
// Curve *c = Create_Curve(num, MSH_SEGM_DISCRETE, 1,
@@ -2140,7 +2144,7 @@ GFace* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int split,
}
else{
Surface *s = Create_Surface(num, MSH_SURF_COMPOUND);
- for(int i= 0; i< faces.size(); i++)
+ for(unsigned int i= 0; i < faces.size(); i++)
s->compound.push_back(faces[i]->tag());
std::list<GEdge*> edges = gfc->edges();
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 50dd7a5e99aa97917ad770e438b1a34fc31d440a..0cd2d3692c7474cf1be4396e7e1cfe8c0ddedef2 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -144,7 +144,7 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
int numf = gm->getMaxElementaryNumber(2)+1;
Surface *s = Create_Surface(numf, MSH_SURF_PLAN);
List_T *temp = List_Create(nLoops, nLoops, sizeof(int));
- for (unsigned int i=0; i< nLoops; i++){
+ for (int i = 0; i < nLoops; i++){
int numl = vecLoops[i]->Num;
List_Add(temp, &numl);
}
@@ -295,7 +295,7 @@ std::vector<GEntity*> GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, i
ep->mesh.NbLayer = 1; //this may be more general for defining different layers
ep->mesh.hLayer.clear();
ep->mesh.hLayer.push_back(hLayer);
- ep->mesh.NbElmLayer.clear();
+ ep->mesh.NbElmLayer.clear();
ep->mesh.NbElmLayer.push_back(nbLayers);
ep->mesh.ExtrudeMesh = true;
if (CTX::instance()->mesh.recombineAll){
@@ -304,8 +304,8 @@ std::vector<GEntity*> GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, i
}
else ep->mesh.Recombine = false;
ep->geo.Source = e->tag();
-
- int type = BOUNDARY_LAYER;
+
+ int type = BOUNDARY_LAYER;
double T0=0., T1=0., T2=0.;
double A0=0., A1=0., A2=0.;
double X0=0., X1=0., X2=0.,alpha=0.;
@@ -347,7 +347,7 @@ std::vector<GEntity*> GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, i
ep,
list_out);
- //create GEntities
+ //create GEntities
gm->importGEOInternals();
//return the new created entity
@@ -368,7 +368,7 @@ std::vector<GEntity*> GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, i
List_Read(list_out, j, &el);
GEdge *gel = gm->getEdgeByTag(el.Num);
extrudedEntities.push_back((GEntity*)gel);
- }
+ }
}
else if(e->dim()==2){
Shape s;
@@ -384,7 +384,7 @@ std::vector<GEntity*> GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, i
List_Read(list_out, j, &sl);
GFace *gfl = gm->getFaceByTag(sl.Num);
extrudedEntities.push_back((GEntity*)gfl);
- }
+ }
}
return extrudedEntities;
@@ -409,7 +409,7 @@ std::vector<GEntity*> GeoFactory::extrudeBoundaryLayer(GModel *gm, GEntity *e, i
// if(gf->getNativeType() == GEntity::GmshModel &&
// gf->geomType() == GEntity::BoundaryLayerSurface){
// ExtrudeParams *ep = gf->meshAttributes.extrude;
- // if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY
+ // if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY
// && std::abs(ep->geo.Source) == e->tag())
// newEnt = gf;
// }
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 0127d3ea354d5b17f398937c7380ea75dad04b8b..aa29cbba7af3e0411e7ff694dba739997a4076bb 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -172,6 +172,8 @@ static MElement *createElementMSH(GModel *m, int num, int typeMSH, int physical,
return e;
}
+#if (FAST_ELEMENTS==0)
+
static bool getElementsByNum(int elemNum[], std::map<int, std::vector<MElement*> > &elements,
bool erase, MElement *elems[], int nbElem = 1)
{
@@ -233,6 +235,8 @@ static void getDomains(int dom1Num, int dom2Num, int type,
}
}
+#endif
+
int GModel::readMSH(const std::string &name)
{
FILE *fp = fopen(name.c_str(), "rb");
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index c80e72785a95b69b722e690f81b7dc48e99f9dfa..e5c9f00d10ae4710f272bc6854ca88291f72bacc 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -439,7 +439,7 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
return V;
}
-// Interpolation transfinie sur un quadrangle :
+// Transfinite interpolation on a quadrangle :
// f(u,v) = (1-u)c4(v) + u c2(v) + (1-v)c1(u) + v c3(u)
// - [ (1-u)(1-v)s1 + u(1-v)s2 + uv s3 + (1-u)v s4 ]
@@ -461,14 +461,10 @@ static Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
s1.Pos.Y, s2.Pos.Y, s3.Pos.Y, s4.Pos.Y, u, v);
V.Pos.Z = TRAN_QUA(c1.Pos.Z, c2.Pos.Z, c3.Pos.Z, c4.Pos.Z,
s1.Pos.Z, s2.Pos.Z, s3.Pos.Z, s4.Pos.Z, u, v);
- //printf("TRANQUA %f %f %f %f %f %f %f %f\n", c1.Pos.Z, c2.Pos.Z, c3.Pos.Z, c4.Pos.Z,
- // s1.Pos.Z, s2.Pos.Z, s3.Pos.Z, s4.Pos.Z);
- //printf("V.Pos.Z %f\n", V.Pos.Z);
return (V);
}
-// Interpolation transfinie sur un triangle : TRAN_QUA avec s1=s4=c4
-// f(u,v) = u c2 (v) + (1-v) c1(u) + v c3(u) - u(1-v) s2 - uv s3
+// Transfinite interpolation on a triangle :
//
// s3(1,1)
// +
@@ -480,6 +476,12 @@ static Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
// +----------+
// s1(0,0) s2(1,0)
+#if 0
+
+// Old-style: TRAN_QUA with s1=s4=c4
+//
+// f(u,v) = u c2 (v) + (1-v) c1(u) + v c3(u) - u(1-v) s2 - uv s3
+//
// u = v = 0 -----> x = c1(0) = s1 --> OK
// u = 1 ; v = 0 -----> x = c2(0) + c1(1) - s2 = s2 --> OK
// u = 1 ; v = 1 -----> x = c2(1) + c3(1) - s3 = s3 --> OK
@@ -487,21 +489,8 @@ static Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
// u = 1 --> c2(v) + (1-v) s2 + v s3 -(1-v) s2 - v s3 --> x = c2(v) --> OK
// u = 0 --> (1-v) s1 + v s1 = s1 !!! not terrible (singular)
-// Transfinite approximation on a triangle
-
-// f(u,v) = (1-u) (c1(u-v) + c3(1-v) - s1) +
-// (u-v) (c2(v) + c1(u) - s2) +
-// v (c3(1-u) + c2(1-u+v) - s3)
-//
-// u = v = 0 --> f(0,0) = s1 + s1 - s1 = s1
-// u = v = 1 --> f(1,1) = s3 + s3 - s3 = s3
-// u = 1 ; v = 0 --> f(1,1) = s2 + s2 - s2 = s2
-// v = 0 --> (1-u)c1(u) + u c1(u) = c1(u) --> COOL
-
#define TRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) u*c2+(1.-v)*c1+v*c3-(u*(1.-v)*s2+u*v*s3);
-#define TRAN_TRIB(c1,c1b,c2,c2b,c3,c3b,s1,s2,s3,u,v) (1.-u)*(c1+c3b-s1)+(u-v)*(c2+c1b-s2)+v*(c3+c2b-s3);
-
static Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
Vertex s1, Vertex s2, Vertex s3, double u, double v)
{
@@ -517,6 +506,22 @@ static Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
return V;
}
+#endif
+
+// New-style:
+//
+// f(u,v) = (1-u) (c1(u-v) + c3(1-v) - s1) +
+// (u-v) (c2(v) + c1(u) - s2) +
+// v (c3(1-u) + c2(1-u+v) - s3)
+//
+// u = v = 0 --> f(0,0) = s1 + s1 - s1 = s1
+// u = v = 1 --> f(1,1) = s3 + s3 - s3 = s3
+// u = 1 ; v = 0 --> f(1,1) = s2 + s2 - s2 = s2
+// v = 0 --> (1-u)c1(u) + u c1(u) = c1(u) --> COOL
+
+#define TRAN_TRIB(c1,c1b,c2,c2b,c3,c3b,s1,s2,s3,u,v)\
+ (1.-u)*(c1+c3b-s1)+(u-v)*(c2+c1b-s2)+v*(c3+c2b-s3);
+
static Vertex TransfiniteTriB(Vertex c1, Vertex c1b, Vertex c2,
Vertex c2b, Vertex c3, Vertex c3b,
Vertex s1, Vertex s2, Vertex s3,
diff --git a/Geo/STensor3.cpp b/Geo/STensor3.cpp
index ae0cf3e6cf069baadfc49b240ef98a5e143a4fb0..4e7204aaba7ad4fea953b9e089ff9eea22927e48 100644
--- a/Geo/STensor3.cpp
+++ b/Geo/STensor3.cpp
@@ -8,17 +8,18 @@
void SMetric3::print (const char *s) const
{
- printf(" metric %s : %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E \n",s,
- (*this)(0,0),(*this)(1,1),(*this)(2,2),
- (*this)(0,1),(*this)(0,2),(*this)(1,2));
+ printf(" metric %s : %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E \n", s,
+ (*this)(0,0), (*this)(1,1), (*this)(2,2),
+ (*this)(0,1), (*this)(0,2), (*this)(1,2));
}
void STensor3::print (const char *s) const
{
- printf(" tensor %s : \n %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n",s,
- (*this)(0,0),(*this)(0,1),(*this)(0,2),
- (*this)(1,0),(*this)(1,1),(*this)(1,2),
- (*this)(2,0),(*this)(2,1),(*this)(2,2));
+ printf(" tensor %s : \n"
+ " %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n", s,
+ (*this)(0,0), (*this)(0,1), (*this)(0,2),
+ (*this)(1,0), (*this)(1,1), (*this)(1,2),
+ (*this)(2,0), (*this)(2,1), (*this)(2,2));
}
SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
@@ -61,20 +62,17 @@ SMetric3 intersection_conserve_mostaniso (const SMetric3 &m1, const SMetric3 &m2
fullVector<double> S1(3);
m1.eig(V1,S1,true);
double lambda1_min = std::min(std::min(fabs(S1(0)),fabs(S1(1))),fabs(S1(2)));
- double lambda1_max = std::max(std::max(fabs(S1(0)),fabs(S1(1))),fabs(S1(2)));
+ //double lambda1_max = std::max(std::max(fabs(S1(0)),fabs(S1(1))),fabs(S1(2)));
fullMatrix<double> V2(3,3);
fullVector<double> S2(3);
m2.eig(V2,S2,true);
double lambda2_min = std::min(std::min(fabs(S2(0)),fabs(S2(1))),fabs(S2(2)));
- double lambda2_max = std::max(std::max(fabs(S2(0)),fabs(S2(1))),fabs(S2(2)));
+ //double lambda2_max = std::max(std::max(fabs(S2(0)),fabs(S2(1))),fabs(S2(2)));
- double ratio1 = lambda1_min/lambda1_max;
- double ratio2 = lambda2_min/lambda2_max;
-
- if (lambda1_min<lambda2_min)
- return intersection_conserveM1(m1,m2);
+ if (lambda1_min < lambda2_min)
+ return intersection_conserveM1(m1, m2);
else
- return intersection_conserveM1(m2,m1);
+ return intersection_conserveM1(m2, m1);
}
// (1-t) * m1 + t * m2
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 4812288758b0a5da6e6ed1fd23693dd3315c965f..ec8ecff864434fd2c0678bf35224a16ff0dd223e 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1660,7 +1660,7 @@ class AttractorField : public Field
annDeallocPts(zeronodes);
delete kdtree;
}
-
+
std::vector<SPoint3> points;
std::vector<SPoint2> uvpoints;
std::vector<int> offset;
@@ -1674,17 +1674,19 @@ class AttractorField : public Field
SBoundingBox3d bb = f->bounds();
SVector3 dd = bb.max() - bb.min();
double maxDist = dd.norm() / n_nodes_by_edge ;
- bool success = f->fillPointCloud(maxDist, &points, &uvpoints);
+ f->fillPointCloud(maxDist, &points, &uvpoints);
offset.push_back(points.size());
}
}
}
- int totpoints =
- nodes_id.size() +
+ int totpoints =
+ nodes_id.size() +
(n_nodes_by_edge-2) * edges_id.size() +
- ((points.size()) ? points.size() : n_nodes_by_edge * n_nodes_by_edge * faces_id.size());
+ ((points.size()) ? points.size() :
+ n_nodes_by_edge * n_nodes_by_edge * faces_id.size());
+
printf("%d points found in points clouds (%d edges)\n",totpoints, edges_id.size());
if(totpoints){
@@ -1757,7 +1759,7 @@ class AttractorField : public Field
}
else {
GFace *f = GModel::current()->getFaceByTag(*it);
- if(f) {
+ if(f) {
if (points.size()){
for(int j = offset[count]; j < offset[count+1];j++) {
zeronodes[k][0] = points[j].x();
@@ -1783,7 +1785,7 @@ class AttractorField : public Field
}
}
}
- }
+ }
}
}
kdtree = new ANNkd_tree(zeronodes, totpoints, 3);
@@ -1847,7 +1849,7 @@ BoundaryLayerField::BoundaryLayerField()
double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge)
{
-
+
if (update_needed){
for(std::list<int>::iterator it = nodes_id.begin();
it != nodes_id.end(); ++it) {
@@ -1863,7 +1865,7 @@ double BoundaryLayerField::operator() (double x, double y, double z, GEntity *ge
}
update_needed = false;
}
-
+
double dist = 1.e22;
AttractorField *cc;
for (std::list<AttractorField*>::iterator it = _att_fields.begin();
@@ -1910,7 +1912,7 @@ void BoundaryLayerField::operator() (AttractorField *cc, double dist,
double beta = CTX::instance()->mesh.smoothRatio;
if (pp.first.dim ==0){
GVertex *v = GModel::current()->getVertexByTag(pp.first.ent);
- SVector3 t1;
+ SVector3 t1;
if (dist < thickness){
t1 = SVector3(1,0,0);
}
@@ -1952,7 +1954,7 @@ void BoundaryLayerField::operator() (AttractorField *cc, double dist,
double oneOverD2_min = .5/(b*b) * (1. + sqrt (1. + ( 4.*cmin*cmin*b*b*b*b/ (h*h*beta*beta))));
double oneOverD2_max = .5/(b*b) * (1. + sqrt (1. + ( 4.*cmax*cmax*b*b*b*b/ (h*h*beta*beta))));
double dmin = sqrt(1./oneOverD2_min);
- double dmax = sqrt(1./oneOverD2_max);
+ double dmax = sqrt(1./oneOverD2_max);
dmin = std::min(dmin,dmax*tgt_aniso_ratio);
metr = buildMetricTangentToSurface(dirMin,dirMax,dmin,dmax,lc_n);
return;
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index d55529ace458ce68c6708f4042fb7a0d32c76883..d23bc96d5c337fbb6b3c37597c8889195d80dc54 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -948,7 +948,7 @@ void signedDistancesPointsLine(std::vector<double> &distances,
distances.resize(pts.size());
closePts.clear();
closePts.resize(pts.size());
- for (int i=0; i<pts.size(); i++) {
+ for (unsigned int i = 0; i < pts.size(); i++) {
double d;
SPoint3 closePt;
const SPoint3 &p = pts[i];
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 75a030c7d26673fc90f0b733246fdc211ea00529..0e1062f0a152de37e4d43c14fba80d5a4ba4f24e 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -15,7 +15,9 @@
#include "GaussIntegration.h"
#include <limits>
-static void printClosure(polynomialBasis::clCont &fullClosure, std::vector<int> &closureRef,
+/*
+static void printClosure(polynomialBasis::clCont &fullClosure,
+ std::vector<int> &closureRef,
polynomialBasis::clCont &closures)
{
for (unsigned int i = 0; i < closures.size(); i++) {
@@ -38,6 +40,7 @@ static void printClosure(polynomialBasis::clCont &fullClosure, std::vector<int>
printf("\n");
}
}
+*/
int polynomialBasis::getTag(int parentTag, int order, bool serendip)
{
@@ -1240,7 +1243,8 @@ static void generateFaceClosureTetFull(polynomialBasis::clCont &closureFull,
}
}
-void rotateHex(int iFace, int iRot, int iSign, double uI, double vI, double &uO, double &vO, double &wO)
+void rotateHex(int iFace, int iRot, int iSign, double uI, double vI,
+ double &uO, double &vO, double &wO)
{
if (iSign<0){
double tmp=uI;
@@ -1310,23 +1314,27 @@ static void checkClosure(int tag) {
}
*/
-static void generateFaceClosureHex(polynomialBasis::clCont &closure, int order, bool serendip, const fullMatrix<double> &points)
+static void generateFaceClosureHex(polynomialBasis::clCont &closure, int order,
+ bool serendip, const fullMatrix<double> &points)
{
closure.clear();
- const polynomialBasis &fsFace = *polynomialBases::find(polynomialBasis::getTag(TYPE_QUA, order, serendip));
+ const polynomialBasis &fsFace = *polynomialBases::find
+ (polynomialBasis::getTag(TYPE_QUA, order, serendip));
for (int iRotate = 0; iRotate < 4; iRotate++){
for (int iSign = 1; iSign >= -1; iSign -= 2){
for (int iFace = 0; iFace < 6; iFace++) {
polynomialBasis::closure cl;
cl.type = fsFace.type;
cl.resize(fsFace.points.size1());
- for (int iNode = 0; iNode < cl.size(); ++iNode) {
+ for (unsigned int iNode = 0; iNode < cl.size(); ++iNode) {
double u,v,w;
- rotateHex(iFace, iRotate, iSign, fsFace.points(iNode, 0), fsFace.points(iNode, 1), u, v, w);
+ rotateHex(iFace, iRotate, iSign, fsFace.points(iNode, 0),
+ fsFace.points(iNode, 1), u, v, w);
cl[iNode] = 0;
double D = std::numeric_limits<double>::max();
for (int jNode = 0; jNode < points.size1(); ++jNode) {
- double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) + pow2(points(jNode, 2) - w);
+ double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) +
+ pow2(points(jNode, 2) - w);
if (d < D) {
cl[iNode] = jNode;
D = d;
@@ -1341,7 +1349,8 @@ static void generateFaceClosureHex(polynomialBasis::clCont &closure, int order,
static void generateFaceClosureHexFull(polynomialBasis::clCont &closure,
std::vector<int> &closureRef,
- int order, bool serendip, const fullMatrix<double> &points)
+ int order, bool serendip,
+ const fullMatrix<double> &points)
{
closure.clear();
int clId = 0;
@@ -1352,11 +1361,13 @@ static void generateFaceClosureHexFull(polynomialBasis::clCont &closure,
cl.resize(points.size1());
for (int iNode = 0; iNode < points.size1(); ++iNode) {
double u,v,w;
- rotateHexFull(iFace, iRotate, iSign, points(iNode, 0), points(iNode, 1), points(iNode, 2), u, v, w);
+ rotateHexFull(iFace, iRotate, iSign, points(iNode, 0), points(iNode, 1),
+ points(iNode, 2), u, v, w);
int J = 0;
double D = std::numeric_limits<double>::max();
for (int jNode = 0; jNode < points.size1(); ++jNode) {
- double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) + pow2(points(jNode, 2) - w);
+ double d = pow2(points(jNode, 0) - u) + pow2(points(jNode, 1) - v) +
+ pow2(points(jNode, 2) - w);
if (d < D) {
J = jNode;
D = d;
@@ -1682,7 +1693,7 @@ const polynomialBasis *polynomialBases::find(int tag)
generateClosureOrder0(F.closures,24);
generateClosureOrder0(F.fullClosures, 24);
F.closureRef.resize(24, 0);
- }
+ }
else {
generateFaceClosureTet(F.closures, F.order);
generateFaceClosureTetFull(F.fullClosures, F.closureRef, F.order, F.serendip);