From 07c95f33d63e18754fac5952aaea8d5cf8608d44 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-Fran=C3=A7ois=20Remacle=20=28students=29?=
<jean-francois.remacle@uclouvain.be>
Date: Mon, 19 Oct 2009 16:58:44 +0000
Subject: [PATCH] clean up + fix bug in recurElement (gbricteux)
---
Solver/SElement.cpp | 62 ++--
Solver/SElement.h | 36 +--
Solver/elasticitySolver.cpp | 97 ++++---
Solver/elasticitySolver.h | 14 +-
Solver/elasticityTerm.cpp | 87 +++---
Solver/elasticityTerm.h | 18 +-
Solver/femTerm.h | 36 +--
Solver/groupOfElements.cpp | 10 +-
Solver/groupOfElements.h | 38 +--
Solver/linearSystemCSR.cpp | 268 +++++++++---------
contrib/DiscreteIntegration/Integration3D.cpp | 178 +++---------
contrib/DiscreteIntegration/recurCut.cpp | 28 +-
utils/api_demos/Makefile | 4 +-
utils/api_demos/mainElasticity.cpp | 3 +-
14 files changed, 393 insertions(+), 486 deletions(-)
diff --git a/Solver/SElement.cpp b/Solver/SElement.cpp
index 5be09e28ca..8c9edcfd34 100644
--- a/Solver/SElement.cpp
+++ b/Solver/SElement.cpp
@@ -17,35 +17,35 @@ etc.
*/
-simpleFunction<double> * SElement::_enrichement_s = 0,* SElement::_enrichement_t = 0;
+simpleFunction<double> *SElement::_enrichement_s = 0, *SElement::_enrichement_t = 0;
-void SElement::gradNodalFunctions ( double u, double v, double w, double invjac[3][3],double Grads[][3],
- simpleFunction<double> *_enrichement)
+void SElement::gradNodalFunctions (double u, double v, double w, double invjac[3][3], double Grads[][3],
+ simpleFunction<double> *_enrichement)
{
double grads[256][3];
_e->getGradShapeFunctions(u, v, w, grads);
int nbNodes = getNumNodalShapeFunctions();
for (int j = 0; j < nbNodes; j++){
- Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
+ Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
invjac[0][2] * grads[j][2];
- Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
+ Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
invjac[1][2] * grads[j][2];
- Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
+ Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
invjac[2][2] * grads[j][2];
}
-
+
if (_enrichement){
const int N = getNumNodalShapeFunctions();
SPoint3 p;
double sf[256];
_e->getShapeFunctions(u, v, w, sf);
// FIXME : re-use sf for computing coordinates
- _e->pnt(u,v,w,p);
- double E = (*_enrichement_s)(p.x(),p.y(),p.z());
- double dEdx,dEdy,dEdz;
- _enrichement_s->gradient(p.x(),p.y(),p.z(),dEdx,dEdy,dEdz);
- for (int i=0;i<N;i++){
+ _e->pnt(u, v, w, p);
+ double E = (*_enrichement_s)(p.x(), p.y(), p.z());
+ double dEdx, dEdy, dEdz;
+ _enrichement_s->gradient(p.x(), p.y(), p.z(), dEdx, dEdy, dEdz);
+ for (int i = 0; i < N; i++){
Grads[i][0] = Grads[i][0] * E + dEdx * sf[i];
Grads[i][1] = Grads[i][1] * E + dEdy * sf[i];
Grads[i][2] = Grads[i][2] * E + dEdz * sf[i];
@@ -53,51 +53,51 @@ void SElement::gradNodalFunctions ( double u, double v, double w, double invjac[
}
}
-void SElement::nodalFunctions ( double u, double v, double w, double s[],
- simpleFunction<double> *_enrichement)
+void SElement::nodalFunctions (double u, double v, double w, double s[],
+ simpleFunction<double> *_enrichement)
{
_e->getShapeFunctions(u, v, w, s);
if (_enrichement){
const int N = getNumNodalShapeFunctions();
SPoint3 p;
// FIXME : re-use s for computing coordinates
- _e->pnt(u,v,w,p);
- double E = (*_enrichement_s)(p.x(),p.y(),p.z());
- for (int i=0;i<N;i++){
+ _e->pnt(u, v, w, p);
+ double E = (*_enrichement_s)(p.x(), p.y(), p.z());
+ for (int i = 0; i < N; i++){
s[i] *= E;
}
}
}
-void SElement::gradNodalShapeFunctions ( double u, double v, double w, double invjac[3][3],
- double grads[][3])
+void SElement::gradNodalShapeFunctions (double u, double v, double w, double invjac[3][3],
+ double grads[][3])
{
- gradNodalFunctions (u,v,w,invjac,grads,_enrichement_s);
+ gradNodalFunctions (u, v, w, invjac, grads, _enrichement_s);
}
-void SElement::gradNodalTestFunctions ( double u, double v, double w, double invjac[3][3],
- double grads[][3])
+void SElement::gradNodalTestFunctions (double u, double v, double w, double invjac[3][3],
+ double grads[][3])
{
- gradNodalFunctions (u,v,w,invjac,grads,_enrichement_t);
+ gradNodalFunctions (u, v, w, invjac, grads, _enrichement_t);
}
-void SElement::nodalShapeFunctions ( double u, double v, double w, double s[])
+void SElement::nodalShapeFunctions (double u, double v, double w, double s[])
{
- nodalFunctions (u,v,w,s,_enrichement_s);
+ nodalFunctions (u, v, w, s, _enrichement_s);
}
-void SElement::nodalTestFunctions ( double u, double v, double w, double s[])
+void SElement::nodalTestFunctions (double u, double v, double w, double s[])
{
- nodalFunctions (u,v,w,s,_enrichement_t);
+ nodalFunctions (u, v, w, s, _enrichement_t);
}
-int SElement::getNumNodalShapeFunctions ( ) const {
- if (_e->getParent())return _e->getParent()->getNumVertices();
+int SElement::getNumNodalShapeFunctions () const {
+ if (_e->getParent()) return _e->getParent()->getNumVertices();
return _e->getNumVertices();
}
-int SElement::getNumNodalTestFunctions ( ) const {
- if (_e->getParent())return _e->getParent()->getNumVertices();
+int SElement::getNumNodalTestFunctions () const {
+ if (_e->getParent()) return _e->getParent()->getNumVertices();
return _e->getNumVertices();
}
diff --git a/Solver/SElement.h b/Solver/SElement.h
index d921d4ff5e..49ac68d718 100644
--- a/Solver/SElement.h
+++ b/Solver/SElement.h
@@ -22,29 +22,29 @@ class SElement
MElement *_e;
// store discrete function space and other data here
// ...
- static simpleFunction<double> *_enrichement_s,*_enrichement_t;
+ static simpleFunction<double> *_enrichement_s, *_enrichement_t;
// gradient of functions (possibly enriched)
- void nodalFunctions ( double u, double v, double w, double s[],
- simpleFunction<double> *_enrichement);
- void gradNodalFunctions ( double u, double v, double w, double invjac[3][3],
- double grad[][3], simpleFunction<double> *_enrichment);
+ void nodalFunctions (double u, double v, double w, double s[],
+ simpleFunction<double> *_enrichement);
+ void gradNodalFunctions (double u, double v, double w, double invjac[3][3],
+ double grad[][3], simpleFunction<double> *_enrichment);
public:
- SElement(MElement *e) : _e(e) {}
+ SElement(MElement *e) : _e(e) {}
~SElement(){}
- MElement *getMeshElement() const { return _e; }
- static void setShapeEnrichement(simpleFunction<double>*f){_enrichement_s=f;}
- static void setTestEnrichement(simpleFunction<double>*f){_enrichement_t=f;}
- static const simpleFunction<double>* getShapeEnrichement(){return _enrichement_s;}
- static const simpleFunction<double>* getTestEnrichement(){return _enrichement_t;}
+ MElement *getMeshElement() const { return _e; }
+ static void setShapeEnrichement(simpleFunction<double>*f) {_enrichement_s = f;}
+ static void setTestEnrichement(simpleFunction<double>*f) {_enrichement_t = f;}
+ static const simpleFunction<double>* getShapeEnrichement() {return _enrichement_s;}
+ static const simpleFunction<double>* getTestEnrichement() {return _enrichement_t;}
int getNumNodalShapeFunctions() const;
- inline MVertex *getVertex(int i) const {return _e->getParent() ? _e->getParent()->getVertex(i) : _e->getVertex(i) ; }
+ inline MVertex *getVertex(int i) const {return _e->getParent() ? _e->getParent()->getVertex(i) : _e->getVertex(i);}
int getNumNodalTestFunctions() const;
- void nodalShapeFunctions ( double u, double v, double w, double s[]);
- void gradNodalShapeFunctions ( double u, double v, double w, double invjac[3][3],
- double grad[][3]);
- void nodalTestFunctions ( double u, double v, double w, double s[]);
- void gradNodalTestFunctions ( double u, double v, double w, double invjac[3][3],
- double grad[][3]);
+ void nodalShapeFunctions (double u, double v, double w, double s[]);
+ void gradNodalShapeFunctions (double u, double v, double w, double invjac[3][3],
+ double grad[][3]);
+ void nodalTestFunctions (double u, double v, double w, double s[]);
+ void gradNodalTestFunctions (double u, double v, double w, double invjac[3][3],
+ double grad[][3]);
};
#endif
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 5b257e316a..0301f64de8 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -18,29 +18,30 @@
#include "PView.h"
#include "PViewData.h"
-PView* elasticitySolver::buildDisplacementView (const std::string &postFileName){
+PView* elasticitySolver::buildDisplacementView (const std::string &postFileName){
std::map<int, std::vector<double> > data;
std::set<MVertex*> v;
- for (int i=0;i<elasticFields.size();++i){
+ for (int i = 0; i < elasticFields.size(); ++i){
groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
- for ( ; it != elasticFields[i].g->end() ; ++it ){
+ for ( ; it != elasticFields[i].g->end(); ++it){
SElement se(*it);
- for (int j=0;j<se.getNumNodalShapeFunctions();++j){
- v.insert(se.getVertex(j));
+ for (int j = 0; j < se.getNumNodalShapeFunctions(); ++j){
+ v.insert(se.getVertex(j));
}
}
}
std::set<MVertex*>::iterator it = v.begin();
- for ( ; it!=v.end();++it){
+ for ( ; it != v.end(); ++it){
std::vector<double> d;
- d.push_back(0);d.push_back(0);d.push_back(0);
- for (int i=0;i<elasticFields.size();++i){
- const double E = (elasticFields[i]._enrichment) ? (*elasticFields[i]._enrichment)((*it)->x(),(*it)->y(),(*it)->z()) : 1.0;
+ d.push_back(0); d.push_back(0); d.push_back(0);
+ for (int i = 0; i < elasticFields.size(); ++i){
+ const double E = (elasticFields[i]._enrichment) ?
+ (*elasticFields[i]._enrichment)((*it)->x(), (*it)->y(), (*it)->z()) : 1.0;
// printf("%g %d \n",pAssembler->getDofValue(*it,0,elasticFields[i]._tag),elasticFields[i]._tag);
- d[0] += pAssembler->getDofValue(*it,0,elasticFields[i]._tag) * E;
- d[1] += pAssembler->getDofValue(*it,1,elasticFields[i]._tag) * E;
- d[2] += pAssembler->getDofValue(*it,2,elasticFields[i]._tag) * E;
+ d[0] += pAssembler->getDofValue(*it, 0, elasticFields[i]._tag) * E;
+ d[1] += pAssembler->getDofValue(*it, 1, elasticFields[i]._tag) * E;
+ d[2] += pAssembler->getDofValue(*it, 2, elasticFields[i]._tag) * E;
}
data[(*it)->getNum()] = d;
}
@@ -76,7 +77,7 @@ static double vonMises(dofManager<double,double> *a, MElement *e,
e->interpolateGrad(valx, u, v, w, gradux);
e->interpolateGrad(valy, u, v, w, graduy);
e->interpolateGrad(valz, u, v, w, graduz);
-
+
double eps[6] = {gradux[0], graduy[1], graduz[2],
0.5 * (gradux[1] + graduy[0]),
0.5 * (gradux[2] + graduz[0]),
@@ -94,15 +95,15 @@ static double vonMises(dofManager<double,double> *a, MElement *e,
double s[9] = {sxx, sxy, sxz,
sxy, syy, syz,
sxz, syz, szz};
-
+
return ComputeVonMises(s);
}
-
+
void elasticitySolver::setMesh(const std::string &meshFileName)
{
pModel = new GModel();
pModel->readMSH(meshFileName.c_str());
- _dim = pModel->getNumRegions() ? 3 : 2;
+ _dim = pModel->getNumRegions() ? 3 : 2;
}
void elasticitySolver::readInputFile(const std::string &fn)
@@ -112,7 +113,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
while(!feof(f)){
if(fscanf(f, "%s", what) != 1) return;
// printf("%s\n", what);
- if (!strcmp(what,"ElasticDomain")){
+ if (!strcmp(what, "ElasticDomain")){
elasticField field;
int physical;
if(fscanf(f, "%d %lf %lf", &physical, &field._E, &field._nu) != 3) return;
@@ -121,7 +122,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
field.g = new groupOfElements (_dim, physical);
elasticFields.push_back(field);
}
- else if (!strcmp(what,"Void")){
+ else if (!strcmp(what, "Void")){
// elasticField field;
// create enrichment ...
// create the group ...
@@ -176,10 +177,14 @@ void elasticitySolver::readInputFile(const std::string &fn)
if(fscanf(f, "%s", name) != 1) return;
setMesh(name);
}
+ else {
+ Msg::Error("Invalid input : %s", what);
+ return;
+ }
}
fclose(f);
}
-
+
void elasticitySolver::solve()
{
#if defined(HAVE_TAUCS)
@@ -191,54 +196,54 @@ void elasticitySolver::solve()
lsys->setNoisy(2);
#endif
pAssembler = new dofManager<double, double>(lsys);
-
+
std::map<int, std::vector<GEntity*> > groups[4];
pModel->getPhysicalGroups(groups);
-
+
// we first do all fixations. the behavior of the dofManager is to
// give priority to fixations : when a dof is fixed, it cannot be
// numbered afterwards
-
+
for (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
it != nodalDisplacements.end(); ++it){
MVertex *v = pModel->getMeshVertexByTag(it->first.first);
- pAssembler->fixVertex(v , it->first.second, _tag, it->second);
- printf("-- Fixing node %3d comp %1d to %8.5f\n",
+ pAssembler->fixVertex(v, it->first.second, _tag, it->second);
+ printf("-- Fixing node %3d comp %1d to %8.5f\n",
it->first.first, it->first.second, it->second);
}
-
+
for (std::map<std::pair<int, int>, double>::iterator it = edgeDisplacements.begin();
it != edgeDisplacements.end(); ++it){
elasticityTerm El(pModel, 1, 1, _tag);
- El.dirichletNodalBC(it->first.first, 1, it->first.second, _tag,
+ El.dirichletNodalBC(it->first.first, 1, it->first.second, _tag,
simpleFunction<double>(it->second), *pAssembler);
- printf("-- Fixing edge %3d comp %1d to %8.5f\n",
+ printf("-- Fixing edge %3d comp %1d to %8.5f\n",
it->first.first, it->first.second, it->second);
}
-
+
for (std::map<std::pair<int, int>, double>::iterator it = faceDisplacements.begin();
it != faceDisplacements.end(); ++it){
elasticityTerm El(pModel, 1, 1, _tag);
- El.dirichletNodalBC(it->first.first, 2, it->first.second, _tag,
+ El.dirichletNodalBC(it->first.first, 2, it->first.second, _tag,
simpleFunction<double>(it->second), *pAssembler);
- printf("-- Fixing face %3d comp %1d to %8.5f\n",
+ printf("-- Fixing face %3d comp %1d to %8.5f\n",
it->first.first, it->first.second, it->second);
}
-
+
// we number the nodes : when a node is numbered, it cannot be numbered
// again with another number. so we loop over elements
- for (int i=0;i<elasticFields.size();++i){
+ for (int i = 0; i < elasticFields.size(); ++i){
groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
- for ( ; it != elasticFields[i].g->end() ; ++it ){
+ for ( ; it != elasticFields[i].g->end(); ++it){
SElement se(*it);
- for (int j=0;j<se.getNumNodalShapeFunctions();++j){
- pAssembler->numberVertex(se.getVertex(j), 0, elasticFields[i]._tag);
- pAssembler->numberVertex(se.getVertex(j), 1, elasticFields[i]._tag);
- pAssembler->numberVertex(se.getVertex(j), 2, elasticFields[i]._tag);
+ for (int j = 0; j < se.getNumNodalShapeFunctions(); ++j){
+ pAssembler->numberVertex(se.getVertex(j), 0, elasticFields[i]._tag);
+ pAssembler->numberVertex(se.getVertex(j), 1, elasticFields[i]._tag);
+ pAssembler->numberVertex(se.getVertex(j), 2, elasticFields[i]._tag);
}
}
}
-
+
// Now we start the assembly process
// First build the force vector
// Nodes at geometrical vertices
@@ -247,11 +252,11 @@ void elasticitySolver::solve()
int iVertex = it->first;
SVector3 f = it->second;
std::vector<GEntity*> ent = groups[0][iVertex];
- for (unsigned int i = 0; i < ent.size(); i++){
+ for (unsigned int i = 0; i < ent.size(); i++){
pAssembler->assemble(ent[i]->mesh_vertices[0], 0, _tag, f.x());
pAssembler->assemble(ent[i]->mesh_vertices[0], 1, _tag, f.y());
pAssembler->assemble(ent[i]->mesh_vertices[0], 2, _tag, f.z());
- printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n",
+ printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n",
ent[i]->mesh_vertices[0]->getNum(), iVertex, f.x(), f.y(), f.z());
}
}
@@ -279,7 +284,7 @@ void elasticitySolver::solve()
El.neumannNodalBC(iFace, 2, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n", iFace, f.x(), f.y(), f.z());
}
-
+
// volume forces
for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
it != volumeForces.end(); ++it){
@@ -294,16 +299,16 @@ void elasticitySolver::solve()
// El.addToRightHandSide(*pAssembler, ent[i]);
}
}
-
+
// assembling the stifness matrix
- for (int i=0;i<elasticFields.size();i++){
+ for (int i = 0; i < elasticFields.size(); i++){
SElement::setShapeEnrichement(elasticFields[i]._enrichment);
- for (int j=0;j<elasticFields.size();j++){
+ for (int j = 0; j < elasticFields.size(); j++){
elasticityTerm El(0, elasticFields[i]._E, elasticFields[i]._nu,
- elasticFields[i]._tag,elasticFields[j]._tag);
+ elasticFields[i]._tag, elasticFields[j]._tag);
SElement::setTestEnrichement(elasticFields[j]._enrichment);
// printf("%d %d\n",elasticFields[i]._tag,elasticFields[j]._tag);
- El.addToMatrix(*pAssembler,*elasticFields[i].g, *elasticFields[j].g);
+ El.addToMatrix(*pAssembler, *elasticFields[i].g, *elasticFields[j].g);
}
}
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index eba195515e..10fa60f342 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -20,8 +20,8 @@ struct elasticField {
int _tag; // tag for the dofManager
groupOfElements *g; // support for this field
simpleFunction<double> *_enrichment; // XFEM enrichment
- double _E,_nu; // specific elastic datas (should be somewhere else)
- elasticField () : g(0),_enrichment(0){}
+ double _E, _nu; // specific elastic datas (should be somewhere else)
+ elasticField () : g(0), _enrichment(0){}
};
// an elastic solver ...
@@ -31,7 +31,7 @@ class elasticitySolver{
int _dim, _tag;
dofManager<double, double> *pAssembler;
// young modulus and poisson coefficient per physical
- std::vector<elasticField > elasticFields;
+ std::vector<elasticField> elasticFields;
// imposed nodal forces
std::map<int, SVector3> nodalForces;
// imposed line forces
@@ -49,10 +49,10 @@ class elasticitySolver{
public:
elasticitySolver(int tag) : _tag(tag) {}
void addNodalForces (int iNode, const SVector3 &f)
- {
+ {
nodalForces[iNode] = f;
}
- void addNodalDisplacement(int iNode, int dir, double val)
+ void addNodalDisplacement(int iNode, int dir, double val)
{
nodalDisplacements[std::make_pair(iNode, dir)] = val;
}
@@ -60,8 +60,8 @@ class elasticitySolver{
// {
// elasticConstants[physical] = std::make_pair(e, nu);
// }
- void setMesh(const std::string &meshFileName);
- virtual void solve();
+ void setMesh(const std::string &meshFileName);
+ virtual void solve();
void readInputFile(const std::string &meshFileName);
PView *buildDisplacementView(const std::string &postFileName);
// PView *buildVonMisesView(const std::string &postFileName);
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index 5b42a5227a..62351800c0 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -21,21 +21,21 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
double Grads[100][3];
double GradsT[100][3];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
-
+
m.set_all(0.);
-
+
double FACT = _E / (1 + _nu);
double C11 = FACT * (1 - _nu) / (1 - 2 * _nu);
double C12 = FACT * _nu / (1 - 2 * _nu);
double C44 = (C11 - C12) / 2;
const double C[6][6] =
- { {C11, C12, C12, 0, 0, 0},
- {C12, C11, C12, 0, 0, 0},
- {C12, C12, C11, 0, 0, 0},
+ { {C11, C12, C12, 0, 0, 0},
+ {C12, C11, C12, 0, 0, 0},
+ {C12, C12, C11, 0, 0, 0},
{ 0, 0, 0, C44, 0, 0},
- { 0, 0, 0, 0, C44, 0},
+ { 0, 0, 0, 0, C44, 0},
{ 0, 0, 0, 0, 0, C44} };
-
+
fullMatrix<double> H(6, 6);
fullMatrix<double> B(6, 3 * nbNodes);
fullMatrix<double> BTH(3 * nbNodes, 6);
@@ -43,7 +43,7 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
for (int i = 0; i < 6; i++)
for (int j = 0; j < 6; j++)
H(i, j) = C[i][j];
-
+
for (int i = 0; i < npts; i++){
const double u = GP[i].pt[0];
const double v = GP[i].pt[1];
@@ -55,57 +55,56 @@ void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
B.set_all(0.);
BT.set_all(0.);
-
+
if (se->getShapeEnrichement() == se->getTestEnrichement()){
for (int j = 0; j < nbNodes; j++){
- // printf(" GR(j) = %12.5E,%12.5E,%12.5E\n", Grads[j][0],Grads[j][1],Grads[j][2]);
-
- BT(j, 0) = B(0, j) = Grads[j][0];
- BT(j, 3) = B(3, j) = Grads[j][1];
- BT(j, 4) = B(4, j) = Grads[j][2];
-
- BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
- BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
- BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
-
- BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
- BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
- BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
+ //printf(" GR(j) = %12.5E,%12.5E,%12.5E\n", Grads[j][0],Grads[j][1],Grads[j][2]);
+
+ BT(j, 0) = B(0, j) = Grads[j][0];
+ BT(j, 3) = B(3, j) = Grads[j][1];
+ BT(j, 4) = B(4, j) = Grads[j][2];
+
+ BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
+ BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
+ BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
+
+ BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
+ BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
+ BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
}
}
else{
printf("coucou AAAAAAAAAAAAARGH \n");
se->gradNodalTestFunctions (u, v, w, invjac,GradsT);
for (int j = 0; j < nbNodes; j++){
- BT(j, 0) = Grads[j][0]; B(0, j) = GradsT[j][0];
- BT(j, 3) = Grads[j][1]; B(3, j) = GradsT[j][1];
- BT(j, 4) = Grads[j][2]; B(4, j) = GradsT[j][2];
-
- BT(j + nbNodes, 1) = Grads[j][1]; B(1, j + nbNodes) = GradsT[j][1];
- BT(j + nbNodes, 3) = Grads[j][0]; B(3, j + nbNodes) = GradsT[j][0];
- BT(j + nbNodes, 5) = Grads[j][2]; B(5, j + nbNodes) = GradsT[j][2];
-
- BT(j + 2 * nbNodes, 2) = Grads[j][2]; B(2, j + 2 * nbNodes) = GradsT[j][2];
- BT(j + 2 * nbNodes, 4) = Grads[j][0]; B(4, j + 2 * nbNodes) = GradsT[j][0];
- BT(j + 2 * nbNodes, 5) = Grads[j][1]; B(5, j + 2 * nbNodes) = GradsT[j][1];
+ BT(j, 0) = Grads[j][0]; B(0, j) = GradsT[j][0];
+ BT(j, 3) = Grads[j][1]; B(3, j) = GradsT[j][1];
+ BT(j, 4) = Grads[j][2]; B(4, j) = GradsT[j][2];
+
+ BT(j + nbNodes, 1) = Grads[j][1]; B(1, j + nbNodes) = GradsT[j][1];
+ BT(j + nbNodes, 3) = Grads[j][0]; B(3, j + nbNodes) = GradsT[j][0];
+ BT(j + nbNodes, 5) = Grads[j][2]; B(5, j + nbNodes) = GradsT[j][2];
+
+ BT(j + 2 * nbNodes, 2) = Grads[j][2]; B(2, j + 2 * nbNodes) = GradsT[j][2];
+ BT(j + 2 * nbNodes, 4) = Grads[j][0]; B(4, j + 2 * nbNodes) = GradsT[j][0];
+ BT(j + 2 * nbNodes, 5) = Grads[j][1]; B(5, j + 2 * nbNodes) = GradsT[j][1];
}
}
BTH.set_all(0.);
- BTH.gemm(BT, H);
+ BTH.gemm(BT, H);
m.gemm(BTH, B, weight * detJ, 1.);
}
return;
- for (int i=0;i<3 * nbNodes;i++){
- for (int j=0;j<3 * nbNodes;j++){
- printf("%g ",m(i,j));
+ for (int i = 0; i < 3 * nbNodes; i++){
+ for (int j = 0; j < 3 * nbNodes; j++){
+ printf("%g ", m(i, j));
}
printf("\n");
}
printf("\n");
printf("\n");
-
}
void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
@@ -118,9 +117,9 @@ void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
double jac[3][3];
double ff[256];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
-
+
m.scale(0.);
-
+
for (int i = 0; i < npts; i++){
const double u = GP[i].pt[0];
const double v = GP[i].pt[1];
@@ -129,9 +128,9 @@ void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
const double detJ = e->getJacobian(u, v, w, jac);
se->nodalTestFunctions(u, v, w, ff);
for (int j = 0; j < nbNodes; j++){
- m(j) += _volumeForce.x() *ff[j] * weight * detJ *.5;
- m(j + nbNodes) += _volumeForce.y() *ff[j] * weight * detJ *.5;
- m(j + 2 * nbNodes) += _volumeForce.z() *ff[j] * weight * detJ *.5;
+ m(j) += _volumeForce.x() * ff[j] * weight * detJ * .5;
+ m(j + nbNodes) += _volumeForce.y() * ff[j] * weight * detJ * .5;
+ m(j + 2 * nbNodes) += _volumeForce.z() * ff[j] * weight * detJ * .5;
}
- }
+ }
}
diff --git a/Solver/elasticityTerm.h b/Solver/elasticityTerm.h
index ae826d930c..1c0292cf88 100644
--- a/Solver/elasticityTerm.h
+++ b/Solver/elasticityTerm.h
@@ -18,20 +18,20 @@ class elasticityTerm : public femTerm<double, double> {
int _iFieldR, _iFieldC;
SVector3 _volumeForce;
public:
- void setFieldC(int i){_iFieldC = i;}
- void setFieldR(int i){_iFieldR = i;}
+ void setFieldC(int i){_iFieldC = i;}
+ void setFieldR(int i){_iFieldR = i;}
// element matrix size : 3 dofs per vertex
- virtual int sizeOfR(SElement *se) const
+ virtual int sizeOfR(SElement *se) const
{
- return 3 * se->getMeshElement()->getNumVertices();
+ return 3 * se->getMeshElement()->getNumVertices();
}
- virtual int sizeOfC(SElement *se) const
- {
- return 3 * se->getMeshElement()->getNumVertices();
+ virtual int sizeOfC(SElement *se) const
+ {
+ return 3 * se->getMeshElement()->getNumVertices();
}
// order dofs in the local matrix :
// dx1, dx2 ... dxn, dy1, ..., dyn, dz1, ... dzn
- Dof getLocalDofR(SElement *se, int iRow) const
+ Dof getLocalDofR(SElement *se, int iRow) const
{
MElement *e = se->getMeshElement();
int iCompR = iRow / e->getNumVertices();
@@ -39,7 +39,7 @@ class elasticityTerm : public femTerm<double, double> {
return Dof(e->getVertex(ithLocalVertex)->getNum(),
Dof::createTypeWithTwoInts(iCompR, _iFieldR));
}
- Dof getLocalDofC(SElement *se, int iCol) const
+ Dof getLocalDofC(SElement *se, int iCol) const
{
MElement *e = se->getMeshElement();
int iCompC = iCol / e->getNumVertices();
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index 77ba2b55da..5c35bbe4c0 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -24,7 +24,7 @@ class femTerm {
GModel *_gm;
public:
femTerm(GModel *gm) : _gm(gm) {}
- virtual ~femTerm(){}
+ virtual ~femTerm() {}
// return the number of columns of the element matrix
virtual int sizeOfC(SElement *se) const = 0;
// return the number of rows of the element matrix
@@ -34,28 +34,28 @@ class femTerm {
virtual Dof getLocalDofR(SElement *se, int iRow) const = 0;
// default behavior: symmetric
virtual Dof getLocalDofC(SElement *se, int iCol) const
- {
- return getLocalDofR(se, iCol);
+ {
+ return getLocalDofR(se, iCol);
}
// compute the elementary matrix
virtual void elementMatrix(SElement *se, fullMatrix<dataMat> &m) const = 0;
- virtual void elementVector(SElement *se, fullVector<dataVec> &m) const
+ virtual void elementVector(SElement *se, fullVector<dataVec> &m) const
{
m.scale(0.0);
}
// add the contribution from all the elements in the intersection
// of two element groups L and C
- void addToMatrix(dofManager<dataVec, dataMat> &dm,
- groupOfElements &L,
- groupOfElements &C) const
+ void addToMatrix(dofManager<dataVec, dataMat> &dm,
+ groupOfElements &L,
+ groupOfElements &C) const
{
groupOfElements::elementContainer::const_iterator it = L.begin();
for ( ; it != L.end() ; ++it){
MElement *eL = *it;
- if ( &C == &L || C.find(eL) ){
- SElement se(eL);
- addToMatrix(dm, &se);
+ if (&C == &L || C.find(eL)){
+ SElement se(eL);
+ addToMatrix(dm, &se);
}
}
}
@@ -69,8 +69,8 @@ class femTerm {
elementMatrix(se, localMatrix);
addToMatrix(dm, localMatrix, se);
}
- void addToMatrix(dofManager<dataVec, dataMat> &dm,
- fullMatrix<dataMat> &localMatrix,
+ void addToMatrix(dofManager<dataVec, dataMat> &dm,
+ fullMatrix<dataMat> &localMatrix,
SElement *se) const
{
const int nbR = localMatrix.size1();
@@ -116,7 +116,7 @@ class femTerm {
}
void dirichletNodalBC(int physical, int dim, int comp, int field,
const simpleFunction<dataVec> &e,
- dofManager<dataVec,dataMat> &dm)
+ dofManager<dataVec, dataMat> &dm)
{
std::vector<MVertex *> v;
GModel *m = _gm;
@@ -124,14 +124,14 @@ class femTerm {
for (unsigned int i = 0; i < v.size(); i++)
dm.fixVertex(v[i], comp, field, e(v[i]->x(), v[i]->y(), v[i]->z()));
}
- void neumannNodalBC(int physical, int dim, int comp,int field,
+ void neumannNodalBC(int physical, int dim, int comp, int field,
const simpleFunction<dataVec> &fct,
dofManager<dataVec, dataMat> &dm)
{
std::map<int, std::vector<GEntity*> > groups[4];
GModel *m = _gm;
m->getPhysicalGroups(groups);
- std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical);
+ std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical);
if (it == groups[dim].end()) return;
double jac[3][3];
double sf[256];
@@ -143,7 +143,7 @@ class femTerm {
int nbNodes = e->getNumVertices();
int npts;
IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
for (int ip = 0; ip < npts; ip++){
const double u = GP[ip].pt[0];
const double v = GP[ip].pt[1];
@@ -161,10 +161,10 @@ class femTerm {
}
}
- void addToRightHandSide(dofManager<dataVec, dataMat> &dm, groupOfElements &C) const
+ void addToRightHandSide(dofManager<dataVec, dataMat> &dm, groupOfElements &C) const
{
groupOfElements::elementContainer::const_iterator it = C.begin();
- for ( ; it != C.end() ; ++it){
+ for ( ; it != C.end(); ++it){
MElement *eL = *it;
SElement se(eL);
int nbR = sizeOfR(&se);
diff --git a/Solver/groupOfElements.cpp b/Solver/groupOfElements.cpp
index 8fcaea61f5..e409dbc663 100644
--- a/Solver/groupOfElements.cpp
+++ b/Solver/groupOfElements.cpp
@@ -5,11 +5,11 @@
groupOfElements::groupOfElements(GFace*gf)
{
elementFilterTrivial filter;
- addElementary(gf,filter);
-}
+ addElementary(gf, filter);
+}
void groupOfElements::addElementary(GEntity *ge, const elementFilter &filter){
- for (int j=0;j<ge->getNumMeshElements();j++){
+ for (int j = 0; j < ge->getNumMeshElements(); j++){
MElement *e = ge->getMeshElement(j);
if (filter(e)){
insert(e);
@@ -22,8 +22,8 @@ void groupOfElements::addPhysical(int dim, int physical,
std::map<int, std::vector<GEntity*> > groups[4];
GModel::current()->getPhysicalGroups(groups);
std::vector<GEntity*> ent = groups[dim][physical];
- for (int i=0;i<ent.size();i++){
- addElementary(ent[i],filter);
+ for (int i = 0; i < ent.size(); i++){
+ addElementary(ent[i], filter);
}
}
diff --git a/Solver/groupOfElements.h b/Solver/groupOfElements.h
index eae4e9eb58..535861f0aa 100644
--- a/Solver/groupOfElements.h
+++ b/Solver/groupOfElements.h
@@ -7,18 +7,18 @@
class elementFilter {
public:
- virtual bool operator () (MElement *) const = 0;
+ virtual bool operator() (MElement *) const = 0;
};
class elementFilterTrivial : public elementFilter {
public:
- bool operator () (MElement *) const {return true;}
+ bool operator() (MElement *) const {return true;}
};
class groupOfElements {
public:
- typedef std::set<MElement*> elementContainer;
- typedef std::set<MVertex*> vertexContainer;
+ typedef std::set<MElement*> elementContainer;
+ typedef std::set<MVertex*> vertexContainer;
private:
vertexContainer _vertices;
elementContainer _elements;
@@ -26,9 +26,9 @@ private:
public:
groupOfElements (int dim, int physical) {
addPhysical (dim, physical);
- }
+ }
- groupOfElements (GFace*);
+ groupOfElements (GFace*);
void addPhysical(int dim, int physical) {
elementFilterTrivial filter;
@@ -39,38 +39,38 @@ public:
void addPhysical(int dim, int physical, const elementFilter &);
- vertexContainer::const_iterator vbegin () const {
+ vertexContainer::const_iterator vbegin() const {
return _vertices.begin();
}
- vertexContainer::const_iterator vend () const {
+ vertexContainer::const_iterator vend() const {
return _vertices.end();
}
-
- elementContainer::const_iterator begin () const {
+
+ elementContainer::const_iterator begin() const {
return _elements.begin();
}
- elementContainer::const_iterator end () const {
+ elementContainer::const_iterator end() const {
return _elements.end();
}
- size_t size () const {
+ size_t size() const {
return _elements.size();
}
// FIXME : NOT VERY ELEGANT !!!
bool find (MElement *e) const {
- if (e->getParent() && _parents.find(e->getParent()) != _parents.end())return true;
- return (_elements.find(e) != _elements.end()) ;
+ if (e->getParent() && _parents.find(e->getParent()) != _parents.end()) return true;
+ return (_elements.find(e) != _elements.end());
}
inline void insert (MElement *e) {
_elements.insert(e);
if (e->getParent()){
- _parents.insert(e->getParent());
- for (int i=0;i<e->getParent()->getNumVertices();i++){
- _vertices.insert(e->getParent()->getVertex(i));
+ _parents.insert(e->getParent());
+ for (int i = 0; i < e->getParent()->getNumVertices(); i++){
+ _vertices.insert(e->getParent()->getVertex(i));
}
}
else{
- for (int i=0;i<e->getNumVertices();i++){
- _vertices.insert(e->getVertex(i));
+ for (int i = 0; i < e->getNumVertices(); i++){
+ _vertices.insert(e->getVertex(i));
}
}
}
diff --git a/Solver/linearSystemCSR.cpp b/Solver/linearSystemCSR.cpp
index 65d98a9e4c..3f83b68c98 100644
--- a/Solver/linearSystemCSR.cpp
+++ b/Solver/linearSystemCSR.cpp
@@ -11,8 +11,8 @@
#include "GmshMessage.h"
#include "linearSystemCSR.h"
-#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
-#define SWAPI(a,b) tempi=(a);(a)=(b);(b)=tempi;
+#define SWAP(a, b) temp = (a); (a) = (b); (b) = temp;
+#define SWAPI(a, b) tempi = (a); (a) = (b); (b) = tempi;
static void *CSRMalloc(size_t size)
{
@@ -25,7 +25,7 @@ static void *CSRMalloc(size_t size)
static void *CSRRealloc(void *ptr, size_t size)
{
if (!size) return(NULL);
- ptr = realloc(ptr,size);
+ ptr = realloc(ptr, size);
return(ptr);
}
@@ -49,20 +49,20 @@ static void CSRList_Realloc(CSRList_T *liste, int n)
static CSRList_T *CSRList_Create(int n, int incr, int size)
{
CSRList_T *liste;
-
- if (n <= 0) n = 1 ;
+
+ if (n <= 0) n = 1 ;
if (incr <= 0) incr = 1;
-
+
liste = (CSRList_T *)CSRMalloc(sizeof(CSRList_T));
-
+
liste->nmax = 0;
liste->incr = incr;
liste->size = size;
liste->n = 0;
liste->isorder = 0;
liste->array = NULL;
-
- CSRList_Realloc(liste,n);
+
+ CSRList_Realloc(liste, n);
return(liste);
}
@@ -77,7 +77,7 @@ static void CSRList_Delete(CSRList_T *liste)
void CSRList_Add(CSRList_T *liste, void *data)
{
liste->n++;
- CSRList_Realloc(liste,liste->n);
+ CSRList_Realloc(liste, liste->n);
liste->isorder = 0;
memcpy(&liste->array[(liste->n - 1) * liste->size], data, liste->size);
}
@@ -99,165 +99,165 @@ void linearSystemCSR<double>::allocate(int nbRows)
delete _b;
delete[] something;
}
-
+
if(nbRows == 0){
- _a = 0;
- _ai = 0;
- _ptr = 0;
- _jptr = 0;
+ _a = 0;
+ _ai = 0;
+ _ptr = 0;
+ _jptr = 0;
_b = 0;
_x = 0;
sorted = false;
something = 0;
return;
}
-
+
_a = CSRList_Create(nbRows, nbRows, sizeof(double));
_ai = CSRList_Create(nbRows, nbRows, sizeof(INDEX_TYPE));
_ptr = CSRList_Create(nbRows, nbRows, sizeof(INDEX_TYPE));
_jptr = CSRList_Create(nbRows + 1, nbRows, sizeof(INDEX_TYPE));
-
+
something = new char[nbRows];
-
+
for (int i = 0; i < nbRows; i++) something[i] = 0;
-
+
_b = new std::vector<double>(nbRows);
_x = new std::vector<double>(nbRows);
}
const int NSTACK = 50;
-const unsigned int M_sort2 = 7;
+const unsigned int M_sort2 = 7;
static void free_ivector(int *v, long nl, long nh)
{
// free an int vector allocated with ivector()
- free((char*)(v+nl-1));
+ free((char*)(v + nl - 1));
}
static int *ivector(long nl, long nh)
{
- // allocate an int vector with subscript range v[nl..nh]
- int *v;
- v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
+ // allocate an int vector with subscript range v[nl..nh]
+ int *v;
+ v = (int *)malloc((size_t) ((nh - nl + 2) * sizeof(int)));
if (!v) fprintf(stderr, "allocation failure in ivector()\n");
- return v-nl+1;
+ return v - nl + 1;
}
-static int cmpij(INDEX_TYPE ai,INDEX_TYPE aj,INDEX_TYPE bi,INDEX_TYPE bj)
+static int cmpij(INDEX_TYPE ai, INDEX_TYPE aj, INDEX_TYPE bi, INDEX_TYPE bj)
{
- if(ai<bi)return -1;
- if(ai>bi)return 1;
- if(aj<bj)return -1;
- if(aj>bj)return 1;
+ if(ai < bi) return -1;
+ if(ai > bi) return 1;
+ if(aj < bj) return -1;
+ if(aj > bj) return 1;
return 0;
}
template <class scalar>
static void _sort2_xkws(unsigned long n, double arr[], INDEX_TYPE ai[], INDEX_TYPE aj[])
{
- unsigned long i,ir=n,j,k,l=1;
- int *istack,jstack=0;
+ unsigned long i, ir = n, j, k, l = 1;
+ int *istack, jstack = 0;
INDEX_TYPE tempi;
- scalar a,temp;
- int b,c;
-
- istack=ivector(1,NSTACK);
+ scalar a, temp;
+ int b, c;
+
+ istack = ivector(1, NSTACK);
for (;;) {
- if (ir-l < M_sort2) {
- for (j=l+1;j<=ir;j++) {
- a=arr[j -1];
- b=ai[j -1];
- c=aj[j -1];
- for (i=j-1;i>=1;i--) {
- if (cmpij(ai[i -1],aj[i -1],b,c) <= 0) break;
- arr[i+1 -1]=arr[i -1];
- ai[i+1 -1]=ai[i -1];
- aj[i+1 -1]=aj[i -1];
+ if (ir - l < M_sort2) {
+ for (j = l + 1; j <= ir; j++) {
+ a = arr[j -1];
+ b = ai[j -1];
+ c = aj[j -1];
+ for (i = j - 1; i >= 1; i--) {
+ if (cmpij(ai[i - 1], aj[i - 1], b, c) <= 0) break;
+ arr[i + 1 - 1] = arr[i - 1];
+ ai[i + 1 - 1] = ai[i - 1];
+ aj[i + 1 - 1] = aj[i - 1];
}
- arr[i+1 -1]=a;
- ai[i+1 -1]=b;
- aj[i+1 -1]=c;
+ arr[i + 1 - 1] = a;
+ ai[i + 1 - 1] = b;
+ aj[i + 1 - 1] = c;
}
if (!jstack) {
- free_ivector(istack,1,NSTACK);
+ free_ivector(istack, 1, NSTACK);
return;
}
- ir=istack[jstack];
- l=istack[jstack-1];
+ ir = istack[jstack];
+ l = istack[jstack - 1];
jstack -= 2;
- }
+ }
else {
- k=(l+ir) >> 1;
- SWAP(arr[k -1],arr[l+1 -1])
- SWAPI(ai[k -1],ai[l+1 -1])
- SWAPI(aj[k -1],aj[l+1 -1])
- if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[ir -1],aj[ir -1])>0){
- SWAP(arr[l+1 -1],arr[ir -1])
- SWAPI(ai[l+1 -1],ai[ir -1])
- SWAPI(aj[l+1 -1],aj[ir -1])
- }
- if (cmpij(ai[l -1],aj[l -1],ai[ir -1],aj[ir -1])>0){
- SWAP(arr[l -1],arr[ir -1])
- SWAPI(ai[l -1],ai[ir -1])
- SWAPI(aj[l -1],aj[ir -1])
- }
- if (cmpij(ai[l+1 -1],aj[l+1 -1],ai[l -1],aj[l -1])>0){
- SWAP(arr[l+1 -1],arr[l -1])
- SWAPI(ai[l+1 -1],ai[l -1])
- SWAPI(aj[l+1 -1],aj[l -1])
- }
- i=l+1;
- j=ir;
- a=arr[l -1];
- b=ai[l -1];
- c=aj[l -1];
+ k = (l + ir) >> 1;
+ SWAP(arr[k - 1], arr[l + 1 - 1])
+ SWAPI(ai[k - 1], ai[l + 1 - 1])
+ SWAPI(aj[k - 1], aj[l + 1 - 1])
+ if (cmpij(ai[l + 1 - 1], aj[l + 1 - 1], ai[ir - 1], aj[ir - 1]) > 0){
+ SWAP(arr[l + 1 - 1], arr[ir - 1])
+ SWAPI(ai[l + 1 - 1], ai[ir - 1])
+ SWAPI(aj[l + 1 - 1], aj[ir - 1])
+ }
+ if (cmpij(ai[l - 1], aj[l - 1], ai[ir - 1], aj[ir - 1]) > 0){
+ SWAP(arr[l - 1], arr[ir - 1])
+ SWAPI(ai[l - 1], ai[ir - 1])
+ SWAPI(aj[l - 1], aj[ir - 1])
+ }
+ if (cmpij(ai[l + 1 - 1], aj[l + 1 - 1], ai[l - 1], aj[l - 1]) > 0){
+ SWAP(arr[l + 1 - 1], arr[l - 1])
+ SWAPI(ai[l + 1 - 1], ai[l - 1])
+ SWAPI(aj[l + 1 - 1], aj[l - 1])
+ }
+ i = l + 1;
+ j = ir;
+ a = arr[l - 1];
+ b = ai[l - 1];
+ c = aj[l - 1];
for (;;) {
- do i++; while (cmpij(ai[i -1],aj[i -1],b,c) < 0);
- do j--; while (cmpij(ai[j -1],aj[j -1],b,c) > 0);
+ do i++; while (cmpij(ai[i -1], aj[i -1], b, c) < 0);
+ do j--; while (cmpij(ai[j -1], aj[j -1], b, c) > 0);
if (j < i) break;
- SWAP(arr[i -1],arr[j -1])
- SWAPI(ai[i -1],ai[j -1])
- SWAPI(aj[i -1],aj[j -1])
- }
- arr[l -1]=arr[j -1];
- arr[j -1]=a;
- ai[l -1]=ai[j -1];
- ai[j -1]=b;
- aj[l -1]=aj[j -1];
- aj[j -1]=c;
+ SWAP(arr[i - 1], arr[j - 1])
+ SWAPI(ai[i - 1], ai[j - 1])
+ SWAPI(aj[i - 1], aj[j - 1])
+ }
+ arr[l - 1] = arr[j - 1];
+ arr[j - 1] = a;
+ ai[l - 1] = ai[j - 1];
+ ai[j - 1] = b;
+ aj[l - 1] = aj[j - 1];
+ aj[j - 1] = c;
jstack += 2;
if (jstack > NSTACK) {
Msg::Fatal("NSTACK too small while sorting the columns of the matrix");
throw;
}
- if (ir-i+1 >= j-l) {
- istack[jstack]=ir;
- istack[jstack-1]=i;
- ir=j-1;
- }
+ if (ir - i + 1 >= j - l) {
+ istack[jstack] = ir;
+ istack[jstack - 1] = i;
+ ir = j - 1;
+ }
else {
- istack[jstack]=j-1;
- istack[jstack-1]=l;
- l=i;
+ istack[jstack] = j - 1;
+ istack[jstack - 1] = l;
+ l = i;
}
}
}
}
template <class scalar>
-void sortColumns_(int NbLines,
- int nnz,
- INDEX_TYPE *ptr,
- INDEX_TYPE *jptr,
- INDEX_TYPE *ai,
- scalar *a)
+void sortColumns_(int NbLines,
+ int nnz,
+ INDEX_TYPE *ptr,
+ INDEX_TYPE *jptr,
+ INDEX_TYPE *ai,
+ scalar *a)
{
// replace pointers by lines
int *count = new int[NbLines];
-
- for(int i=0;i<NbLines;i++){
+
+ for(int i = 0; i < NbLines; i++){
count[i] = 0;
- INDEX_TYPE _position = jptr[i];
+ INDEX_TYPE _position = jptr[i];
while(1){
count[i]++;
INDEX_TYPE _position_temp = _position;
@@ -265,20 +265,20 @@ void sortColumns_(int NbLines,
ptr[_position_temp] = i;
if (_position == 0) break;
}
- }
- _sort2_xkws<double>(nnz,a,ptr,ai);
+ }
+ _sort2_xkws<double>(nnz, a, ptr, ai);
jptr[0] = 0;
- for(int i=1;i<=NbLines;i++){
- jptr[i] = jptr[i-1] + count[i-1];
+ for(int i = 1; i <= NbLines; i++){
+ jptr[i] = jptr[i - 1] + count[i - 1];
}
-
- for(int i=0;i<NbLines;i++){
- for (int j= jptr[i] ; j<jptr[i+1]-1 ; j++){
- ptr[j] = j+1;
+
+ for(int i = 0; i < NbLines; i++){
+ for (int j = jptr[i]; j < jptr[i + 1] - 1; j++){
+ ptr[j] = j + 1;
}
- ptr[jptr[i+1]] = 0;
+ ptr[jptr[i + 1]] = 0;
}
-
+
delete[] count;
}
@@ -293,17 +293,17 @@ int linearSystemCSRGmm<double>::systemSolve()
sortColumns_(_b->size(),
CSRList_Nbr(_a),
(INDEX_TYPE *) _ptr->array,
- (INDEX_TYPE *) _jptr->array,
- (INDEX_TYPE *) _ai->array,
+ (INDEX_TYPE *) _jptr->array,
+ (INDEX_TYPE *) _ai->array,
(double*) _a->array);
sorted = true;
-
- gmm::csr_matrix_ref<double*,INDEX_TYPE *,INDEX_TYPE *, 0>
+
+ gmm::csr_matrix_ref<double*, INDEX_TYPE *, INDEX_TYPE *, 0>
ref((double*)_a->array, (INDEX_TYPE *) _ai->array,
(INDEX_TYPE *)_jptr->array, _b->size(), _b->size());
- gmm::csr_matrix<double,0> M;
+ gmm::csr_matrix<double, 0> M;
M.init_with(ref);
-
+
gmm::ildltt_precond<gmm::csr_matrix<double, 0> > P(M, 10, 1.e-10);
gmm::iteration iter(_prec);
iter.set_noisy(_noisy);
@@ -327,12 +327,12 @@ int linearSystemCSRTaucs<double>::systemSolve()
sortColumns_(_b->size(),
CSRList_Nbr(_a),
(INDEX_TYPE *) _ptr->array,
- (INDEX_TYPE *) _jptr->array,
- (INDEX_TYPE *) _ai->array,
+ (INDEX_TYPE *) _jptr->array,
+ (INDEX_TYPE *) _ai->array,
(double*) _a->array);
}
sorted = true;
-
+
taucs_ccs_matrix myVeryCuteTaucsMatrix;
myVeryCuteTaucsMatrix.n = myVeryCuteTaucsMatrix.m = _b->size();
//myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)_ptr->array;
@@ -340,23 +340,23 @@ int linearSystemCSRTaucs<double>::systemSolve()
myVeryCuteTaucsMatrix.rowind = (INDEX_TYPE*)_ai->array;
myVeryCuteTaucsMatrix.colptr = (INDEX_TYPE*)_jptr->array;
myVeryCuteTaucsMatrix.values.d = (double*)_a->array;
- myVeryCuteTaucsMatrix.flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
-
- char* options[] = { "taucs.factor.LLT=true", NULL };
+ myVeryCuteTaucsMatrix.flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
+
+ char* options[] = { "taucs.factor.LLT=true", NULL };
clock_t t1 = clock();
- int result = taucs_linsolve(&myVeryCuteTaucsMatrix,
- NULL,
- 1,
+ int result = taucs_linsolve(&myVeryCuteTaucsMatrix,
+ NULL,
+ 1,
&(*_x)[0],
&(*_b)[0],
options,
- NULL);
+ NULL);
clock_t t2 = clock();
- Msg::Info("TAUCS has solved %d unknowns in %8.3f seconds",
+ Msg::Info("TAUCS has solved %d unknowns in %8.3f seconds",
_b->size(), (double)(t2 - t1) / CLOCKS_PER_SEC);
if (result != TAUCS_SUCCESS){
Msg::Error("Taucs Was Not Successfull %d", result);
- }
+ }
return 1;
}
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 517334a52d..338dca15e4 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -63,12 +63,14 @@ inline double distance(const DI_CuttingPoint &p1, const DI_CuttingPoint &p2) {
inline DI_Point middle (const DI_Point &p1, const DI_Point &p2) {
return DI_Point ((p1.x() + p2.x()) / 2, (p1.y() + p2.y()) / 2, (p1.z() + p2.z()) / 2);
}
-inline DI_Point middle (const DI_Point &p1, const DI_Point &p2, const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
+inline DI_Point middle (const DI_Point &p1, const DI_Point &p2,
+ const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
return DI_Point ((p1.x() + p2.x()) / 2, (p1.y() + p2.y()) / 2, (p1.z() + p2.z()) / 2, e, RPNi);
}
// barycentre
-inline DI_Point barycenter (const DI_Element *el, const DI_Element *e, const std::vector<const gLevelset *> RPN) {
+inline DI_Point barycenter (const DI_Element *el,
+ const DI_Element *e, const std::vector<const gLevelset *> RPN) {
double x[3] = {0., 0., 0.};
int n;
for (n = 0; n < el->nbVert(); n++) {
@@ -233,7 +235,8 @@ double qualityTri(const DI_Point &p0, const DI_Point &p1, const DI_Point &p2) {
}
// Return the best quality for a quad cut into 2 triangles
-int bestQuality(const DI_Point &p1, const DI_Point &p2, const DI_Point &p3, const DI_Point &p4, DI_Triangle &t1, DI_Triangle &t2) {
+int bestQuality(const DI_Point &p1, const DI_Point &p2, const DI_Point &p3, const DI_Point &p4,
+ DI_Triangle &t1, DI_Triangle &t2) {
int cut = 0;
double qual11 = qualityTri(p1, p2, p3);
double qual12 = qualityTri(p1, p3, p4);
@@ -448,7 +451,8 @@ int bestQuality (const DI_Point &p0, const DI_Point &p1, const DI_Point &p2,
}
// computes the intersection between a level set and a linear edge
-DI_Point Newton(const DI_Point &p1, const DI_Point &p2, const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
+DI_Point Newton(const DI_Point &p1, const DI_Point &p2,
+ const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
double eps = -p1.ls() / (p2.ls() - p1.ls());
DI_Point p(p1.x() + eps * (p2.x() - p1.x()), p1.y() + eps * (p2.y() - p1.y()), p1.z() + eps * (p2.z() - p1.z()));
double pls = e->evalLs(p.x(), p.y(), p.z());
@@ -468,8 +472,10 @@ DI_Point Newton(const DI_Point &p1, const DI_Point &p2, const DI_Element *e, con
return p;
}
-// compute the position of the middle of a quadratic edge (intersection between the bisector of the linear edge and the levelset)
-DI_Point quadMidNode(const DI_Point &p1, const DI_Point &p2, const DI_Point &pf, const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
+// compute the position of the middle of a quadratic edge
+//(intersection between the bisector of the linear edge and the levelset)
+DI_Point quadMidNode(const DI_Point &p1, const DI_Point &p2, const DI_Point &pf,
+ const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
// middle of the edge between the 2 cutting points
DI_Point midEN((p1.x() + p2.x()) / 2., (p1.y() + p2.y()) / 2., (p1.z() + p2.z()) / 2.);
midEN.addLs(e);
@@ -493,10 +499,10 @@ DI_Point quadMidNode(const DI_Point &p1, const DI_Point &p2, const DI_Point &pf,
return Newton(midEN, pt, e, RPNi);
}
-// --------------------------------------------------------------------------------------------------------------------------------------------------
-// --------------------------------------------------------------------------------------------------------------------------------------------------
+// ------------------------------------------------------------------------------------------------
+// ------------------------------------------------------------------------------------------------
-// DI_Point methods ----------------------------------------------------------------------------------------------
+// DI_Point methods -------------------------------------------------------------------------------
DI_Point::DI_Point (double x, double y, double z, gLevelset &ls) : x_(x), y_(y), z_(z) {
addLs(ls(x, y, z));
}
@@ -537,7 +543,7 @@ bool DI_Point::equal(const DI_Point &p) const {
return (fabs(x() - p.x()) < EQUALITY_TOL && fabs(y() - p.y()) < EQUALITY_TOL && fabs(z() - p.z()) < EQUALITY_TOL);
}
-// DI_IntegrationPoint methods -----------------------------------------------------------------------------------------------------
+// DI_IntegrationPoint methods --------------------------------------------------------------------
void DI_IntegrationPoint::computeLs (const DI_Element *e, const std::vector<const gLevelset*> RPN) {
int prim = 0; std::vector<double> Ls;
for(int l = 0; l < (int)RPN.size(); l++) {
@@ -553,7 +559,7 @@ void DI_IntegrationPoint::computeLs (const DI_Element *e, const std::vector<cons
setLs(Ls.back());
}
-// DI_CuttingPoint methods -----------------------------------------------------------------------------------------------------
+// DI_CuttingPoint methods ------------------------------------------------------------------------
void DI_CuttingPoint::chooseLs (const gLevelset *Lsi) {
if(Ls.size() < 2) return;
double ls1 = Ls[Ls.size() - 2], ls2 = Ls[Ls.size() - 1];
@@ -565,8 +571,9 @@ bool DI_CuttingPoint::equal (const DI_CuttingPoint &p) const {
return (fabs(x() - p.x()) < EQUALITY_TOL && fabs(y() - p.y()) < EQUALITY_TOL && fabs(z() - p.z()) < EQUALITY_TOL);
}
-// DI_Element methods ------------------------------------------------------------------------------------------------
-DI_Element::DI_Element(const DI_Element &cp) : lsTag_(cp.lsTag()), polOrder_(cp.getPolynomialOrder()), integral_(cp.integral()) {
+// DI_Element methods -----------------------------------------------------------------------------
+DI_Element::DI_Element(const DI_Element &cp) : lsTag_(cp.lsTag()), polOrder_(cp.getPolynomialOrder()),
+ integral_(cp.integral()) {
//printf("constructeur de copie d'element %d : ",cp.type()); cout << &cp << "->" << this << endl;
pts_ = new DI_Point* [cp.nbVert()];
for(int i = 0; i < cp.nbVert(); i++)
@@ -705,7 +712,8 @@ void DI_Element::clearLs() {
for(int i = 0; i < nbMid(); i++)
(mid_[i]->Ls).clear();
}
-bool DI_Element::addQuadEdge (int edge, DI_Point *xm, const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
+bool DI_Element::addQuadEdge (int edge, DI_Point *xm,
+ const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
/*if(edge >= nbEdg()) {printf("wrong number (%d) for quadratic edge for a ", edge); print(); return false;}
int s1, s2; vert(edge, s1, s2);
bool quad0 = isQuad();
@@ -800,8 +808,9 @@ void DI_Element::mappingEl (DI_Element *el) const {
}
el->computeIntegral();
}
-void DI_Element::integrationPoints (const int polynomialOrder, const DI_Element *loc, const DI_Element *e,
- const std::vector<const gLevelset *> RPN, std::vector<DI_IntegrationPoint> &ip) const {
+void DI_Element::integrationPoints (const int polynomialOrder, const DI_Element *loc,
+ const DI_Element *e, const std::vector<const gLevelset *> RPN,
+ std::vector<DI_IntegrationPoint> &ip) const {
std::vector<DI_IntegrationPoint> ip_ref;
getRefIntegrationPoints(polynomialOrder, ip_ref);
for (int j = 0; j < (int)ip_ref.size(); j++) {
@@ -819,7 +828,8 @@ void DI_Element::integrationPoints (const int polynomialOrder, const DI_Element
ip.push_back(ip_ref[j]);
}
}
-void DI_Element::getCuttingPoints (const DI_Element *e, const std::vector<const gLevelset *> RPNi, std::vector<DI_CuttingPoint> &cp) const {
+void DI_Element::getCuttingPoints (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
+ std::vector<DI_CuttingPoint> &cp) const {
int s1, s2;
for(int i = 0; i < nbEdg(); i++){
vert(i, s1, s2);
@@ -1009,7 +1019,7 @@ void DI_Element::printls() const {
printf("tag=%d\n", lsTag());
}
-// DI_Line methods -----------------------------------------------------------------------------------------------------------
+// DI_Line methods --------------------------------------------------------------------------------
void DI_Line::computeIntegral() {
switch (getPolynomialOrder()) {
case 1 :
@@ -1055,19 +1065,8 @@ void DI_Line::getGradShapeFunctions (const double u, const double v, const doubl
printf("Order %d line function space not implemented ", order); print();
}
}
-/*double ln_detJ1 (const double &x, const double &x0, const double &x1) {
- return (-x0 + x1) / 2.; // detJ is constant in a linear line
-}
-double ln_detJ2 (const double &x, const double &x0, const double &x1, const double &x2) {
- return x0 * (2. * x - 1.) / 2. + x1 * (2. * x + 1.) / 2. + x2 * (-2. * x);
-}
-double DI_Line::detJ (const double &xP, const double &yP, const double &zP) const {
- assert(yP == 0 && zP == 0 && y(0) == 0 && y(1) == 0 && y2(0) == 0 && z(0) == 0 && z(1) == 0 && z2(0) == 0);
- if(!isQuad()) return ln_detJ1(xP, x(0), x(1));
- return ln_detJ2(xP, x(0), x(1), x2(0));
-}*/
-// DI_Triangle methods -----------------------------------------------------------------------------------------------------------
+// DI_Triangle methods ----------------------------------------------------------------------------
void DI_Triangle::computeIntegral() {
switch (getPolynomialOrder()) {
case 1 :
@@ -1102,7 +1101,8 @@ void DI_Triangle::getShapeFunctions (double u, double v, double w, double s[], i
printf("Order %d triangle function space not implemented ", order); print();
}
}
-void DI_Triangle::getGradShapeFunctions (const double u, const double v, const double w, double s[][3], int ord) const {
+void DI_Triangle::getGradShapeFunctions (const double u, const double v, const double w,
+ double s[][3], int ord) const {
int order = (ord == -1) ? getPolynomialOrder() : ord;
switch (order) {
case 1 :
@@ -1137,24 +1137,8 @@ void DI_Triangle::getGradShapeFunctions (const double u, const double v, const d
double DI_Triangle::quality() const {
return qualityTri(pt(0), pt(1), pt(2));
}
-/*double tr_detJ1 (const double &x, const double &y,
- const double &x0, const double &y0, const double &x1, const double &y1, const double &x2, const double &y2) {
- return (x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0); // detJ is constant in a linear triangle
-}
-double tr_detJ2 (const double &x, const double &y,
- const double &x0, const double &y0, const double &x1, const double &y1, const double &x2, const double &y2,
- const double &x3, const double &y3, const double &x4, const double &y4, const double &x5, const double &y5) {
- return (x0*(-3.+4.*x+4.*y)+x1*(4.*x-1.)+x3*4.*(1.-2.*x-y)+x4*4.*y-x5*4.*y) * (y0*(-3.+4.*x+4.*y)+y2*(4.*y-1.)-y3*4.*x+y4*4.*x+y5*4.*(1.-x-2.*y))
- -(x0*(-3.+4.*x+4.*y)+x2*(4.*y-1.)-x3*4.*x+x4*4.*x+x5*4.*(1.-x-2.*y)) * (y0*(-3.+4.*x+4.*y)+y1*(4.*x-1.)+y3*4.*(1.-2.*x-y)+y4*4.*y-y5*4.*y);
-}
-double DI_Triangle::detJ (const double &xP, const double &yP, const double &zP) const {
- if(!(z(0) == 0 && z(1) == 0 && z(2) == 0 && z2(0) == 0 && z2(1) == 0 && z2(2) == 0)) print();
- assert(zP == 0 && z(0) == 0 && z(1) == 0 && z(2) == 0 && z2(0) == 0 && z2(1) == 0 && z2(2) == 0);
- if(!isQuad()) return tr_detJ1(xP, yP, x(0), y(0), x(1), y(1), x(2), y(2));
- return tr_detJ2(xP, yP, x(0), y(0), x(1), y(1), x(2), y(2), x2(0), y2(0), x2(1), y2(1), x2(2), y2(2));
-}*/
-// DI_Quad methods -------------------------------------------------------------------------------------------------------------------
+// DI_Quad methods --------------------------------------------------------------------------------
void DI_Quad::computeIntegral() {
switch (getPolynomialOrder()) {
case 1 :
@@ -1236,32 +1220,8 @@ void DI_Quad::getGradShapeFunctions (const double u, const double v, const doubl
printf("Order %d quadrangle function space not implemented ", order); print();
}
}
-/*double q_detJ1 (const double &x, const double &y,
- const double &x0, const double &y0, const double &x1, const double &y1,
- const double &x2, const double &y2, const double &x3, const double &y3) {
- return (x0*(-1.+y)+x1*(1.-y)+x2*(1.+y)+x3*(-1.-y))*(y0*(x-1.)+y1*(-1.-x)+y2*(1.+x)+y3*(1.-x))/16.
- -(x0*(x-1.)+x1*(-1.-x)+x2*(1.+x)+x3*(1.-x))*(y0*(-1.+y)+y1*(1.-y)+y2*(1.+y)+y3*(-1.-y))/16.;
-}
-double q_detJ2 (const double &x, const double &y,
- const double &x0, const double &y0, const double &x1, const double &y1, const double &x2, const double &y2,
- const double &x3, const double &y3, const double &x4, const double &y4, const double &x5, const double &y5,
- const double &x6, const double &y6, const double &x7, const double &y7, const double &x8, const double &y8) {
- return (x0*y*(1-y)*(1-2*x)-x1*y*(1-y)*(1+2*x)+x2*y*(1+y)*(1+2*x)-x3*y*(1+y)*(1-2*x)+
- x4*4*y*(1-y)*x+x5*2*(1-y)*(1+y)*(1+2*x)-x6*4*y*(1+y)*x-x7*2*(1-y)*(1+y)*(1-2*x)-x8*8*(1-y)*(1+y)*x)
- *(y0*x*(1-x)*(1-2*y)-y1*x*(1+x)*(1-2*y)+y2*x*(1+x)*(1+2*y)-y3*x*(1-x)*(1+2*y)-
- y4*2*(1-x)*(1+x)*(1-2*y)-y5*4*x*(1+x)*y+y6*2*(1-x)*(1+x)*(1+2*y)+y7*4*x*(1-x)*y-y8*8*(1-x)*(1+x)*y)/16.
- -(x0*x*(1-x)*(1-2*y)-x1*x*(1+x)*(1-2*y)+x2*x*(1+x)*(1+2*y)-x3*x*(1-x)*(1+2*y)-
- x4*2*(1-x)*(1+x)*(1-2*y)-x5*4*x*(1+x)*y+x6*2*(1-x)*(1+x)*(1+2*y)+x7*4*x*(1-x)*y-x8*8*(1-x)*(1+x)*y)
- *(y0*y*(1-y)*(1-2*x)-y1*y*(1-y)*(1+2*x)+y2*y*(1+y)*(1+2*x)-y3*y*(1+y)*(1-2*x)+
- y4*4*y*(1-y)*x+y5*2*(1-y)*(1+y)*(1+2*x)-y6*4*y*(1+y)*x-y7*2*(1-y)*(1+y)*(1-2*x)-y8*8*(1-y)*(1+y)*x)/16.;
-}
-double DI_Quad::detJ (const double &xP, const double &yP, const double &zP) const {
- assert (zP == 0 && z(0) == 0 && z(1) == 0 && z(2) == 0 && z(3) == 0 && z2(0) == 0 && z2(1) == 0 && z2(2) == 0 && z2(3) == 0 && z2(4) == 0);
- if(!isQuad()) return q_detJ1(xP, yP, x(0), y(0), x(1), y(1), x(2), y(2), x(3), y(3));
- return q_detJ2(xP, yP, x(0), y(0), x(1), y(1), x(2), y(2), x(3), y(3), x2(0), y2(0), x2(1), y2(1), x2(2), y2(2), x2(3), y2(3), x2(4), y2(4));
-}*/
-// DI_Tetra methods ----------------------------------------------------------------------------------------------------------------------
+// DI_Tetra methods -------------------------------------------------------------------------------
void DI_Tetra::computeIntegral() {
integral_ = TetraVol(pt(0), pt(1), pt(2), pt(3));
}
@@ -1339,44 +1299,7 @@ double DI_Tetra::quality() const {
return qualityTet(x(0), y(0), z(0), x(1), y(1), z(1), x(2), y(2), z(2), x(3), y(3), z(3));
}
-/*double tet_detJ1 (const double &x, const double &y, const double &z,
- const double &x0, const double &y0, const double &z0, const double &x1, const double &y1, const double &z1,
- const double &x2, const double &y2, const double &z2, const double &x3, const double &y3, const double &z3) {
- return det3(x1-x0, x2-x0, x3-x0, y1-y0, y2-y0, y3-y0, z1-z0, z2-z0, z3-z0); // detJ is constant in a linear tetrahedron
-}
-double tet_detJ2 (const double &x, const double &y, const double &z,
- const double &x0, const double &y0, const double &z0, const double &x1, const double &y1, const double &z1,
- const double &x2, const double &y2, const double &z2, const double &x3, const double &y3, const double &z3,
- const double &x4, const double &y4, const double &z4, const double &x5, const double &y5, const double &z5,
- const double &x6, const double &y6, const double &z6, const double &x7, const double &y7, const double &z7,
- const double &x8, const double &y8, const double &z8, const double &x9, const double &y9, const double &z9) {
- return det3(x0*(-3+4*x+4*y+4*z)+x1*(4*x-1)+x4*4*(1-2*x-y-z)-x5*4*y-x6*4*z+x7*4*y+x9*4*z,
- x0*(-3+4*x+4*y+4*z)+x2*(4*y-1)-x4*4*x+x5*4*(1-x-2*y-z)-x6*4*z+x7*4*x+x8*4*z,
- x0*(-3+4*x+4*y+4*z)+x3*(4*z-1)-x4*4*x-x5*4*y+x6*4*(1-x-y-2*z)+x8*4*y+x9*4*x,
- y0*(-3+4*x+4*y+4*z)+y1*(4*x-1)+y4*4*(1-2*x-y-z)-y5*4*y-y6*4*z+y7*4*y+y9*4*z,
- y0*(-3+4*x+4*y+4*z)+y2*(4*y-1)-y4*4*x+y5*4*(1-x-2*y-z)-y6*4*z+y7*4*x+y8*4*z,
- y0*(-3+4*x+4*y+4*z)+y3*(4*z-1)-y4*4*x-y5*4*y+y6*4*(1-x-y-2*z)+y8*4*y+y9*4*x,
- z0*(-3+4*x+4*y+4*z)+z1*(4*x-1)+z4*4*(1-2*x-y-z)-z5*4*y-z6*4*z+z7*4*y+z9*4*z,
- z0*(-3+4*x+4*y+4*z)+z2*(4*y-1)-z4*4*x+z5*4*(1-x-2*y-z)-z6*4*z+z7*4*x+z8*4*z,
- z0*(-3+4*x+4*y+4*z)+z3*(4*z-1)-z4*4*x-z5*4*y+z6*4*(1-x-y-2*z)+z8*4*y+z9*4*x);
-}
-double DI_Tetra::detJ (const double &xP, const double &yP, const double &zP) const {
- if(!isQuad()) return tet_detJ1(xP, yP, zP, x(0), y(0), z(0), x(1), y(1), z(1), x(2), y(2), z(2), x(3), y(3), z(3));
- return tet_detJ2(xP, yP, zP, x(0), y(0), z(0), x(1), y(1), z(1), x(2), y(2), z(2), x(3), y(3), z(3),
- x2(0), y2(0), z2(0), x2(1), y2(1), z2(1), x2(2), y2(2), z2(2), x2(3), y2(3), z2(3), x2(4), y2(4), z2(4), x2(5), y2(5), z2(5));
-}*/
-/*void DI_Tetra::quad(const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
- mid_ = new DI_Point*[6];
- mid_[0] = new DI_Point((x(0)+x(1))/2., (y(0)+y(1))/2., (z(0)+z(1))/2., e, RPNi);
- mid_[1] = new DI_Point((x(0)+x(2))/2., (y(0)+y(2))/2., (z(0)+z(2))/2., e, RPNi);
- mid_[2] = new DI_Point((x(0)+x(3))/2., (y(0)+y(3))/2., (z(0)+z(3))/2., e, RPNi);
- mid_[3] = new DI_Point((x(1)+x(2))/2., (y(1)+y(2))/2., (z(1)+z(2))/2., e, RPNi);
- mid_[4] = new DI_Point((x(2)+x(3))/2., (y(2)+y(3))/2., (z(2)+z(3))/2., e, RPNi);
- mid_[5] = new DI_Point((x(3)+x(1))/2., (y(3)+y(1))/2., (z(3)+z(1))/2., e, RPNi);
- quad_ = true;
-}*/
-
-// Hexahedron methods ------------------------------------------------------------------------------------------------------------
+// Hexahedron methods -----------------------------------------------------------------------------
void DI_Hexa::computeIntegral() {
integral_ = TetraVol(pt(0), pt(1), pt(3), pt(4)) + TetraVol(pt(1), pt(4), pt(5), pt(7))
+ TetraVol(pt(1), pt(3), pt(4), pt(7)) + TetraVol(pt(2), pt(5), pt(6), pt(7))
@@ -1461,30 +1384,6 @@ void DI_Hexa::getGradShapeFunctions (const double u, const double v, const doubl
printf("Order %d hexahedron function space not implemented ", order); print();
}
}
-/*void DI_Hexa::quad(const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
- mid_ = new DI_Point*[19];
- mid_[0] = new DI_Point((x(0)+x(1))/2., (y(0)+y(1))/2., (z(0)+z(1))/2., e, RPNi);
- mid_[1] = new DI_Point((x(1)+x(2))/2., (y(1)+y(2))/2., (z(1)+z(2))/2., e, RPNi);
- mid_[2] = new DI_Point((x(2)+x(3))/2., (y(2)+y(3))/2., (z(2)+z(3))/2., e, RPNi);
- mid_[3] = new DI_Point((x(3)+x(0))/2., (y(3)+y(0))/2., (z(3)+z(0))/2., e, RPNi);
- mid_[4] = new DI_Point((x(0)+x(4))/2., (y(0)+y(4))/2., (z(0)+z(4))/2., e, RPNi);
- mid_[5] = new DI_Point((x(1)+x(5))/2., (y(1)+y(5))/2., (z(1)+z(5))/2., e, RPNi);
- mid_[6] = new DI_Point((x(2)+x(6))/2., (y(2)+y(6))/2., (z(2)+z(6))/2., e, RPNi);
- mid_[7] = new DI_Point((x(3)+x(7))/2., (y(3)+y(7))/2., (z(3)+z(7))/2., e, RPNi);
- mid_[8] = new DI_Point((x(4)+x(5))/2., (y(4)+y(5))/2., (z(4)+z(5))/2., e, RPNi);
- mid_[9] = new DI_Point((x(5)+x(6))/2., (y(5)+y(6))/2., (z(5)+z(6))/2., e, RPNi);
- mid_[10] = new DI_Point((x(6)+x(7))/2., (y(6)+y(7))/2., (z(6)+z(7))/2., e, RPNi);
- mid_[11] = new DI_Point((x(7)+x(4))/2., (y(7)+y(4))/2., (z(7)+z(4))/2., e, RPNi);
- mid_[12] = new DI_Point((x(0)+x(1)+x(2)+x(3))/4., (y(0)+y(1)+y(2)+y(3))/4., (z(0)+z(1)+z(2)+z(3))/4., e, RPNi);
- mid_[13] = new DI_Point((x(0)+x(1)+x(4)+x(5))/4., (y(0)+y(1)+y(4)+y(5))/4., (z(0)+z(1)+z(4)+z(5))/4., e, RPNi);
- mid_[14] = new DI_Point((x(1)+x(2)+x(5)+x(6))/4., (y(1)+y(2)+y(5)+y(6))/4., (z(1)+z(2)+z(5)+z(6))/4., e, RPNi);
- mid_[15] = new DI_Point((x(2)+x(3)+x(6)+x(7))/4., (y(2)+y(3)+y(6)+y(7))/4., (z(2)+z(3)+z(6)+z(7))/4., e, RPNi);
- mid_[16] = new DI_Point((x(0)+x(3)+x(4)+x(7))/4., (y(0)+y(3)+y(4)+y(7))/4., (z(0)+z(3)+z(4)+z(7))/4., e, RPNi);
- mid_[17] = new DI_Point((x(4)+x(5)+x(6)+x(7))/4., (y(4)+y(5)+y(6)+y(7))/4., (z(4)+z(5)+z(6)+z(7))/4., e, RPNi);
- mid_[18] = new DI_Point((x(0)+x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7))/8., (y(0)+y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7))/8.,
- (z(0)+z(1)+z(2)+z(3)+z(4)+z(5)+z(6)+z(7))/8., e, RPNi);
- quad_ = true;
-}*/
// ----------------------------------------------------------------------------
// -------------------------- SPLITTING ---------------------------------------
@@ -1633,7 +1532,8 @@ void DI_Quad::splitIntoTriangles(std::vector<DI_Triangle> &triangles) const
triangles.push_back(DI_Triangle(pt(1), pt(2), pt(3), lsTag())); // 123
}
-// Split a reference tetrahedron cut by a level set (defined in a hex) into sub tetrahedra and create triangles on the surface of the level set
+// Split a reference tetrahedron cut by a level set (defined in a hex) into sub tetrahedra
+// and create triangles on the surface of the level set
void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
std::vector<DI_Tetra> &subTetras, std::vector<DI_Triangle> &surfTriangles,
std::vector<DI_CuttingPoint> &cuttingPoints, std::vector<DI_QualError> &QError) const
@@ -1662,7 +1562,8 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset
int COUNT = 0; // number of edges cut
int IND[4]; // ids of edges cut
- int tab [6][4] = {{0, 1, 2, 3}, {0, 2, 3, 1}, {0, 3, 1, 2}, {1, 2, 0, 3}, {2, 3, 0, 1}, {3, 1, 0, 2}}; // edges nodes and opposite edges nodes
+ // edges nodes and opposite edges nodes
+ int tab [6][4] = {{0, 1, 2, 3}, {0, 2, 3, 1}, {0, 3, 1, 2}, {1, 2, 0, 3}, {2, 3, 0, 1}, {3, 1, 0, 2}};
// compute the intersections between the edges of the tetra and the level set.
for(int i = 0; i < 6; i++){
@@ -2738,7 +2639,8 @@ bool DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
std::vector<const gLevelset *> RPN, RPNi;
Ls.getRPN(RPN);
- DI_Hexa hh(-1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, -1, 1, 1, 8.); // reference hexa
+ // reference hexa
+ DI_Hexa hh(-1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, -1, 1, 1, 8.);
std::vector<DI_Hexa> hh_subHexas;
std::vector<DI_Tetra> hh_subTetras;
std::vector<DI_Quad> hh_surfQuads;
diff --git a/contrib/DiscreteIntegration/recurCut.cpp b/contrib/DiscreteIntegration/recurCut.cpp
index 14cabddd82..e93eb2d90e 100644
--- a/contrib/DiscreteIntegration/recurCut.cpp
+++ b/contrib/DiscreteIntegration/recurCut.cpp
@@ -37,13 +37,13 @@ void recurCut(RecurElement *re, int maxlevel, int level)
RecurElement *re0 = new RecurElement(&DI_Triangle(p1, p12, p13));
recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re;
- RecurElement *re1 = new RecurElement(&DI_Triangle(p2, p12, p23));
+ RecurElement *re1 = new RecurElement(&DI_Triangle(p2, p23, p12));
recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re;
RecurElement *re2 = new RecurElement(&DI_Triangle(p3, p13, p23));
recurCut(re2, maxlevel, level);
re->sub[2] = re2; re2->super = re;
- RecurElement *re3 = new RecurElement(&DI_Triangle(p12, p13, p23));
+ RecurElement *re3 = new RecurElement(&DI_Triangle(p12, p23, p13));
recurCut(re3, maxlevel, level);
re->sub[3] = re3; re3->super = re;
}
@@ -87,22 +87,22 @@ void recurCut(RecurElement *re, int maxlevel, int level)
RecurElement *re0 = new RecurElement(&DI_Tetra(p1, p12, p13, p14));
recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re;
- RecurElement *re1 = new RecurElement(&DI_Tetra(p2, p12, p23, p24));
+ RecurElement *re1 = new RecurElement(&DI_Tetra(p2, p23, p12, p24));
recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re;
RecurElement *re2 = new RecurElement(&DI_Tetra(p3, p13, p23, p34));
recurCut(re2, maxlevel, level);
re->sub[2] = re2; re2->super = re;
- RecurElement *re3 = new RecurElement(&DI_Tetra(p4, p14, p24, p34));
+ RecurElement *re3 = new RecurElement(&DI_Tetra(p4, p14, p34, p24));
recurCut(re3, maxlevel, level);
re->sub[3] = re3; re3->super = re;
RecurElement *re4 = new RecurElement(&DI_Tetra(p12, p14, p24, p34));
recurCut(re4, maxlevel, level);
re->sub[4] = re4; re4->super = re;
- RecurElement *re5 = new RecurElement(&DI_Tetra(p12, p23, p24, p34));
+ RecurElement *re5 = new RecurElement(&DI_Tetra(p12, p23, p34, p24));
recurCut(re5, maxlevel, level);
re->sub[5] = re5; re5->super = re;
- RecurElement *re6 = new RecurElement(&DI_Tetra(p12, p13, p23, p34));
+ RecurElement *re6 = new RecurElement(&DI_Tetra(p12, p13, p34, p23));
recurCut(re6, maxlevel, level);
re->sub[6] = re6; re6->super = re;
RecurElement *re7 = new RecurElement(&DI_Tetra(p12, p13, p14, p34));
@@ -140,25 +140,25 @@ void recurCut(RecurElement *re, int maxlevel, int level)
RecurElement *re0 = new RecurElement(&DI_Hexa(p1, p12, p1234, p14, p15, p1256, p12345678, p1458));
recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re;
- RecurElement *re1 = new RecurElement(&DI_Hexa(p2, p23, p1234, p12, p26, p2367, p12345678, p1256));
+ RecurElement *re1 = new RecurElement(&DI_Hexa(p12, p2, p23, p1234, p1256, p26, p2367, p12345678));
recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re;
- RecurElement *re2 = new RecurElement(&DI_Hexa(p3, p34, p1234, p23, p37, p3478, p12345678, p2367));
+ RecurElement *re2 = new RecurElement(&DI_Hexa(p1234, p23, p3, p34, p12345678, p2367, p37, p3478));
recurCut(re2, maxlevel, level);
re->sub[2] = re2; re2->super = re;
- RecurElement *re3 = new RecurElement(&DI_Hexa(p4, p14, p1234, p34, p48, p1458, p12345678, p3478));
+ RecurElement *re3 = new RecurElement(&DI_Hexa(p14, p1234, p34, p4, p1458, p12345678, p3478, p48));
recurCut(re3, maxlevel, level);
re->sub[3] = re3; re3->super = re;
- RecurElement *re4 = new RecurElement(&DI_Hexa(p5, p58, p5678, p56, p15, p1458, p12345678, p1256));
+ RecurElement *re4 = new RecurElement(&DI_Hexa(p15, p1256, p12345678, p1458, p5, p56, p5678, p58));
recurCut(re4, maxlevel, level);
re->sub[4] = re4; re4->super = re;
- RecurElement *re5 = new RecurElement(&DI_Hexa(p6, p56, p5678, p67, p26, p1256, p12345678, p2367));
+ RecurElement *re5 = new RecurElement(&DI_Hexa(p1256, p26, p2367, p12345678, p56, p6, p67, p5678));
recurCut(re5, maxlevel, level);
re->sub[5] = re5; re5->super = re;
- RecurElement *re6 = new RecurElement(&DI_Hexa(p7, p67, p5678, p78, p37, p2367, p12345678, p3478));
+ RecurElement *re6 = new RecurElement(&DI_Hexa(p12345678, p2367, p37, p3478, p5678, p67, p7, p78));
recurCut(re6, maxlevel, level);
re->sub[6] = re6; re6->super = re;
- RecurElement *re7 = new RecurElement(&DI_Hexa(p8, p78, p5678, p58, p48, p3478, p12345678, p1458));
+ RecurElement *re7 = new RecurElement(&DI_Hexa(p1458, p12345678, p3478, p48, p58, p5678, p78, p8));
recurCut(re7, maxlevel, level);
re->sub[7] = re7; re7->super = re;
}
@@ -186,7 +186,7 @@ bool signChange (RecurElement *re, const DI_Element *e, const std::vector<const
else {
for(unsigned int p = 0; p < cp.size(); p++)
cp[p].chooseLs(Lsi);
- re->el->chooseLs(Lsi);
+ if (re->super) re->el->chooseLs(Lsi);
}
}
re->el->clearLs();
diff --git a/utils/api_demos/Makefile b/utils/api_demos/Makefile
index 0c1c04bbe8..5449a0e003 100644
--- a/utils/api_demos/Makefile
+++ b/utils/api_demos/Makefile
@@ -7,7 +7,7 @@ OCC = /Users/remacle/SOURCES/opencascade/lib/libTKSTEP.a /Users/remacle/SOURCES/
mainElasticity: mainElasticity.cpp
g++ ${FLAGS} -o mainElasticity ${INC} mainElasticity.cpp\
- ${GMSH} -llapack -lblas -lm -L/Users/remacle/SOURCES/gmsh/contrib/taucs/lib -ltaucs
+ ${GMSH} -llapack -lblas -lm -lmetis -L/home/gaetan/develop/gmsh/contrib/taucs/lib/linux -ltaucs
mainCartesian: mainCartesian.cpp
g++ ${FLAGS} -o mainCartesian ${INC} mainCartesian.cpp\
@@ -31,7 +31,7 @@ mainPost: mainPost.cpp
mainLevelset: mainLevelset.cpp
g++ ${FLAGS} -o mainLevelset ${INC} mainLevelset.cpp\
- ${GMSH} -llapack -lblas
+ ${GMSH} -llapack -lblas -lm -lmetis -L/home/gaetan/develop/gmsh/contrib/taucs/lib/linux -ltaucs
clean:
rm -rf mainSimple mainAntTweakBar mainGlut mainPost mainLevelset\
diff --git a/utils/api_demos/mainElasticity.cpp b/utils/api_demos/mainElasticity.cpp
index 543a73b312..03fb945db9 100644
--- a/utils/api_demos/mainElasticity.cpp
+++ b/utils/api_demos/mainElasticity.cpp
@@ -12,6 +12,7 @@ int main (int argc, char* argv[]){
// globals are still present in Gmsh
GmshInitialize(argc, argv);
+ GmshSetOption("General","Terminal",1.);
// instanciate a solver
elasticitySolver mySolver (1000);
@@ -20,7 +21,7 @@ int main (int argc, char* argv[]){
mySolver.readInputFile(argv[1]);
// solve the problem
- mySolver.solve();
+ mySolver.solve(); printf("problem solved\n");
PView *pv = mySolver.buildDisplacementView("displacement");
pv->getData()->writeMSH("disp.msh", false);
--
GitLab