diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 354b91f7765a4d30f04f9e0cfd5f9e65f144cfb7..6a62a8e1f8d65b4f7eac758862042f52dfe8eed2 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -30,6 +30,7 @@
#include "GmshMessage.h"
#include "linearSystem.h"
#include "Options.h"
+#include "elasticitySolver.h"
#if defined(HAVE_OPENGL)
#include "drawContext.h"
@@ -409,6 +410,7 @@ binding::binding()
#if defined(HAVE_SOLVER)
function::registerBindings(this);
linearSystemCSRGmm<double>::registerBindings(this);
+ elasticitySolverRegisterBindings(this);
#endif
}
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 9b82c5526aeac43bebf0d7bba72ef7e317e9f0ea..e0b4a2bbb6b8085eb0b601c0940115b46d3a260a 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -26,7 +26,6 @@
#include "GEntity.h"
#include "dofManager.h"
#include "laplaceTerm.h"
-#include "distanceTerm.h"
#include "crossConfTerm.h"
#include "convexCombinationTerm.h"
#include "diagBCTerm.h"
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index d6ba6b25d4b8b91af3c721db933d25bd699778c4..1a4c13a7d9556b5b62b6914be84e1818fad96495 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -694,6 +694,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
fprintf(fp, "%d\n", numElements);
int num = elementStartNum;
std::map<MElement*, int> parentsNum;
+
_elementIndexCache.clear();
// points
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index f885578e7d1e8101dc51a1d0d2679ea8cface95f..3d2cb90f24431058ab1960145c6993b461176904 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -203,8 +203,8 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
norme(b);
prodve(a, b, c);
norme(c);
- jac[1][0] = b[0]; jac[1][1] = b[1]; jac[1][2] = b[2];
- jac[2][0] = c[0]; jac[2][1] = c[1]; jac[2][2] = c[2];
+ jac[0][1] = b[0]; jac[1][1] = b[1]; jac[2][1] = b[2];
+ jac[0][2] = c[0]; jac[1][2] = c[1]; jac[2][2] = c[2];
break;
}
case 2:
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index 1cac248c5a4cf236139c9b1f042494cfe9a96caf..f3a78985f55a07fa7da08407c775419fe1648af2 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -24,6 +24,7 @@ class backgroundMesh : public simpleFunction<double>
std::map<MVertex*,double> _sizes;
std::map<MVertex*,MVertex*> _3Dto2D;
std::map<MVertex*,MVertex*> _2Dto3D;
+ std::map<MVertex*,double> _distance;
static backgroundMesh * _current;
backgroundMesh(GFace *);
~backgroundMesh();
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 9cd56ca934d26f855ab66c0bd13940783ba9d023..66a1339928bb1cb070fcf1421de10f026bb782c6 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -668,6 +668,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
const double MINE_ = 0.67, MAXE_ = 1.4;
while (1){
+
// we count the number of local mesh modifs.
int nb_split = 0;
int nb_smooth = 0;
@@ -697,18 +698,29 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
double minE = MINE_;
double t1 = Cpu();
splitEdgePass(gf, m, maxE, nb_split);
+
double t2 = Cpu();
swapEdgePass(gf, m, nb_swap);
swapEdgePass(gf, m, nb_swap);
swapEdgePass(gf, m, nb_swap);
+
+ // if (computeNodalSizeField){
+ // char name[256]; sprintf(name,"iter%d_SPLIT.pos",IT);
+ // outputScalarField(m.triangles, name, 0);
+ // }
double t3 = Cpu();
collapseEdgePass(gf, m, minE, MAXNP, nb_collaps);
+
double t4 = Cpu();
double t5 = Cpu();
smoothVertexPass(gf, m, nb_smooth, false);
double t6 = Cpu();
swapEdgePass ( gf, m, nb_swap);
double t7 = Cpu();
+ // if (computeNodalSizeField){
+ // char name[256]; sprintf(name,"iter%d_COLLAPSE.pos",IT);
+ // outputScalarField(m.triangles, name, 0);
+ // }
// clean up the mesh
t_spl += t2 - t1;
t_sw += t3 - t2;
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 2d79ea40227b83c24ff7fb419f21e127de9cfa1e..a6ee7a905095451d823d2e8f276fb85fbad21e13 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -656,6 +656,20 @@ static void insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
}
+
+static void insertManyPoints(GFace *gf,
+ std::list<SPoint2> &points,
+ std::vector<double> &Us,
+ std::vector<double> &Vs,
+ std::vector<double> &vSizes,
+ std::vector<double> &vSizesBGM,
+ std::vector<SMetric3> &vMetricsBGM,
+ std::set<MTri3*,compareTri3Ptr> &AllTris,
+ std::set<MTri3*,compareTri3Ptr> *ActiveTris = 0){
+
+}
+
+
void bowyerWatson(GFace *gf)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
@@ -677,6 +691,11 @@ void bowyerWatson(GFace *gf)
int ITER = 0;
while (1){
+ // if(ITER % 1== 0){
+ // char name[245];
+ // sprintf(name,"del2d%d-ITER%4d.pos",gf->tag(),ITER);
+ // _printTris (name, AllTris, Us,Vs,false);
+ // }
MTri3 *worst = *AllTris.begin();
if (worst->isDeleted()){
delete worst->tri();
@@ -688,7 +707,7 @@ void bowyerWatson(GFace *gf)
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
vSizes.size(), worst->getRadius());
double center[2],metric[3],r2;
- if (worst->getRadius() < 0.5 * sqrt(2.0)) break;
+ if (worst->getRadius() < /*1.333333/(sqrt(3.0))*/0.5 * sqrt(2.0)) break;
circUV(worst->tri(), Us, Vs, center, gf);
MTriangle *base = worst->tri();
double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
@@ -710,11 +729,6 @@ void bowyerWatson(GFace *gf)
insertAPoint(gf, AllTris.begin(), center, metric, Us, Vs, vSizes,
vSizesBGM, vMetricsBGM, AllTris);
}
- // if(ITER % 10== 0){
- // char name[245];
- // sprintf(name,"del2d%d-ITER%d.pos",gf->tag(),ITER);
- // _printTris (name, AllTris, Us,Vs,false);
- // }
}
transferDataStructure(gf, AllTris, Us, Vs);
}
@@ -802,6 +816,15 @@ void bowyerWatsonFrontal(GFace *gf)
// insert points
while (1){
+ /*
+ if(ITER % 1== 0){
+ char name[245];
+ sprintf(name,"delfr2d%d-ITER%4d.pos",gf->tag(),ITER);
+ _printTris (name, AllTris, Us,Vs,false);
+ sprintf(name,"delfr2dA%d-ITER%4d.pos",gf->tag(),ITER);
+ _printTris (name, ActiveTris, Us,Vs,false);
+ }
+ */
if (!ActiveTris.size())break;
MTri3 *worst = (*ActiveTris.begin());
ActiveTris.erase(ActiveTris.begin());
@@ -882,3 +905,12 @@ void bowyerWatsonFrontal(GFace *gf)
// _printTris (name, AllTris, Us, Vs,true);
transferDataStructure(gf, AllTris, Us, Vs);
}
+
+void addBoundaryLayers(GFace *gf) {
+ // first compute the distance function u on the existing mesh
+ // then compute the dual function v that has gradients orthogonal everywhere
+ // build a set of points in the boundary layer
+ // connect everybody with delaunay
+
+}
+
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index c258266eed4899632b4968e30cd0f49b39c37d6b..ad53fd9bd3933b31a26667a1b0b5fd1baa710870 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -9,6 +9,7 @@
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemGMM.h"
+#include "linearSystemFull.h"
#include "Numeric.h"
#include "GModel.h"
#include "functionSpace.h"
@@ -208,6 +209,8 @@ void elasticitySolver::readInputFile(const std::string &fn)
void elasticitySolver::solve()
{
+ linearSystemFull<double> *lsys = new linearSystemFull<double>;
+ /*
#if defined(HAVE_TAUCS)
linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#elif defined(HAVE_PETSC)
@@ -216,6 +219,7 @@ void elasticitySolver::solve()
linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
lsys->setNoisy(2);
#endif
+ */
if (pAssembler) delete pAssembler;
pAssembler = new dofManager<double>(lsys);
@@ -277,6 +281,15 @@ void elasticitySolver::solve()
Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
}
+ for (int i=0;i<pAssembler->sizeOfR();i++){
+ for (int j=0;j<pAssembler->sizeOfR();j++){
+ double d;lsys->getFromMatrix(i,j,d);
+ printf("%12.5E ",d);
+ }
+ double d;lsys->getFromRightHandSide(i,d);
+ printf(" | %12.5E\n",d);
+ }
+
printf("nDofs=%d\n",pAssembler->sizeOfR());
printf("nFixed=%d\n",pAssembler->sizeOfF());
printf("-- done assembling!\n");
@@ -493,3 +506,25 @@ PView* elasticitySolver::buildVonMisesView(const std::string &postFileName)
}
#endif
+
+#include "Bindings.h"
+void elasticitySolverRegisterBindings(binding *b)
+{
+ classBinding *cb;
+ cb = b->addClass<elasticitySolver> ("elasticitySolver");
+ cb->setDescription("A class that enables to solve elasticity problems");
+
+ methodBinding *cm;
+ cm = cb->addMethod("read", &elasticitySolver::read);
+ cm->setDescription ("reads an input file");
+ cm->setArgNames("fileName",NULL);
+
+ cm = cb->addMethod("solve", &elasticitySolver::solve);
+ cm->setDescription ("solves the problem");
+
+ cm = cb->setConstructor<elasticitySolver,int>();
+ cm->setDescription ("A new elasticitySolver. The parameter is the unknowns tag");
+ cm->setArgNames("tag",NULL);
+
+}
+
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index b70f3c25383b7b45574ffc45b16cf3dbad0476b7..f5866f1d4e6bcda0bc98200e6b3acfb0c3cc1d64 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -82,8 +82,9 @@ class elasticitySolver
if (pAssembler) delete pAssembler;
}
void readInputFile(const std::string &meshFileName);
+ void read(char *s) {readInputFile(s);}
virtual void setMesh(const std::string &meshFileName);
- virtual void solve();
+ void solve();
virtual PView *buildDisplacementView(const std::string &postFileName);
virtual PView *buildLagrangeMultiplierView(const std::string &posFileName);
virtual PView *buildElasticEnergyView(const std::string &postFileName);
@@ -94,5 +95,7 @@ class elasticitySolver
// (const std::string &errorFileName, const elasticityData &ref, double, int);
};
+class binding;
+void elasticitySolverRegisterBindings(binding *b);
#endif
diff --git a/Solver/linearSystemFull.h b/Solver/linearSystemFull.h
index eb8c17839b76b9a2259d2a5b72b947ef52b843c7..641f0b02a76975aa81dae7d2a0aa512fd39137b1 100644
--- a/Solver/linearSystemFull.h
+++ b/Solver/linearSystemFull.h
@@ -84,6 +84,7 @@ class linearSystemFull : public linearSystem<scalar> {
{
if (_b->size())
_a->luSolve(*_b, *_x);
+ // _x->print("X in solve");
return 1;
}
};
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index 1c76468a9f22925d7d23865489e709d8b5a309af..98a466872d295229079df1d1e6fb00460ec6099d 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -68,6 +68,7 @@ template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,
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());
shapeFcts.getKeys(e, R);
testFcts.getKeys(e, C);
// std::cout << "assembling normal test function ; lagrange trial function : " << std::endl;