diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index bad7035122383a1a7dff22d77084c63cb9485576..da9cb2e04153d8254328fd611905b15b84422830 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -266,7 +266,6 @@ double MElement::getJacobian(double u, double v, double w, double jac[3][3])
jac[j][2] += v->z() * gg[j];
}
}
-
return _computeDeterminantAndRegularize(this, jac);
}
@@ -290,7 +289,6 @@ double MElement::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
jac[j][2] += v->z() * gsf(i, j);
}
}
-
return _computeDeterminantAndRegularize(this, jac);
}
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 2725ef0588daf13d5948c884970b815daf01bea7..34c8afd7e2cbd0a3c511622b003534361fc9d076 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -921,7 +921,7 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
bool own = (eParent && !e->ownsParent()) ? false : true;
if(poly[0].size()) {
if(eType == MSH_TRI_B || eType == MSH_POLYG_B)
- p1 = new MPolygonBorder(poly[0], ++numEle, ePart,
+ p1 = new MPolygonBorder(poly[0], ++numEle, ePart, own, parent,
copy->getDomain(0), copy->getDomain(1));
else p1 = new MPolygon(poly[0], ++numEle, ePart, own, parent);
own = false;
@@ -936,7 +936,7 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
}
if(poly[1].size()) {
if(eType == MSH_TRI_B || eType == MSH_POLYG_B)
- p2 = new MPolygonBorder(poly[1], ++numEle, ePart,
+ p2 = new MPolygonBorder(poly[1], ++numEle, ePart, own, parent,
copy->getDomain(0), copy->getDomain(1));
else p2 = new MPolygon(poly[1], ++numEle, ePart, own, parent);
elements[8][elementary].push_back(p2);
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index e18bb4a8f7f13d589c56be4b6cd5a8d07577093f..147c5715e911bc30e2dbe573325c6a781308833b 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -391,15 +391,15 @@ class MPolygonBorder : public MPolygon {
MElement* _domains[2];
IntPt *_intpt;
public:
- MPolygonBorder(std::vector<MTriangle*> v, int num = 0, int part = 0,
- MElement* d1 = NULL, MElement* d2 = NULL)
- : MPolygon(v, num, part), _intpt(0)
+ MPolygonBorder(std::vector<MTriangle*> v, int num = 0, int part = 0, bool own = false,
+ MElement *p = NULL, MElement *d1 = NULL, MElement *d2 = NULL)
+ : MPolygon(v, num, part, own, p), _intpt(0)
{
_domains[0] = d1; _domains[1] = d2;
}
- MPolygonBorder(std::vector<MVertex*> v, int num = 0, int part = 0,
- MElement* d1 = NULL, MElement* d2 = NULL)
- : MPolygon(v, num, part), _intpt(0)
+ MPolygonBorder(std::vector<MVertex*> v, int num = 0, int part = 0, bool own = false,
+ MElement *p = NULL, MElement* d1 = NULL, MElement* d2 = NULL)
+ : MPolygon(v, num, part, own, p), _intpt(0)
{
_domains[0] = d1; _domains[1] = d2;
}
diff --git a/Post/adaptiveData.h b/Post/adaptiveData.h
index 1b2edf5e9df464d81b074398b0d52eb14b293dbb..e4479cb786c26cc9d3c5f1f91cc2121bde131284 100644
--- a/Post/adaptiveData.h
+++ b/Post/adaptiveData.h
@@ -137,14 +137,13 @@ class adaptivePrism {
public:
bool visible;
adaptivePoint *p[6];
- adaptivePrism *e[12];
+ adaptivePrism *e[8];
static std::list<adaptivePrism*> all;
static std::set<adaptivePoint> allPoints;
static int numNodes, numEdges;
public:
- adaptivePrism(adaptivePoint *p1, adaptivePoint *p2,
- adaptivePoint *p3, adaptivePoint *p4,
- adaptivePoint *p5, adaptivePoint *p6)
+ adaptivePrism(adaptivePoint *p1, adaptivePoint *p2, adaptivePoint *p3,
+ adaptivePoint *p4, adaptivePoint *p5, adaptivePoint *p6)
: visible(false)
{
p[0] = p1;
diff --git a/Solver/FuncGradDisc.h b/Solver/FuncGradDisc.h
index a508c3aad250dadb2fa685120b41a1e032b0ef2c..7021490625d2a166ce1c197c4d778198c30434f8 100644
--- a/Solver/FuncGradDisc.h
+++ b/Solver/FuncGradDisc.h
@@ -20,46 +20,40 @@
class FuncGradDisc : public simpleFunctionOnElement<double> {
private :
-
- gLevelset *_ls;
- GModel * _pModel;
-
+ gLevelset *_ls;
+ GModel * _pModel;
public :
- FuncGradDisc(gLevelset *ls,GModel *pModel){
- _ls = ls;
- _pModel = pModel;
-
+ FuncGradDisc(gLevelset *ls, GModel *pModel){
+ _ls = ls;
+ _pModel = pModel;
}
- double operator () (double x,double y,double z) const {
-
+ double operator () (double x, double y, double z) const {
- // --- F2 --- //
- MElement *e = getElement();
- SPoint3 p(x,y,z);
+ // --- F2 --- //
+ MElement *e = getElement();
+ SPoint3 p(x, y, z);
- if (e->getParent()) e = e->getParent();
- double xyz[3] = {x,y,z};
- double uvw[3];
- e->xyz2uvw(xyz,uvw);
- double val[30];
- e->getShapeFunctions(uvw[0], uvw[1], uvw[2], val);
- double f = 0;
- for (int i = 0; i < e->getNumShapeFunctions(); i++)
- {
- MVertex *v = e->getShapeFunctionNode(i);
- //std::cout<<"val[i] :" << val[i] << "\n";
- //std::cout<<"ls(i) :" << (*_ls)(v->x(),v->y(),v->z()) << "\n";
- f = f + std::abs((*_ls)(v->x(), v->y(), v->z())) * val[i];
- }
- f = f - std::abs((*_ls)(x, y, z));
+ if (e->getParent()) e = e->getParent();
+ double xyz[3] = {x, y, z};
+ double uvw[3];
+ e->xyz2uvw(xyz, uvw);
+ double val[30];
+ e->getShapeFunctions(uvw[0], uvw[1], uvw[2], val);
+ double f = 0;
+ for (int i = 0; i < e->getNumShapeFunctions(); i++){
+ MVertex *v = e->getShapeFunctionNode(i);
+ //std::cout<<"val[i] :" << val[i] << "\n";
+ //std::cout<<"ls(i) :" << (*_ls)(v->x(),v->y(),v->z()) << "\n";
+ f = f + std::abs((*_ls)(v->x(), v->y(), v->z())) * val[i];
+ }
+ f = f - std::abs((*_ls)(x, y, z));
- //std::cout<<"val f :" << f << "\n";
- return f;
+ //std::cout<<"val f :" << f << "\n";
+ return f;
-
- // --- F1 --- //
+ // --- F1 --- //
// SPoint3 p(x,y,z);
@@ -83,66 +77,61 @@ class FuncGradDisc : public simpleFunctionOnElement<double> {
void gradient (double x, double y, double z,
- double & dfdx, double & dfdy, double & dfdz) const
- {
-
-
- // ---- F2 ---- //
- MElement *e=getElement();
- SPoint3 p(x,y,z);
- if (e->getParent()) e = e->getParent();
- double xyz[3] = {x,y,z};
- double uvw[3];
- e->xyz2uvw(xyz,uvw);
- double gradsuvw[256][3];
- e->getGradShapeFunctions(uvw[0],uvw[1],uvw[2],gradsuvw);
-
- double jac[3][3];
- double invjac[3][3];
- double dNdx;
- double dNdy;
- double dNdz;
- const double detJ = e->getJacobian(uvw[0], uvw[1], uvw[2], jac);
- inv3x3(jac, invjac);
-
- dfdx = 0;
- dfdy = 0;
- dfdz = 0;
-
- if ((*_ls)(x,y,z)>0)
- {
- for (int i = 0; i < e->getNumShapeFunctions(); i++)
- {
- dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
- dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
- dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
-
- MVertex *v = e->getShapeFunctionNode(i);
- dfdx = dfdx + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdx ;
- dfdx = dfdx - (*_ls)(v->x(), v->y(), v->z()) * dNdx;
- dfdy = dfdy + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdy ;
- dfdy = dfdy - (*_ls)(v->x(), v->y(), v->z()) * dNdy;
- dfdz = dfdz + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdz ;
- dfdz = dfdz - (*_ls)(v->x(), v->y(), v->z()) * dNdz;
- }
- }else
- {
- for (int i = 0; i < e->getNumShapeFunctions(); i++)
- {
- dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
- dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
- dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
-
- MVertex *v = e->getShapeFunctionNode(i);
- dfdx = dfdx + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdx ;
- dfdx = dfdx + (*_ls)(v->x(), v->y(), v->z()) * dNdx;
- dfdy = dfdy + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdy ;
- dfdy = dfdy + (*_ls)(v->x(), v->y(), v->z()) * dNdy;
- dfdz = dfdz + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdz ;
- dfdz = dfdz + (*_ls)(v->x(), v->y(), v->z()) * dNdz;
- }
- }
- }
+ double &dfdx, double &dfdy, double &dfdz) const {
+
+ // ---- F2 ---- //
+ MElement *e=getElement();
+ SPoint3 p(x, y, z);
+ if (e->getParent()) e = e->getParent();
+ double xyz[3] = {x, y, z};
+ double uvw[3];
+ e->xyz2uvw(xyz, uvw);
+ double gradsuvw[256][3];
+ e->getGradShapeFunctions(uvw[0], uvw[1], uvw[2], gradsuvw);
+
+ double jac[3][3];
+ double invjac[3][3];
+ double dNdx;
+ double dNdy;
+ double dNdz;
+ const double detJ = e->getJacobian(uvw[0], uvw[1], uvw[2], jac);
+ inv3x3(jac, invjac);
+
+ dfdx = 0;
+ dfdy = 0;
+ dfdz = 0;
+
+ if ((*_ls)(x, y, z) > 0) {
+ for (int i = 0; i < e->getNumShapeFunctions(); i++) {
+ dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
+ dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
+ dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
+
+ MVertex *v = e->getShapeFunctionNode(i);
+ dfdx = dfdx + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdx ;
+ dfdx = dfdx - (*_ls)(v->x(), v->y(), v->z()) * dNdx;
+ dfdy = dfdy + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdy ;
+ dfdy = dfdy - (*_ls)(v->x(), v->y(), v->z()) * dNdy;
+ dfdz = dfdz + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdz ;
+ dfdz = dfdz - (*_ls)(v->x(), v->y(), v->z()) * dNdz;
+ }
+ }
+ else{
+ for (int i = 0; i < e->getNumShapeFunctions(); i++) {
+ dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
+ dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
+ dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
+
+ MVertex *v = e->getShapeFunctionNode(i);
+ dfdx = dfdx + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdx ;
+ dfdx = dfdx + (*_ls)(v->x(), v->y(), v->z()) * dNdx;
+ dfdy = dfdy + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdy ;
+ dfdy = dfdy + (*_ls)(v->x(), v->y(), v->z()) * dNdy;
+ dfdz = dfdz + std::abs((*_ls)(v->x(), v->y(), v->z())) * dNdz ;
+ dfdz = dfdz + (*_ls)(v->x(), v->y(), v->z()) * dNdz;
+ }
+ }
+ }
// ---- F1 ------ //
@@ -197,7 +186,6 @@ class FuncGradDisc : public simpleFunctionOnElement<double> {
// }
// }
-
};
#endif
diff --git a/Solver/FuncHeaviside.h b/Solver/FuncHeaviside.h
index d576b4cf8eb8e1eef5b37acd086cb95c04b192f6..6a24f67f2db2a18604c2b3a2e00b4629e354ab83 100644
--- a/Solver/FuncHeaviside.h
+++ b/Solver/FuncHeaviside.h
@@ -15,34 +15,30 @@
#include "DILevelset.h"
-class FuncHeaviside : public simpleFunctionOnElement<double> {
+class FuncHeaviside : public simpleFunctionOnElement<double> {
private :
gLevelset *_ls;
public :
- FuncHeaviside(gLevelset *ls){
- _ls = ls;
+ FuncHeaviside(gLevelset *ls) : _ls(ls) { }
+ virtual double operator() (double x, double y, double z) const {
+ if (_ls->isInsideDomain (x, y, z))
+ return 1;
+ else
+ return -1;
}
- virtual double operator () (double x,double y,double z) const {
- if (_ls->isInsideDomain (x,y,z)){
- return 1 ;
- }else{
- return -1;
- }
- }
- virtual double operator () (double x,double y,double z, MElement *e) const {
- if (_ls->isInsideDomain (x,y,z)){
- return 1 ;
- }else{
- return -1;
- }
+ virtual double operator() (double x,double y,double z, MElement *e) const {
+ if (_ls->isInsideDomain (x, y, z))
+ return 1;
+ else
+ return -1;
}
virtual void gradient (double x, double y, double z,
- double & dfdx, double & dfdy, double & dfdz) const
+ double &dfdx, double &dfdy, double &dfdz) const
{ dfdx = dfdy = dfdz = 0.0; }
virtual void gradient (double x, double y, double z,
- double & dfdx, double & dfdy, double & dfdz, MElement *e) const
+ double &dfdx, double &dfdy, double &dfdz, MElement *e) const
{ dfdx = dfdy = dfdz = 0.0; }
};
diff --git a/Solver/convexCombinationTerm.h b/Solver/convexCombinationTerm.h
index b2e5b79e5c13b2747f1bcaaabb03a45bc1739349..d999a2182c1e894c1ef23ba213775c8127d322ff 100644
--- a/Solver/convexCombinationTerm.h
+++ b/Solver/convexCombinationTerm.h
@@ -45,9 +45,9 @@ class convexCombinationTerm : public femTerm<double> {
const double _diff = 1.0;
for (int j = 0; j < nbSF; j++){
for (int k = 0; k < nbSF; k++){
- m(j,k) = -1.*_diff;
+ m(j, k) = -1. * _diff;
}
- m(j,j) = (nbSF - 1) * _diff;
+ m(j, j) = (nbSF - 1) * _diff;
}
}
diff --git a/Solver/crossConfTerm.h b/Solver/crossConfTerm.h
index faffd4c086e1b81da8bc2c3ffcea8a40324043df..2f4dab63b0ffdf99082158af49d49da45e61b8d4 100644
--- a/Solver/crossConfTerm.h
+++ b/Solver/crossConfTerm.h
@@ -95,7 +95,7 @@ class crossConfTerm : public femTerm<double> {
int nbSF = e->getNumShapeFunctions();
fullMatrix<double> *mat;
- mat = new fullMatrix<double>(nbSF,nbSF);
+ mat = new fullMatrix<double>(nbSF, nbSF);
elementMatrix(se, *mat);
fullVector<double> val(nbSF);
@@ -110,7 +110,7 @@ class crossConfTerm : public femTerm<double> {
m.scale(0.);
for (int i = 0; i < nbSF; i++)
for (int j = 0; j < nbSF; j++)
- m(i) += -(*mat)(i,j)*val(j);
+ m(i) += -(*mat)(i, j) * val(j);
}
};
diff --git a/Solver/diagBCTerm.h b/Solver/diagBCTerm.h
index b7bc918e4893e8fd3aa071f5d051b89c3847c830..37bd9d3c6390c363de4ed1bfa7eaf343313560bf 100644
--- a/Solver/diagBCTerm.h
+++ b/Solver/diagBCTerm.h
@@ -43,12 +43,12 @@ class diagBCTerm : public femTerm<double> {
const int nbSF = e->getNumShapeFunctions();
for (int j = 0; j < nbSF; j++){
for (int k = 0; k < nbSF; k++) {
- m(j,k) = 0.0;
- m(k,j) = 0.0;
+ m(j, k) = 0.0;
+ m(k, j) = 0.0;
}
MVertex *v = e->getShapeFunctionNode(j);
- if( v->onWhat()->dim() < 2 ) m(j,j) = 1.0;
- else m(j,j) = 0.0;
+ if( v->onWhat()->dim() < 2 ) m(j, j) = 1.0;
+ else m(j, j) = 0.0;
}
}
};
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 61ae399f5dee82479208ca9d71ec9d3219490fdd..4d1530e07c1c8a8fd3ef0a6b1b247020653d49aa 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -198,12 +198,12 @@ class dofManager{
{
numberDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
}
- virtual inline void getDofValue(std::vector<Dof> &keys,std::vector<dataVec> &Vals)
+ virtual inline void getDofValue(std::vector<Dof> &keys, std::vector<dataVec> &Vals)
{
- int ndofs=keys.size();
+ int ndofs = keys.size();
size_t originalSize = Vals.size();
- Vals.resize(originalSize+ndofs);
- for (int i=0;i<ndofs;++i) getDofValue(keys[i], Vals[originalSize+i]);
+ Vals.resize(originalSize + ndofs);
+ for (int i = 0; i < ndofs; ++i) getDofValue(keys[i], Vals[originalSize+i]);
}
virtual inline void getDofValue(Dof key, dataVec &val) const
{
@@ -211,7 +211,7 @@ class dofManager{
typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
if (it != fixed.end()) {
val = it->second;
- return ;
+ return;
}
}
{
@@ -227,11 +227,11 @@ class dofManager{
if (it != constraints.end()){
dataVec tmp(val);
val = it->second.shift;
- for (unsigned i=0;i<(it->second).linear.size();i++){
+ for (unsigned i = 0; i < (it->second).linear.size(); i++){
std::map<Dof, int>::const_iterator itu = unknown.find
(((it->second).linear[i]).first);
getDofValue(((it->second).linear[i]).first, tmp);
- dofTraits<T>::gemm(val,((it->second).linear[i]).second, tmp, 1, 1);
+ dofTraits<T>::gemm(val, ((it->second).linear[i]).second, tmp, 1, 1);
}
return ;
}
@@ -260,13 +260,13 @@ class dofManager{
if (it != constraints.end()){
v = it->second.shift;
dataVec tmp(v);
- for (unsigned i=0;i<(it->second).linear.size();i++){
+ for (unsigned i = 0; i < (it->second).linear.size(); i++){
std::map<Dof, int>::const_iterator itu = unknown.find
(((it->second).linear[i]).first);
getDofValue(((it->second).linear[i]).first, tmp);
dofTraits<T>::gemm(v, ((it->second).linear[i]).second, tmp, 1, 1);
}
- return ;
+ return;
}
}
}
@@ -280,16 +280,16 @@ class dofManager{
std::map<Dof, int>::iterator itR = unknown.find(R);
if (itR != unknown.end())
{
- typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator itConstraint;
itConstraint = constraints.find(C);
if (itConstraint != constraints.end()){
for (unsigned i = 0; i < (itConstraint->second).linear.size(); i++){
- insertInSparsityPattern(R,(itConstraint->second).linear[i].first);
+ insertInSparsityPattern(R, (itConstraint->second).linear[i].first);
}
}
}
else{ // test function ; (no shift ?)
- typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator itConstraint;
itConstraint = constraints.find(R);
if (itConstraint != constraints.end()){
for (unsigned i = 0; i < (itConstraint->second).linear.size(); i++){
@@ -310,9 +310,10 @@ class dofManager{
_current->insertInSparsityPattern(itR->second, itC->second);
}
else{
- typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
+ typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
if (itFixed != fixed.end()) {
- }else insertInSparsityPatternLinConst(R, C);
+ }
+ else insertInSparsityPatternLinConst(R, C);
}
}
if (itR == unknown.end())
@@ -332,13 +333,14 @@ class dofManager{
_current->addToMatrix(itR->second, itC->second, value);
}
else{
- typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
+ typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C);
if (itFixed != fixed.end()) {
// tmp = -value * itFixed->second
dataVec tmp(itFixed->second);
- dofTraits<T>::gemm(tmp , value, itFixed->second, -1, 0);
+ dofTraits<T>::gemm(tmp, value, itFixed->second, -1, 0);
_current->addToRightHandSide(itR->second, tmp);
- }else assembleLinConst(R, C, value);
+ }
+ else assembleLinConst(R, C, value);
}
}
if (itR == unknown.end())
@@ -375,9 +377,10 @@ class dofManager{
if (itFixed != fixed.end()){
// tmp = -m(i,j) * itFixed->second
dataVec tmp(itFixed->second);
- dofTraits<T>::gemm(tmp , m(i,j), itFixed->second, -1, 0);
+ dofTraits<T>::gemm(tmp, m(i, j), itFixed->second, -1, 0);
_current->addToRightHandSide(NR[i], tmp);
- }else assembleLinConst(R[i], C[j], m(i,j));
+ }
+ else assembleLinConst(R[i], C[j], m(i, j));
}
}
}
@@ -404,13 +407,13 @@ class dofManager{
_current->addToRightHandSide(NR[i], m(i));
}
else{
- typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator itConstraint;
itConstraint = constraints.find(R[i]);
if (itConstraint != constraints.end()){
for (unsigned j = 0; j < (itConstraint->second).linear.size(); j++){
dataMat tmp;
- dofTraits<T>::gemm(tmp,(itConstraint->second).linear[j].second,m(i), 1, 0);
- assemble((itConstraint->second).linear[j].first,tmp);
+ dofTraits<T>::gemm(tmp, (itConstraint->second).linear[j].second, m(i), 1, 0);
+ assemble((itConstraint->second).linear[j].first, tmp);
}
}
}
@@ -437,15 +440,15 @@ class dofManager{
if (itFixed != fixed.end()){
// tmp = -m(i,j) * itFixed->second
dataVec tmp(itFixed->second);
- dofTraits<T>::gemm(tmp , m(i,j), itFixed->second, -1, 0);
+ dofTraits<T>::gemm(tmp, m(i, j), itFixed->second, -1, 0);
_current->addToRightHandSide(NR[i], tmp);
- } else assembleLinConst(R[i],R[j],m(i,j));
+ } else assembleLinConst(R[i], R[j], m(i, j));
}
}
}
else{
for (unsigned int j = 0; j < R.size(); j++){
- assembleLinConst(R[i],R[j],m(i,j));
+ assembleLinConst(R[i], R[j], m(i, j));
}
}
}
@@ -471,13 +474,13 @@ class dofManager{
_current->addToRightHandSide(itR->second, value);
}
else{
- typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator itConstraint;
itConstraint = constraints.find(R);
if (itConstraint != constraints.end()){
for (unsigned j = 0; j < (itConstraint->second).linear.size(); j++){
dataMat tmp;
- dofTraits<T>::gemm(tmp,(itConstraint->second).linear[j].second,value, 1, 0);
- assemble((itConstraint->second).linear[j].first,tmp);
+ dofTraits<T>::gemm(tmp, (itConstraint->second).linear[j].second, value, 1, 0);
+ assemble((itConstraint->second).linear[j].first, tmp);
}
}
}
@@ -506,7 +509,7 @@ class dofManager{
if(it != _linearSystems.end())
_current = it->second;
else{
- Msg::Error("Current matrix %s not found ", name.c_str());
+ Msg::Error("Current matrix %s not found ", name.c_str());
throw;
}
}
@@ -522,20 +525,20 @@ class dofManager{
inline void setLinearConstraint (Dof key, DofAffineConstraint<dataVec> &affineconstraint)
{
constraints[key] = affineconstraint;
- // constraints.insert(std::make_pair(key,affineconstraint));
+ // constraints.insert(std::make_pair(key, affineconstraint));
}
inline void assembleLinConst(const Dof &R, const Dof &C, const dataMat &value)
{
std::map<Dof, int>::iterator itR = unknown.find(R);
if (itR != unknown.end())
{
- typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator itConstraint;
itConstraint = constraints.find(C);
if (itConstraint != constraints.end()){
dataMat tmp(value);
for (unsigned i = 0; i < (itConstraint->second).linear.size(); i++){
dofTraits<T>::gemm(tmp, (itConstraint->second).linear[i].second, value, 1, 0);
- assemble(R,(itConstraint->second).linear[i].first,tmp);
+ assemble(R, (itConstraint->second).linear[i].first, tmp);
}
dataMat tmp2(value);
dofTraits<T>::gemm(tmp2, value, itConstraint->second.shift, -1, 0);
@@ -543,7 +546,7 @@ class dofManager{
}
}
else{ // test function ; (no shift ?)
- typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ typename std::map<Dof, DofAffineConstraint<dataVec> >::iterator itConstraint;
itConstraint = constraints.find(R);
if (itConstraint != constraints.end()){
dataMat tmp(value);
@@ -591,7 +594,7 @@ void dofManager<T>::_parallelFinalize()
nRequest[i] = 0;
for (std::map <Dof, int >::iterator it = ghosted.begin(); it != ghosted.end(); it++)
nRequest[it->second] ++;
- MPI_Alltoall(nRequest,1,MPI_INT,nRequested,1,MPI_INT,MPI_COMM_WORLD);
+ MPI_Alltoall(nRequest, 1, MPI_INT, nRequested, 1, MPI_INT, MPI_COMM_WORLD);
long int **recv0 = new long int*[Msg::GetCommSize()];
int **recv1 = new int*[Msg::GetCommSize()];
long int **send0 = new long int*[Msg::GetCommSize()];
@@ -629,7 +632,7 @@ void dofManager<T>::_parallelFinalize()
index != MPI_UNDEFINED) {
if (status.MPI_TAG == 0) {
for (int j = 0; j < nRequested[index]; j++) {
- std::map<Dof,int>::iterator it = unknown.find
+ std::map<Dof, int>::iterator it = unknown.find
(Dof(recv0[index][j*2], recv0[index][j*2+1]));
if (it == unknown.end ())
Msg::Error ("ghost Dof does not exist on parent process");
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 249bef5b8907d75ef55d05f8adb86515cb1be3e9..0866a5c05b853ecb791622b9818194d6f9d8171e 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -36,24 +36,24 @@ struct elasticField {
struct BoundaryCondition
{
- int _tag; // tag for the dofManager
- enum location{UNDEF,ON_VERTEX,ON_EDGE,ON_FACE,ON_VOLUME};
+ int _tag; // tag for the dofManager
+ enum location{UNDEF, ON_VERTEX, ON_EDGE, ON_FACE, ON_VOLUME};
location onWhat; // on vertices or elements
groupOfElements *g; // support for this BC
- BoundaryCondition() : _tag(0),onWhat(UNDEF),g(0) {}
+ BoundaryCondition() : _tag(0), onWhat(UNDEF), g(0) {}
};
struct dirichletBC : public BoundaryCondition
{
int _comp; // component
simpleFunction<double> *_f;
- dirichletBC ():BoundaryCondition(),_comp(0),_f(0){}
+ dirichletBC ():BoundaryCondition(), _comp(0), _f(0){}
};
struct neumannBC : public BoundaryCondition
{
simpleFunction<SVector3> *_f;
- neumannBC () : BoundaryCondition(),_f(NULL){}
+ neumannBC () : BoundaryCondition(), _f(NULL){}
};
// an elastic solver ...
class elasticitySolver
@@ -74,7 +74,7 @@ class elasticitySolver
std::vector<dirichletBC> allDirichlet;
public:
- elasticitySolver(int tag) : _tag(tag),pAssembler(0),LagSpace(0),LagrangeMultiplierSpace(0) {}
+ elasticitySolver(int tag) : _tag(tag), pAssembler(0), LagSpace(0), LagrangeMultiplierSpace(0) {}
elasticitySolver(GModel *model, int tag);
@@ -82,12 +82,12 @@ class elasticitySolver
void addNeumannBC (int dim, int entityId, const std::vector<double> value);
void addElasticDomain (int tag, double e, double nu);
- #if defined (HAVE_LUA)
+#if defined (HAVE_LUA)
void addDirichletBCLua (int dim, int entityId, int component, std::string luaFunctionName, lua_State *L);
void addNeumannBCLua (int dim, int entityId, std::string luaFunctionName, lua_State *L);
- #endif
+#endif
virtual ~elasticitySolver()
{
diff --git a/Solver/filters.h b/Solver/filters.h
index ea322ead152364db935f7d132bd4c2ef94c0b889..5e05b1745dddacbde373b66a0f2b2de6381f50a5 100644
--- a/Solver/filters.h
+++ b/Solver/filters.h
@@ -19,35 +19,27 @@
class FilterNodeEnriched
{
-
private :
-
std::set<int> *_TagEnrichedVertex;
std::set<int> * _EnrichComp;
public :
-
- FilterNodeEnriched(std::set<int > * TagEnrichedVertex , std::set<int> * EnrichComp)
- {
+ FilterNodeEnriched(std::set<int> *TagEnrichedVertex, std::set<int> *EnrichComp){
_TagEnrichedVertex = TagEnrichedVertex;
_EnrichComp = EnrichComp;
}
- virtual bool operator () (Dof & key) const
- {
+ virtual bool operator() (Dof &key) const{
std::set<int>::iterator it1;
std::set<int>::iterator it2;
- int i1,i2;
- Dof::getTwoIntsFromType(key.getType(), i1,i2);
- it2 = _EnrichComp->find(i1);
+ int i1, i2;
+ Dof::getTwoIntsFromType(key.getType(), i1, i2);
+ it2 = _EnrichComp->find(i1);
it1 = _TagEnrichedVertex->find(key.getEntity());
- if ((it1!=_TagEnrichedVertex->end()) && (it2 != _EnrichComp->end()))
- {
+ if ((it1 != _TagEnrichedVertex->end()) && (it2 != _EnrichComp->end()))
return true;
- }
else return false;
}
-
//std::vector<int> * getEnrichComp(){return _EnrichComp;}
// void SetEnrichedVertex(MElement *elep, std::vector<int> & EnrichedVertex,int &nbdofs)
@@ -69,57 +61,40 @@ class FilterNodeEnriched
class FilterElementsCutByLevelSet
{
-
- private :
-
- std::set<int> _TagEnrichedVertex;
- std::pair<int,int> _LevelSetEntity;
- std::set<int> *_EnrichComp;
-
- public :
-
- FilterElementsCutByLevelSet(std::pair<int,int> LevelSetEntity , std::set<int> * EnrichComp)
- {
-
- _EnrichComp = EnrichComp;
- _LevelSetEntity = LevelSetEntity;
-
- // groupOfElements to get all the elements associate with the level set -- (work with *current GModel)
- groupOfElements *LevelSetElements = new groupOfElements (_LevelSetEntity.first, _LevelSetEntity.second);
-
- // tag enriched vertex determination
- std::set<MElement*>::const_iterator it = LevelSetElements->begin();
- for (; it != LevelSetElements->end(); it++)
- {
- MElement *e = *it;
- if (e->getParent()) // if element got parents
- {
- for (int k = 0; k < e->getParent()->getNumVertices(); ++k)
- { // for all vertices in the element parent
- _TagEnrichedVertex.insert(e->getParent()->getVertex(k)->getNum());
- }
- }
- }
-
- }
-
- virtual bool operator () (Dof & key) const
- {
- std::set<int>::const_iterator it1;
- std::set<int>::const_iterator it2;
- int i1,i2;
- Dof::getTwoIntsFromType(key.getType(), i1,i2);
- it2 = _EnrichComp->find(i1);
- it1 = _TagEnrichedVertex.find(key.getEntity());
- if ((it1!=_TagEnrichedVertex.end()) && (it2 != _EnrichComp->end()))
- {
- return true;
+ private :
+ std::set<int> _TagEnrichedVertex;
+ std::pair<int,int> _LevelSetEntity;
+ std::set<int> *_EnrichComp;
+ public :
+ FilterElementsCutByLevelSet(std::pair<int,int> LevelSetEntity , std::set<int> * EnrichComp){
+ _EnrichComp = EnrichComp;
+ _LevelSetEntity = LevelSetEntity;
+ // groupOfElements to get all the elements associate with the level set -- (work with *current GModel)
+ groupOfElements *LevelSetElements = new groupOfElements (_LevelSetEntity.first, _LevelSetEntity.second);
+ // tag enriched vertex determination
+ std::set<MElement*>::const_iterator it = LevelSetElements->begin();
+ for (; it != LevelSetElements->end(); it++) {
+ MElement *e = *it;
+ if (e->getParent()){ // if element got parents
+ for (int k = 0; k < e->getParent()->getNumVertices(); ++k){ // for all vertices in the element parent
+ _TagEnrichedVertex.insert(e->getParent()->getVertex(k)->getNum());
+ }
}
- else return false;
}
-
-
+ }
+
+ virtual bool operator () (Dof & key) const{
+ std::set<int>::const_iterator it1;
+ std::set<int>::const_iterator it2;
+ int i1, i2;
+ Dof::getTwoIntsFromType(key.getType(), i1, i2);
+ it2 = _EnrichComp->find(i1);
+ it1 = _TagEnrichedVertex.find(key.getEntity());
+ if ((it1!=_TagEnrichedVertex.end()) && (it2 != _EnrichComp->end())){
+ return true;
+ }
+ else return false;
+ }
};
-
#endif
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 76b621281d80f48707de24d297edef7fb7189e35..0fd8d52fa380a7d28346b9bec6c56e21e019d0b6 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -23,13 +23,13 @@ void function::addFunctionReplace(functionReplace &fr)
void function::setArgument(fullMatrix<double> &v, const function *f, int iMap)
{
- if (f==NULL)
+ if (f == NULL)
throw;
arguments.push_back(argument(v, iMap, f));
dependencies.insert(dependency(iMap, f));
for (std::set<dependency>::const_iterator it = f->dependencies.begin();
it != f->dependencies.end(); it++) {
- if (it->iMap> 0 && iMap >0)
+ if (it->iMap > 0 && iMap > 0)
Msg::Error("Consecutive secondary caches");
dependencies.insert(dependency(iMap + it->iMap, it->f));
}
@@ -37,7 +37,7 @@ void function::setArgument(fullMatrix<double> &v, const function *f, int iMap)
functionReplace &replace = *_functionReplaces[i];
for (std::set<dependency>::iterator it = replace._fromParent.begin();
it != replace._fromParent.end(); it++) {
- if (it->iMap> 0 && iMap >0)
+ if (it->iMap > 0 && iMap > 0)
Msg::Error("Consecutive secondary caches");
dependencies.insert(dependency(iMap + it->iMap, it->f));
}
@@ -128,7 +128,7 @@ class functionNormals : public function {
static functionNormals *_instance;
// constructor is private only 1 instance can exists, call get to
// access the instance
- functionNormals():function(0){}
+ functionNormals() : function(0){}
public:
void call(dataCacheMap *m, fullMatrix<double> &sol)
{
@@ -155,7 +155,7 @@ function *function::getNormals()
functionReplace::functionReplace()
{
- _nChildren=0;
+ _nChildren = 0;
}
void functionReplace::get(fullMatrix<double> &v, const function *f, int iMap)
@@ -163,24 +163,24 @@ void functionReplace::get(fullMatrix<double> &v, const function *f, int iMap)
bool allDepFromParent = true;
for (std::set<function::dependency>::const_iterator itDep = f->dependencies.begin();
itDep != f->dependencies.end(); itDep++){
- bool depFromParent = (_replaced.count(*itDep)==0);
+ bool depFromParent = (_replaced.count(*itDep) == 0);
for (std::set<function::dependency>::const_iterator itDep2 = itDep->f->dependencies.begin();
itDep2 != itDep->f->dependencies.end() && depFromParent; itDep2++)
- depFromParent &= (_replaced.count(*itDep2)==0);
+ depFromParent &= (_replaced.count(*itDep2) == 0);
if(depFromParent)
_master->dependencies.insert(*itDep);
else
allDepFromParent = false;
}
function::dependency asDep(iMap, f);
- if (allDepFromParent && _replaced.count(asDep)==0)
+ if (allDepFromParent && _replaced.count(asDep) == 0)
_master->dependencies.insert(asDep);
_toCompute.push_back(function::argument(v, iMap, f));
}
void functionReplace::replace(fullMatrix<double> &v, const function *f, int iMap)
{
- _replaced.insert(function::dependency(iMap,f));
+ _replaced.insert(function::dependency(iMap, f));
_toReplace.push_back(function::argument(v, iMap, f));
}
@@ -201,7 +201,7 @@ void functionReplace::compute()
// dataCacheDouble
dataCacheDouble::dataCacheDouble(dataCacheMap *m, function *f):
- _cacheMap(*m),_value(m->getNbEvaluationPoints(),f->getNbCol()), _function(f)
+ _cacheMap(*m), _value(m->getNbEvaluationPoints(), f->getNbCol()), _function(f)
{
m->addDataCacheDouble(this, f->isInvalitedOnElement());
_valid = false;
@@ -215,12 +215,12 @@ void dataCacheDouble::resize(int nrow)
void dataCacheDouble::_eval()
{
- for(unsigned int i=0;i<_directDependencies.size(); i++){
+ for(unsigned int i = 0; i < _directDependencies.size(); i++){
_function->arguments[i].val->setAsProxy(_directDependencies[i]->get());
}
for (unsigned i = 0; i < _function->_functionReplaces.size(); i++) {
_function->_functionReplaces[i]->currentCache = &functionReplaceCaches[i];
- for (unsigned j = 0; j < functionReplaceCaches[i].toReplace.size() ; j++){
+ for (unsigned j = 0; j < functionReplaceCaches[i].toReplace.size(); j++){
_function->_functionReplaces[i]->_toReplace[j].val->setAsProxy
((*functionReplaceCaches[i].toReplace[j])._value);
}
@@ -266,12 +266,12 @@ dataCacheDouble *dataCacheMap::get(const function *f, dataCacheDouble *caller, b
// do I have a cache for this function ?
dataCacheDouble *&r = _cacheDoubleMap[f];
// can I use the cache of my parent ?
- if(_parent && r==NULL) {
+ if(_parent && r == NULL) {
bool okFromParent = true;
for (std::set<function::dependency>::const_iterator it = f->dependencies.begin();
it != f->dependencies.end(); it++) {
if (it->iMap > _parent->_secondaryCaches.size())
- okFromParent=false;
+ okFromParent = false;
dataCacheMap *m = getSecondaryCache(it->iMap);
if (m->_cacheDoubleMap.find(it->f) != m->_cacheDoubleMap.end()) {
okFromParent = false;
@@ -279,10 +279,10 @@ dataCacheDouble *dataCacheMap::get(const function *f, dataCacheDouble *caller, b
}
}
if (okFromParent)
- r = _parent->get (f,caller);
+ r = _parent->get (f, caller);
}
// no cache found, create a new one
- if (r==NULL) {
+ if (r == NULL) {
if (!createIfNotPresent)
return NULL;
r = new dataCacheDouble (this, (function*)(f));
@@ -295,12 +295,12 @@ dataCacheDouble *dataCacheMap::get(const function *f, dataCacheDouble *caller, b
functionReplaceCache replaceCache;
functionReplace *replace = f->_functionReplaces[i];
dataCacheMap *rMap = newChild();
- for (unsigned i = 0; i < _secondaryCaches.size(); i ++)
+ for (unsigned i = 0; i < _secondaryCaches.size(); i++)
rMap->addSecondaryCache (getSecondaryCache(i+1)->newChild());
- for (int i = 0; i < replace->_nChildren; i ++)
+ for (int i = 0; i < replace->_nChildren; i++)
rMap->addSecondaryCache (rMap->newChild());
for (std::vector<function::argument>::iterator it = replace->_toReplace.begin();
- it!= replace->_toReplace.end(); it++ ) {
+ it!= replace->_toReplace.end(); it++) {
dataCacheMap *m = rMap->getSecondaryCache(it->iMap);
dataCacheDouble *s = new dataCacheDouble(m, (function*)_translate(it->f));
m->_cacheDoubleMap[_translate(it->f)] = s;
@@ -347,11 +347,11 @@ void dataCacheMap::setNbEvaluationPoints(int nbEvaluationPoints)
{
_nbEvaluationPoints = nbEvaluationPoints;
for(std::set<dataCacheDouble*>::iterator it = _allDataCaches.begin();
- it!= _allDataCaches.end(); it++){
+ it != _allDataCaches.end(); it++){
(*it)->resize(nbEvaluationPoints);
}
for(std::list<dataCacheMap*>::iterator it = _children.begin();
- it!= _children.end(); it++) {
+ it != _children.end(); it++) {
(*it)->setNbEvaluationPoints(nbEvaluationPoints);
}
}
@@ -364,27 +364,27 @@ void functionConstant::set(double val)
{
if(getNbCol() != 1)
Msg::Error ("set scalar value on a vectorial constant function");
- _source(0,0) = val;
+ _source(0, 0) = val;
}
void functionConstant::call(dataCacheMap *m, fullMatrix<double> &val)
{
- for (int i=0; i<val.size1(); i++)
- for (int j=0; j<_source.size1(); j++)
- val(i,j)=_source(j,0);
+ for (int i = 0; i < val.size1(); i++)
+ for (int j = 0; j < _source.size1(); j++)
+ val(i, j)=_source(j, 0);
}
functionConstant::functionConstant(std::vector<double> source) : function(source.size())
{
- _source = fullMatrix<double>(source.size(),1);
- for (size_t i=0; i<source.size(); i++)
- _source(i,0) = source[i];
+ _source = fullMatrix<double>(source.size(), 1);
+ for (size_t i = 0; i < source.size(); i++)
+ _source(i, 0) = source[i];
}
functionConstant::functionConstant(double source) : function(1)
{
- _source.resize(1,1);
- _source(0,0) = source;
+ _source.resize(1, 1);
+ _source(0, 0) = source;
}
// functionSum
@@ -394,14 +394,15 @@ class functionSum : public function {
fullMatrix<double> _f0, _f1;
void call(dataCacheMap *m, fullMatrix<double> &val)
{
- for (int i=0; i<val.size1(); i++)
- for (int j=0; j<val.size2(); j++)
- val(i,j)= _f0(i,j) + _f1(i,j);
+ for (int i = 0; i < val.size1(); i++)
+ for (int j = 0; j < val.size2(); j++)
+ val(i, j)= _f0(i, j) + _f1(i, j);
}
functionSum(const function *f0, const function *f1) : function(f0->getNbCol())
{
if (f0->getNbCol() != f1->getNbCol()) {
- Msg::Error("trying to sum 2 functions of different sizes: %d %d\n",f0->getNbCol(),f1->getNbCol());
+ Msg::Error("trying to sum 2 functions of different sizes: %d %d\n",
+ f0->getNbCol(), f1->getNbCol());
throw;
}
setArgument (_f0, f0);
@@ -421,14 +422,15 @@ class functionProd : public function {
fullMatrix<double> _f0, _f1;
void call(dataCacheMap *m, fullMatrix<double> &val)
{
- for (int i=0; i<val.size1(); i++)
- for (int j=0; j<val.size2(); j++)
- val(i,j)= _f0(i,j)*_f1(i,j);
+ for (int i = 0; i < val.size1(); i++)
+ for (int j = 0; j < val.size2(); j++)
+ val(i, j)= _f0(i, j) * _f1(i, j);
}
functionProd(const function *f0, const function *f1) : function(f0->getNbCol())
{
if (f0->getNbCol() != f1->getNbCol()) {
- Msg::Error("trying to compute product of 2 functions of different sizes: %d %d\n",f0->getNbCol(),f1->getNbCol());
+ Msg::Error("trying to compute product of 2 functions of different sizes: %d %d\n",
+ f0->getNbCol(), f1->getNbCol());
throw;
}
setArgument (_f0, f0);
@@ -472,15 +474,15 @@ class functionCatComp : public function {
std::vector<fullMatrix<double> > _fMatrix;
void call(dataCacheMap *m, fullMatrix<double> &val)
{
- for (int i=0; i<val.size1(); i++)
- for (int comp=0; comp < _nComp; comp++)
- val(i,comp)= _fMatrix[comp](i,0);
+ for (int i = 0; i<val.size1(); i++)
+ for (int comp = 0; comp < _nComp; comp++)
+ val(i,comp) = _fMatrix[comp](i, 0);
}
functionCatComp(std::vector<const function *> fArray) : function(fArray.size())
{
- _nComp =fArray.size();
+ _nComp = fArray.size();
_fMatrix.resize(_nComp);
- for (int i=0; i<_nComp; i++)
+ for (int i = 0; i < _nComp; i++)
setArgument (_fMatrix[i], fArray[i]);
}
};
@@ -498,9 +500,9 @@ class functionScale : public function {
double _s;
void call(dataCacheMap *m, fullMatrix<double> &val)
{
- for(int i=0; i<val.size1(); i++)
- for(int j=0; j<val.size2(); j++)
- val(i,j)= _f0(i,j)*_s;
+ for(int i = 0; i < val.size1(); i++)
+ for(int j = 0; j < val.size2(); j++)
+ val(i, j) = _f0(i, j)*_s;
}
functionScale(const function *f0, const double s) : function(f0->getNbCol())
{
@@ -520,7 +522,7 @@ class functionCoordinates : public function {
fullMatrix<double> uvw;
void call (dataCacheMap *m, fullMatrix<double> &xyz)
{
- for (int i=0; i<uvw.size1(); i++) {
+ for (int i = 0; i < uvw.size1(); i++) {
SPoint3 p;
m->getElement()->pnt(uvw(i, 0), uvw(i, 1), uvw(i, 2), p);
xyz(i, 0) = p.x();
@@ -556,33 +558,33 @@ class functionStructuredGridFile : public function {
fullMatrix<double> coord;
public:
int n[3];
- double d[3],o[3];
- double get(int i,int j, int k)
+ double d[3], o[3];
+ double get(int i, int j, int k)
{
- return v[(i*n[1]+j)*n[2]+k];
+ return v[(i * n[1] + j) * n[2] + k];
}
double *v;
void call(dataCacheMap *m, fullMatrix<double> &val)
{
- for (int pt=0; pt<val.size1(); pt++) {
+ for (int pt = 0; pt < val.size1(); pt++) {
double xi[3];
int id[3];
- for(int i=0;i<3;i++){
- id[i]=(int)((coord(pt,i)-o[i])/d[i]);
- int a=id[i];
- id[i]=std::max(0,std::min(n[i]-1,id[i]));
- xi[i]=(coord(pt,i)-o[i]-id[i]*d[i])/d[i];
- xi[i]=std::min(1.,std::max(0.,xi[i]));
+ for(int i = 0; i < 3; i++){
+ id[i] = (int)((coord(pt, i) - o[i]) / d[i]);
+ int a = id[i];
+ id[i] = std::max(0, std::min(n[i] - 1, id[i]));
+ xi[i] = (coord(pt,i) - o[i] - id[i] * d[i]) / d[i];
+ xi[i] = std::min(1., std::max(0., xi[i]));
}
- val(pt,0) =
- get(id[0] ,id[1] ,id[2] )*(1-xi[0])*(1-xi[1])*(1-xi[2])
- +get(id[0] ,id[1] ,id[2]+1 )*(1-xi[0])*(1-xi[1])*( xi[2])
- +get(id[0] ,id[1]+1 ,id[2] )*(1-xi[0])*( xi[1])*(1-xi[2])
- +get(id[0] ,id[1]+1 ,id[2]+1 )*(1-xi[0])*( xi[1])*( xi[2])
- +get(id[0]+1 ,id[1] ,id[2] )*( xi[0])*(1-xi[1])*(1-xi[2])
- +get(id[0]+1 ,id[1] ,id[2]+1 )*( xi[0])*(1-xi[1])*( xi[2])
- +get(id[0]+1 ,id[1]+1 ,id[2] )*( xi[0])*( xi[1])*(1-xi[2])
- +get(id[0]+1 ,id[1]+1 ,id[2]+1 )*( xi[0])*( xi[1])*( xi[2]);
+ val(pt, 0) =
+ get(id[0] , id[1] , id[2] ) * (1-xi[0]) * (1-xi[1]) * (1-xi[2]) +
+ get(id[0] , id[1] , id[2]+1 ) * (1-xi[0]) * (1-xi[1]) * ( xi[2]) +
+ get(id[0] , id[1]+1, id[2] ) * (1-xi[0]) * ( xi[1]) * (1-xi[2]) +
+ get(id[0] , id[1]+1, id[2]+1 ) * (1-xi[0]) * ( xi[1]) * ( xi[2]) +
+ get(id[0]+1, id[1] , id[2] ) * ( xi[0]) * (1-xi[1]) * (1-xi[2]) +
+ get(id[0]+1, id[1] , id[2]+1 ) * ( xi[0]) * (1-xi[1]) * ( xi[2]) +
+ get(id[0]+1, id[1]+1, id[2] ) * ( xi[0]) * ( xi[1]) * (1-xi[2]) +
+ get(id[0]+1, id[1]+1, id[2]+1 ) * ( xi[0]) * ( xi[1]) * ( xi[2]);
}
}
functionStructuredGridFile(const std::string filename, const function *coordFunction)
@@ -592,12 +594,12 @@ class functionStructuredGridFile : public function {
std::ifstream input(filename.c_str());
if(!input)
Msg::Error("cannot open file : %s",filename.c_str());
- if(filename.substr(filename.size()-4,4)!=".bin") {
- input>>o[0]>>o[1]>>o[2]>>d[0]>>d[1]>>d[2]>>n[0]>>n[1]>>n[2];
- int nt = n[0]*n[1]*n[2];
+ if(filename.substr(filename.size()-4, 4) != ".bin") {
+ input >> o[0] >> o[1] >> o[2] >> d[0] >> d[1] >> d[2] >> n[0] >> n[1] >> n[2];
+ int nt = n[0] * n[1] * n[2];
v = new double [nt];
- for (int i=0; i<nt; i++)
- input>>v[i];
+ for (int i = 0; i < nt; i++)
+ input >> v[i];
} else {
input.read((char *)o, 3 * sizeof(double));
input.read((char *)d, 3 * sizeof(double));
@@ -624,7 +626,7 @@ class functionLua : public function {
void call (dataCacheMap *m, fullMatrix<double> &res)
{
lua_getfield(_L, LUA_GLOBALSINDEX, _luaFunctionName.c_str());
- for (int i=0;i< arguments.size(); i++)
+ for (int i = 0; i < arguments.size(); i++)
luaStack<const fullMatrix<double>*>::push(_L, &args[i]);
luaStack<const fullMatrix<double>*>::push(_L, &res);
lua_call(_L, arguments.size()+1, 0);
@@ -642,7 +644,7 @@ class functionLua : public function {
// functionC
void functionC::buildLibraryFromFile(const std::string cfilename, const std::string libfilename) {
- FILE *tmpMake = fopen("_tmpMake","w");
+ FILE *tmpMake = fopen("_tmpMake", "w");
fprintf(tmpMake, "include $(DG_BUILD_DIR)/CMakeFiles/dg.dir/flags.make\n"
"%s : %s\n"
"\tg++ -fPIC -shared -o $@ $(CXX_FLAGS) $(CXX_DEFINES) $<\n",
@@ -655,10 +657,10 @@ void functionC::buildLibraryFromFile(const std::string cfilename, const std::str
void functionC::buildLibrary(std::string code, std::string filename)
{
- FILE *tmpSrc = fopen("_tmpSrc.cpp","w");
+ FILE *tmpSrc = fopen("_tmpSrc.cpp", "w");
fprintf(tmpSrc, "%s\n",code.c_str());
fclose(tmpSrc);
- buildLibraryFromFile("_tmpSrc.cpp", filename);
+ buildLibraryFromFile("_tmpSrc.cpp", filename);
UnlinkFile("_tmpSrc.cpp");
}
void functionC::call (dataCacheMap *m, fullMatrix<double> &val)
@@ -710,11 +712,11 @@ functionC::functionC (std::string file, std::string symbol, int nbCol,
{
#if defined(HAVE_DLOPEN)
args.resize(dependencies.size());
- for(int i=0; i < dependencies.size(); i++) {
+ for(int i = 0; i < dependencies.size(); i++) {
setArgument(args[i], dependencies[i]);
}
void *dlHandler;
- dlHandler = dlopen(file.c_str(),RTLD_NOW);
+ dlHandler = dlopen(file.c_str(), RTLD_NOW);
callback = (void(*)(void))dlsym(dlHandler, symbol.c_str());
if(!callback)
Msg::Error("Cannot get the callback to the compiled C function");
@@ -738,78 +740,78 @@ void function::registerBindings(binding *b)
cb = b->addClass<functionConstant>("functionConstant");
cb->setDescription("A constant (scalar or vector) function");
- mb = cb->setConstructor<functionConstant,std::vector <double> >();
- mb->setArgNames("v",NULL);
+ mb = cb->setConstructor<functionConstant, std::vector <double> >();
+ mb->setArgNames("v", NULL);
mb->setDescription("A new constant function wich values 'v' everywhere. v can be a row-vector.");
cb->setParentClass<function>();
cb = b->addClass<functionSum>("functionSum");
cb->setDescription("A sum of two functions 'a + b'. The arguments a, b must have same dimension.");
- mb = cb->setConstructor<functionSum,const function*, const function*>();
- mb->setArgNames("a","b",NULL);
+ mb = cb->setConstructor<functionSum, const function*, const function*>();
+ mb->setArgNames("a", "b", NULL);
mb->setDescription("Creates a new functionSum instance with given arguments");
cb->setParentClass<function>();
cb = b->addClass<functionProd>("functionProd");
cb->setDescription("A pointwise product of two functions 'a(i,j)*b(i,j)'. The arguments a, b must have same dimension.");
- mb = cb->setConstructor<functionProd,const function*, const function*>();
- mb->setArgNames("a","b",NULL);
+ mb = cb->setConstructor<functionProd, const function*, const function*>();
+ mb->setArgNames("a", "b", NULL);
mb->setDescription("Creates a new functionProd instance with given arguments");
cb->setParentClass<function>();
cb = b->addClass<functionExtractComp>("functionExtractComp");
cb->setDescription("Extracts a given component of the vector valued function.");
- mb = cb->setConstructor<functionExtractComp,const function*, int>();
- mb->setArgNames("function","component",NULL);
+ mb = cb->setConstructor<functionExtractComp, const function*, int>();
+ mb->setArgNames("function","component", NULL);
mb->setDescription("Creates a new functionExtractComp instance with given arguments");
cb->setParentClass<function>();
cb = b->addClass<functionCatComp>("functionCatComp");
cb->setDescription("Creates a vector valued function by concatenating the given scalar functions. Uses only the first component of each function.");
- mb = cb->setConstructor<functionCatComp,std::vector <const function*> >();
- mb->setArgNames("functionArray",NULL);
+ mb = cb->setConstructor<functionCatComp, std::vector <const function*> >();
+ mb->setArgNames("functionArray", NULL);
mb->setDescription("Creates a new functionCatComp instance with given arguments");
cb->setParentClass<function>();
cb = b->addClass<functionScale>("functionScale");
cb->setDescription("Scales a function by a given scalar.");
- mb = cb->setConstructor<functionScale,const function*, double>();
- mb->setArgNames("function","scalar",NULL);
+ mb = cb->setConstructor<functionScale, const function*, double>();
+ mb->setArgNames("function","scalar", NULL);
mb->setDescription("Creates a new functionScale instance with given arguments");
cb->setParentClass<function>();
cb = b->addClass<functionCoordinates>("functionCoordinates");
cb->setDescription("A function to access the coordinates (xyz). This is a "
"single-instance class, use the 'get' member to access the instance.");
- mb = cb->addMethod("get",&functionCoordinates::get);
+ mb = cb->addMethod("get", &functionCoordinates::get);
mb->setDescription("return the unique instance of this class");
cb->setParentClass<function>();
cb = b->addClass<functionSolution>("functionSolution");
cb->setDescription("A function to access the solution. This is a single-instance "
"class, use the 'get' member to access the instance.");
- mb = cb->addMethod("get",&functionSolution::get);
+ mb = cb->addMethod("get", &functionSolution::get);
mb->setDescription("return the unique instance of this class");
cb->setParentClass<function>();
cb = b->addClass<functionNormals>("functionNormals");
cb->setDescription("A function to access the face normals. This is a single-instance "
"class, use the 'get' member to access the instance.");
- mb = cb->addMethod("get",&functionNormals::get);
+ mb = cb->addMethod("get", &functionNormals::get);
mb->setDescription("return the unique instance of this class");
cb->setParentClass<function>();
cb = b->addClass<functionSolutionGradient>("functionSolutionGradient");
cb->setDescription("A function to access the gradient of the solution. This is "
"a single-instance class, use the 'get' member to access the instance.");
- mb = cb->addMethod("get",&functionCoordinates::get);
+ mb = cb->addMethod("get", &functionCoordinates::get);
mb->setDescription("return the unique instance of this class");
cb->setParentClass<function>();
cb = b->addClass<functionStructuredGridFile>("functionStructuredGridFile");
cb->setParentClass<function>();
cb->setDescription("A function to interpolate through data given on a structured grid");
- mb = cb->setConstructor<functionStructuredGridFile,std::string, const function*>();
+ mb = cb->setConstructor<functionStructuredGridFile, std::string, const function*>();
mb->setArgNames("fileName","coordinateFunction",NULL);
mb->setDescription("Tri-linearly interpolate through data in file 'fileName' at "
"coordinate given by 'coordinateFunction'.\nThe file format "
@@ -818,11 +820,11 @@ void function::registerBindings(binding *b)
#if defined(HAVE_DLOPEN)
cb = b->addClass<functionC>("functionC");
cb->setDescription("A function that compile a C code");
- mb = cb->setConstructor<functionC,std::string, std::string,int,std::vector<const function*> >();
- mb->setArgNames("file", "symbol", "nbCol", "arguments",NULL);
+ mb = cb->setConstructor<functionC, std::string, std::string, int, std::vector<const function*> >();
+ mb->setArgNames("file", "symbol", "nbCol", "arguments", NULL);
mb->setDescription(" ");
mb = cb->addMethod("buildLibrary", &functionC::buildLibrary);
- mb->setArgNames("code", "libraryFileName",NULL);
+ mb->setArgNames("code", "libraryFileName", NULL);
mb->setDescription("build a dynamic library from the given code");
cb->setParentClass<function>();
#endif
@@ -830,7 +832,7 @@ void function::registerBindings(binding *b)
#if defined (HAVE_LUA)
cb= b->addClass<functionLua>("functionLua");
cb->setDescription("A function (see the 'function' documentation entry) defined in LUA.");
- mb = cb->setConstructor<functionLua,int,std::string,std::vector<const function*>,lua_State*>();
+ mb = cb->setConstructor<functionLua, int, std::string, std::vector<const function*>, lua_State*>();
mb->setArgNames("d", "f", "dep", NULL);
mb->setDescription("A new functionLua which evaluates a vector of dimension 'd' "
"using the lua function 'f'. This function can take other functions "
diff --git a/Solver/laplaceTerm.h b/Solver/laplaceTerm.h
index 94f44cf59b689c71b7af7acdfa42ae5ef6f49d5e..c1124a851ef345f89e668fb9e398a5ea7b17a507 100644
--- a/Solver/laplaceTerm.h
+++ b/Solver/laplaceTerm.h
@@ -23,7 +23,7 @@ class laplaceTerm : public helmholtzTerm<double> {
int nbSF = e->getNumShapeFunctions();
fullMatrix<double> *mat;
- mat = new fullMatrix<double>(nbSF,nbSF);
+ mat = new fullMatrix<double>(nbSF, nbSF);
elementMatrix(se, *mat);
fullVector<double> val(nbSF);
@@ -38,7 +38,7 @@ class laplaceTerm : public helmholtzTerm<double> {
m.scale(0.);
for (int i = 0; i < nbSF; i++)
for (int j = 0; j < nbSF; j++)
- m(i) += -(*mat)(i,j) * val(j);
+ m(i) += -(*mat)(i, j) * val(j);
}
};
diff --git a/Solver/orthogonalTerm.h b/Solver/orthogonalTerm.h
index 908f6c624dc48206a2f3e70936a04d3a47999ec0..00c4cec5305acc6df375cd2748e88c59faa9dce0 100644
--- a/Solver/orthogonalTerm.h
+++ b/Solver/orthogonalTerm.h
@@ -14,9 +14,9 @@ class orthogonalTerm : public helmholtzTerm<double> {
PView *_pview;
dofManager<double> *_dofView;
bool withDof;
- std::map<MVertex*,double > *_distance_map;
+ std::map<MVertex*, double > *_distance_map;
public:
- orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k, std::map<MVertex*,double > *distance_map)
+ orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k, std::map<MVertex*, double > *distance_map)
: helmholtzTerm<double>(gm, iField, iField, k, 0), _distance_map(distance_map), withDof(false) {}
orthogonalTerm(GModel *gm, int iField, simpleFunction<double> *k, PView *pview)
: helmholtzTerm<double>(gm, iField, iField, k, 0), _pview(pview), withDof(false) {}
@@ -92,7 +92,7 @@ class orthogonalTerm : public helmholtzTerm<double> {
m.scale(0.);
for(int i = 0; i < nbSF; i++)
for(int j = 0; j < nbSF; j++)
- m(i) += -mat(i,j) * val(j);
+ m(i) += -mat(i, j) * val(j);
}
};
diff --git a/Solver/quadratureRules.h b/Solver/quadratureRules.h
index e41f8510068f7729a36ecaad1ef270361053fb0a..880bfbb37fd69ae6a53450488175ec15bbc93a53 100644
--- a/Solver/quadratureRules.h
+++ b/Solver/quadratureRules.h
@@ -21,44 +21,44 @@ class QuadratureBase
{
public :
virtual ~QuadratureBase(){}
- virtual int getIntPoints(MElement *e,IntPt **GP) =0;
+ virtual int getIntPoints(MElement *e, IntPt **GP) =0;
};
class GaussQuadrature : public QuadratureBase
{
public :
- enum IntegCases {Other,Val,Grad,ValVal,GradGrad};
+ enum IntegCases {Other, Val, Grad, ValVal, GradGrad};
private :
int order;
IntegCases info;
public :
- GaussQuadrature(int order_=0):order(order_),info(Other) {}
- GaussQuadrature(IntegCases info_):order(0),info(info_) {}
+ GaussQuadrature(int order_ = 0) : order(order_), info(Other) {}
+ GaussQuadrature(IntegCases info_) : order(0), info(info_) {}
virtual ~GaussQuadrature(){}
- int getIntPoints(MElement *e,IntPt **GP)
+ int getIntPoints(MElement *e, IntPt **GP)
{
int integrationOrder;
int npts;
- int geoorder=e->getPolynomialOrder();
+ int geoorder = e->getPolynomialOrder();
switch(info)
{
case Other :
integrationOrder = order;
break;
case Val :
- integrationOrder=geoorder+1;
+ integrationOrder = geoorder + 1;
break;
case Grad :
- integrationOrder=geoorder;
+ integrationOrder = geoorder;
break;
case ValVal :
- integrationOrder=2*geoorder;
+ integrationOrder = 2 * geoorder;
break;
case GradGrad :
- integrationOrder=3*(geoorder-1)+1;
+ integrationOrder = 3 * (geoorder - 1) + 1;
break;
- default : integrationOrder=1;
+ default : integrationOrder = 1;
}
e->getIntegrationPoints(integrationOrder, &npts, GP);
return npts;
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index 98a466872d295229079df1d1e6fb00460ec6099d..3b3f7bbfac21dc4eebd1968addfd32fd441df2a2 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -22,53 +22,54 @@
-template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler) // symmetric
+template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term, FunctionSpaceBase &space,
+ Iterator itbegin, Iterator itend,
+ QuadratureBase &integrator, Assembler &assembler)
+ // symmetric
{
fullMatrix<typename Assembler::dataMat> localMatrix;
std::vector<Dof> R;
- for (Iterator it = itbegin;it!=itend; ++it)
- {
+ for (Iterator it = itbegin; it != itend; ++it){
MElement *e = *it;
R.clear();
IntPt *GP;
- int npts=integrator.getIntPoints(e,&GP);
- term.get(e,npts,GP,localMatrix);
- space.getKeys(e,R);
+ int npts = integrator.getIntPoints(e, &GP);
+ term.get(e, npts, GP, localMatrix); //localMatrix.print();
+ space.getKeys(e, R);
assembler.assemble(R, localMatrix);
}
}
-template<class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler) // symmetric
+template<class Assembler> void Assemble(BilinearTermBase &term, FunctionSpaceBase &space, MElement *e,
+ QuadratureBase &integrator, Assembler &assembler) // symmetric
{
fullMatrix<typename Assembler::dataMat> localMatrix;
std::vector<Dof> R;
IntPt *GP;
- int npts=integrator.getIntPoints(e,&GP);
- term.get(e,npts,GP,localMatrix);
- space.getKeys(e,R);
+ int npts = integrator.getIntPoints(e, &GP);
+ term.get(e, npts, GP, localMatrix);
+ space.getKeys(e, R);
assembler.assemble(R, localMatrix);
}
-template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,
- FunctionSpaceBase &shapeFcts,
- FunctionSpaceBase &testFcts,
- Iterator itbegin,
- Iterator itend,
- QuadratureBase &integrator,
- Assembler &assembler) // non symmetric
+template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term,
+ FunctionSpaceBase &shapeFcts,
+ FunctionSpaceBase &testFcts,
+ Iterator itbegin, Iterator itend,
+ QuadratureBase &integrator,
+ Assembler &assembler) // non symmetric
{
fullMatrix<typename Assembler::dataMat> localMatrix;
std::vector<Dof> R;
std::vector<Dof> C;
- for (Iterator it = itbegin; it != itend; ++it)
- {
+ for (Iterator it = itbegin; it != itend; ++it){
MElement *e = *it;
R.clear();
C.clear();
IntPt *GP;
int npts = integrator.getIntPoints(e, &GP);
- term.get(e, npts, GP, localMatrix);
- printf("local matrix size = %d %d\n",localMatrix.size1(),localMatrix.size2());
+ term.get(e, npts, GP, localMatrix); //localMatrix.print();
+ printf("local matrix size = %d %d\n", localMatrix.size1(), localMatrix.size2());
shapeFcts.getKeys(e, R);
testFcts.getKeys(e, C);
// std::cout << "assembling normal test function ; lagrange trial function : " << std::endl;
@@ -96,69 +97,73 @@ template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,
-template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler)
+template<class Iterator, class Assembler> void Assemble(LinearTermBase &term, FunctionSpaceBase &space,
+ Iterator itbegin, Iterator itend,
+ QuadratureBase &integrator, Assembler &assembler)
{
fullVector<typename Assembler::dataMat> localVector;
std::vector<Dof> R;
- for (Iterator it = itbegin;it!=itend; ++it)
- {
+ for (Iterator it = itbegin; it != itend; ++it){
MElement *e = *it;
R.clear();
IntPt *GP;
- int npts=integrator.getIntPoints(e,&GP);
- term.get(e,npts,GP,localVector);
- space.getKeys(e,R);
+ int npts = integrator.getIntPoints(e, &GP);
+ term.get(e, npts, GP, localVector); //localVector.print();
+ space.getKeys(e, R);
assembler.assemble(R, localVector);
}
}
-template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler)
+template<class Assembler> void Assemble(LinearTermBase &term, FunctionSpaceBase &space, MElement *e,
+ QuadratureBase &integrator, Assembler &assembler)
{
fullVector<typename Assembler::dataMat> localVector;
std::vector<Dof> R;
IntPt *GP;
- int npts=integrator.getIntPoints(e,&GP);
- term.get(e,npts,GP,localVector);
- space.getKeys(e,R);
+ int npts = integrator.getIntPoints(e, &GP);
+ term.get(e, npts, GP, localVector);
+ space.getKeys(e, R);
assembler.assemble(R, localVector);
}
-template<class Iterator,class dataMat> void Assemble(ScalarTermBase &term,Iterator itbegin,Iterator itend,QuadratureBase &integrator,dataMat & val)
+template<class Iterator, class dataMat> void Assemble(ScalarTermBase &term,
+ Iterator itbegin, Iterator itend,
+ QuadratureBase &integrator, dataMat & val)
{
dataMat localval;
- for (Iterator it = itbegin;it!=itend; ++it)
- {
+ for (Iterator it = itbegin; it != itend; ++it){
MElement *e = *it;
IntPt *GP;
- int npts=integrator.getIntPoints(e,&GP);
- term.get(e,npts,GP,localval);
- val+=localval;
+ int npts = integrator.getIntPoints(e, &GP);
+ term.get(e, npts, GP, localval);
+ val += localval;
}
}
-template<class Iterator,class dataMat> void Assemble(ScalarTermBase &term,MElement *e,QuadratureBase &integrator,dataMat & val)
+template<class Iterator, class dataMat> void Assemble(ScalarTermBase &term, MElement *e,
+ QuadratureBase &integrator, dataMat & val)
{
dataMat localval;
IntPt *GP;
- int npts=integrator.getIntPoints(e,&GP);
- term.get(e,npts,GP,localval);
- val+=localval;
+ int npts = integrator.getIntPoints(e, &GP);
+ term.get(e, npts, GP, localval);
+ val += localval;
}
-template<class Assembler> void FixDofs(Assembler &assembler,std::vector<Dof> &dofs,std::vector<typename Assembler::dataVec> &vals)
+template<class Assembler> void FixDofs(Assembler &assembler, std::vector<Dof> &dofs,
+ std::vector<typename Assembler::dataVec> &vals)
{
- int nbff=dofs.size();
- for (int i=0;i<nbff;++i)
- {
- assembler.fixDof(dofs[i],vals[i]);
+ int nbff = dofs.size();
+ for (int i = 0; i < nbff; ++i){
+ assembler.fixDof(dofs[i], vals[i]);
}
}
class FilterDof
{
public:
- virtual bool operator()(Dof key)=0;
+ virtual bool operator()(Dof key) = 0;
};
class FilterDofTrivial : public FilterDof
@@ -167,41 +172,40 @@ class FilterDofTrivial : public FilterDof
virtual bool operator()(Dof key) {return true;}
};
-class FilterDofComponent :public FilterDof
+class FilterDofComponent : public FilterDof
{
int comp;
public :
- FilterDofComponent(int comp_):comp(comp_) {}
+ FilterDofComponent(int comp_) : comp(comp_) {}
virtual bool operator()(Dof key)
{
- int type=key.getType();
- int icomp,iphys;
+ int type = key.getType();
+ int icomp, iphys;
Dof::getTwoIntsFromType(type, icomp, iphys);
- if (icomp==comp) return true;
+ if (icomp == comp) return true;
return false;
}
};
-template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MElement *e,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter)
+template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space, MElement *e, Assembler &assembler,
+ simpleFunction<typename Assembler::dataVec> &fct,
+ FilterDof &filter)
{
std::vector<MVertex*> tabV;
- int nv=e->getNumVertices();
+ int nv = e->getNumVertices();
std::vector<Dof> R;
- space.getKeys(e,R);
+ space.getKeys(e, R);
tabV.reserve(nv);
- for (int i=0;i<nv;++i) tabV.push_back(e->getVertex(i));
+ for (int i = 0; i < nv; ++i)
+ tabV.push_back(e->getVertex(i));
- for (std::vector<Dof>::iterator itd=R.begin();itd!=R.end();++itd)
- {
- Dof key=*itd;
- if (filter(key))
- {
- for (int i=0;i<nv;++i)
- {
- if (tabV[i]->getNum()==key.getEntity())
- {
- assembler.fixDof(key, fct(tabV[i]->x(),tabV[i]->y(),tabV[i]->z()));
+ for (std::vector<Dof>::iterator itd = R.begin(); itd != R.end(); ++itd){
+ Dof key = *itd;
+ if (filter(key)){
+ for (int i = 0;i < nv; ++i){
+ if (tabV[i]->getNum() == key.getEntity()){
+ assembler.fixDof(key, fct(tabV[i]->x(), tabV[i]->y(), tabV[i]->z()));
break;
}
}
@@ -209,53 +213,51 @@ template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MElement *e
}
}
-template<class Iterator,class Assembler> void FixNodalDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter)
+template<class Iterator, class Assembler> void FixNodalDofs(FunctionSpaceBase &space, Iterator itbegin,
+ Iterator itend, Assembler &assembler,
+ simpleFunction<typename Assembler::dataVec> &fct,
+ FilterDof &filter)
{
- for (Iterator it=itbegin;it!=itend;++it)
- {
- FixNodalDofs(space,*it,assembler,fct,filter);
- }
+ for (Iterator it = itbegin; it != itend; ++it)
+ FixNodalDofs(space, *it, assembler, fct, filter);
}
-template<class Iterator,class Assembler> void FixVoidNodalDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
+template<class Iterator, class Assembler> void FixVoidNodalDofs(FunctionSpaceBase &space, Iterator itbegin,
+ Iterator itend, Assembler &assembler)
{
FilterDofTrivial filter;
simpleFunction<typename Assembler::dataVec> fct(0.);
FixNodalDofs(space, itbegin, itend, assembler, fct, filter);
}
-template<class Iterator,class Assembler> void NumberDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
+template<class Iterator, class Assembler> void NumberDofs(FunctionSpaceBase &space, Iterator itbegin,
+ Iterator itend, Assembler &assembler)
{
- for (Iterator it=itbegin;it!=itend;++it)
- {
- MElement *e=*it;
+ for (Iterator it = itbegin; it != itend; ++it){
+ MElement *e = *it;
std::vector<Dof> R;
- space.getKeys(e,R);
- int nbdofs=R.size();
- for (int i=0;i<nbdofs;++i) assembler.numberDof(R[i]);
-
+ space.getKeys(e, R);
+ int nbdofs = R.size();
+ for (int i = 0; i < nbdofs; ++i) assembler.numberDof(R[i]);
}
}
//// Mean HangingNodes
-//template <class Assembler> void FillHangingNodes(FunctionSpaceBase &space,std::map<int,std::vector <int> > &HangingNodes, Assembler &assembler, int &field, int &dim)
+//template <class Assembler> void FillHangingNodes(FunctionSpaceBase &space, std::map<int,std::vector <int> > &HangingNodes, Assembler &assembler, int &field, int &dim)
//{
-// std::map<int,std::vector <int> >::iterator ith;
+// std::map<int, std::vector <int> >::iterator ith;
// ith = HangingNodes.begin();
// int compt = 1;
-// while (ith!=HangingNodes.end())
-// {
+// while (ith != HangingNodes.end()){
// float fac;
-// fac = 1.0/(ith->second).size();
-// for (int j=0;j<dim;j++)
-// {
+// fac = 1.0 / (ith->second).size();
+// for (int j = 0; j < dim; j++){
// DofAffineConstraint<double> constraint;
// int type = Dof::createTypeWithTwoInts(j, field);
-// Dof hgnd(ith->first,type);
-// for (int i=0;i<(ith->second).size();i++)
-// {
-// Dof key((ith->second)[i],type);
-// std::pair<Dof, double > linDof(key,fac);
+// Dof hgnd(ith->first, type);
+// for (int i = 0; i < (ith->second).size(); i++){
+// Dof key((ith->second)[i], type);
+// std::pair<Dof, double > linDof(key, fac);
// constraint.linear.push_back(linDof);
// }
// constraint.shift = 0;
diff --git a/Solver/solverField.h b/Solver/solverField.h
index bdd6f57b0d3c58e4e189e8e21e8e09f48ae2470a..4e31277784b700b4eaca1e8338b8fd9231e2f321 100644
--- a/Solver/solverField.h
+++ b/Solver/solverField.h
@@ -31,7 +31,7 @@ class SolverField : public FunctionSpace<T> // being able to use it instead of a
dofManager<double> *dm;
FunctionSpace<T> *fs;
public:
- SolverField(dofManager<double> *dm_, FunctionSpace<T> *fs_) : dm(dm_),fs(fs_) {}
+ SolverField(dofManager<double> *dm_, FunctionSpace<T> *fs_) : dm(dm_), fs(fs_) {}
virtual int getNumKeys(MVertex *ver) { return 1;}
virtual int getNumKeys(MElement *ele) { return 1;}
private:
@@ -44,60 +44,60 @@ class SolverField : public FunctionSpace<T> // being able to use it instead of a
std::vector<Dof> D;
std::vector<ValType> SFVals;
std::vector<double> DMVals;
- fs->getKeys(ele,D);
- dm->getDofValue(D,DMVals);
- fs->f(ele,u,v,w,SFVals);
- val=ValType();
- for (unsigned int i=0;i<D.size();++i)
- val+=SFVals[i]*DMVals[i];
+ fs->getKeys(ele, D);
+ dm->getDofValue(D, DMVals);
+ fs->f(ele, u, v, w, SFVals);
+ val = ValType();
+ for (unsigned int i = 0; i < D.size(); ++i)
+ val += SFVals[i] * DMVals[i];
}
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
ValType val;
- f(ele,u,v,w,val);
+ f(ele, u, v, w, val);
vals.push_back(val);
}
- virtual void gradf(MElement *ele, double u, double v, double w,GradType &grad)
+ virtual void gradf(MElement *ele, double u, double v, double w, GradType &grad)
{
std::vector<Dof> D;
std::vector<GradType> SFGrads;
std::vector<double> DMVals;
- fs->getKeys(ele,D);
- dm->getDofValue(D,DMVals);
- fs->gradf(ele,u,v,w,SFGrads);
- grad=GradType();
- for (unsigned int i=0;i<D.size();++i)
- grad+= SFGrads[i] * DMVals[i];
+ fs->getKeys(ele, D);
+ dm->getDofValue(D, DMVals);
+ fs->gradf(ele, u, v, w, SFGrads);
+ grad = GradType();
+ for (unsigned int i = 0; i < D.size(); ++i)
+ grad += SFGrads[i] * DMVals[i];
}
// A quoi sert cette fonction ?? (Evalue le hessien au noeuds
-/* virtual void hessf(MElement *ele, double u, double v, double w,HessType &hess)
+/* virtual void hessf(MElement *ele, double u, double v, double w, HessType &hess)
{
// Pas besoin des dof etc
// std::vector<Dof> D;
std::vector<HessType> SFHess;
// std::vector<double> DMVals;
-// fs->getKeys(ele,D);
-// dm->getDofValue(D,DMVals);
- fs->hessf(ele,u,v,w,SFHess);
+// fs->getKeys(ele, D);
+// dm->getDofValue(D, DMVals);
+ fs->hessf(ele, u, v, w, SFHess);
-// hess=HessType;
-// for (int i=0;i<D.size();++i)
-// hess+= SFHess[i] * DMVals[i];
+// hess = HessType;
+// for (int i = 0; i < D.size(); ++i)
+// hess += SFHess[i] * DMVals[i];
}*/
- virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+ virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
{
GradType grad;
- gradf(ele,u,v,w,grad);
+ gradf(ele, u, v, w, grad);
grads.push_back(grad);
}
- virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess)
+ virtual void hessfuvw(MElement *ele, double u, double v, double w, std::vector<HessType> &hess)
{
//HessType hes;
- fs->hessfuvw(ele,u,v,w,hess);
+ fs->hessfuvw(ele, u, v, w, hess);
}
};
@@ -108,7 +108,7 @@ class Formulation
std::vector<FunctionSpace<double>* > scalarfs;
std::vector<FunctionSpace<SVector3>* > vectorfs;
std::vector<groupOfElements* > groups;
- std::vector<std::pair<MElement*,std::vector<groupOfElements&> > > links;
+ std::vector<std::pair<MElement*, std::vector<groupOfElements&> > > links;
dofManager<double> *dm; //
};