diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index f2d2cefc793f61889587163e545f568e1534e2d5..df6946665f545c133d38cafdc9bedf6ead597be2 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -521,11 +521,9 @@ bool GFaceCompound::parametrize() const
bool withoutFolding = parametrize_conformal_spectral() ;
//bool withoutFolding = parametrize_conformal();
if ( withoutFolding == false ){
-
- double alpha = 0.0;
- Msg::Warning("$$$ Parametrization switched to combination map: A*conf+(1-A)*harm with A=%g", alpha);
- parametrize(ITERU,HARMONIC, alpha);
- parametrize(ITERV,HARMONIC, alpha);
+ Msg::Warning("$$$ Parametrization switched to harmonic map");
+ parametrize(ITERU,HARMONIC);
+ parametrize(ITERV,HARMONIC);
//buildOct(); exit(1);
}
}
@@ -1124,9 +1122,22 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom, double al
}
bool GFaceCompound::parametrize_conformal_spectral() const
+{
+
+#if !defined(HAVE_PETSC) && !defined(HAVE_SLEPC)
+{
+
+ Msg::Error("-----------------------------------------------------------------------------!");
+ Msg::Error("Gmsh should be compiled with petsc and slepc for using the conformal map !");
+ Msg::Error("Switch to harmonic map or see doc on the wiki for installing petsc and slepc !");
+ Msg::Error("https://geuz.org/trac/gmsh/wiki/STLRemeshing (username:gmsh,passwd:gmsh) !");
+ Msg::Error("-----------------------------------------------------------------------------!");
+ Msg::Exit(1);
+
+}
+#else
{
-#if defined(HAVE_PETSC)
std::vector<MVertex*> ordered;
std::vector<double> coords;
bool success = orderVertices(_U0, ordered, coords);
@@ -1236,9 +1247,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
}
else return false;
-#else
- return false;
-
+}
#endif
}
bool GFaceCompound::parametrize_conformal() const
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index cf6c0b4f5ca4f7877750f5b6687c3bb2731c9a9e..f5c3987e385b20ae15cf15b56e8cbd3f4e8a75e9 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -22,6 +22,7 @@
#if defined(HAVE_POST)
#include "PView.h"
#include "PViewData.h"
+#include "function.h"
#endif
static void printLinearSystem(linearSystemCSRTaucs<double> *lsys)
@@ -67,8 +68,9 @@ void elasticitySolver::setMesh(const std::string &meshFileName)
void elasticitySolver::solve()
{
- linearSystemFull<double> *lsys = new linearSystemFull<double>;
- /*
+
+ //linearSystemFull<double> *lsys = new linearSystemFull<double>;
+
#if defined(HAVE_TAUCS)
linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#elif defined(HAVE_PETSC)
@@ -77,7 +79,7 @@ void elasticitySolver::solve()
linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
lsys->setNoisy(2);
#endif
- */
+
assemble(lsys);
lsys->systemSolve();
printf("-- done solving!\n");
@@ -422,13 +424,13 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
printf("nDofs=%d\n",pAssembler->sizeOfR());
printf("nFixed=%d\n",pAssembler->sizeOfF());
- printf("-- done assembling!\n");
+
}
#if defined(HAVE_POST)
static void deformation(dofManager<double> *a, MElement *e,
- double u, double v, double w, int _tag, double *eps){
+ double u, double v, double w, int _tag, double *eps){
double valx[256];
double valy[256];
double valz[256];
@@ -454,8 +456,8 @@ static void deformation(dofManager<double> *a, MElement *e,
static double vonMises(dofManager<double> *a, MElement *e,
double u, double v, double w,
- double E, double nu, int _tag)
-{
+ double E, double nu, int _tag){
+
double valx[256];
double valy[256];
double valz[256];
@@ -485,9 +487,7 @@ static double vonMises(dofManager<double> *a, MElement *e,
double sxz = A * eps[4];
double syz = A * eps[5];
- double s[9] = {sxx, sxy, sxz,
- sxy, syy, syz,
- sxz, syz, szz};
+ double s[9] = {sxx, sxy, sxz, sxy, syy, syz,sxz, syz, szz};
return ComputeVonMises(s);
}
@@ -540,6 +540,71 @@ PView* elasticitySolver::buildDisplacementView (const std::string postFileName)
return pv;
}
+PView* elasticitySolver::buildStressesView (const std::string postFileName)
+{
+ std::cout << "build stresses view"<< std::endl;
+ std::map<int, std::vector<double> > data;
+ GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
+ for (unsigned int i = 0; i < elasticFields.size(); ++i)
+ {
+ double E = elasticFields[i]._E;
+ double nu = elasticFields[i]._nu;
+ SolverField<SVector3> Field(pAssembler, LagSpace);
+ IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+ BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+ for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
+ {
+ MElement *e=*it;
+ int nbVertex = e->getNumVertices();
+ std::vector<SVector3> val(nbVertex);
+
+ double valx[256];
+ double valy[256];
+ double valz[256];
+ for (int k = 0; k < nbVertex; k++){
+ MVertex *v = e->getVertex(k);
+ MPoint p(v);
+ Field.f(&p,0,0,0,val[k]);
+ valx[k] =val[k](0);
+ valy[k] =val[k](1);
+ valz[k] =val[k](2);
+ }
+
+ double gradux[3];
+ double graduy[3];
+ double graduz[3];
+ double u=0.33, v=0.33, w=0.0;
+ 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]),
+ 0.5 * (graduy[2] + graduz[1])};
+
+
+ double A = E / (1. + nu);
+ double B = A * (nu / (1. - 2 * nu));
+ double trace = eps[0] + eps[1] + eps[2] ;
+ double sxx = A * eps[0] + B * trace;
+ double syy = A * eps[1] + B * trace;
+ double szz = A * eps[2] + B * trace;
+ double sxy = A * eps[3];
+ double sxz = A * eps[4];
+ double syz = A * eps[5];
+
+ std::vector<double> vec(9);
+ vec[0]=sxx; vec[1]=sxy; vec[2]=sxz; vec[3]=sxy; vec[4]=syy; vec[5]=syz; vec[6]=sxz; vec[7]=syz; vec[8]=szz;
+
+ data[e->getNum()]=vec;
+ }
+ }
+ PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
+ return pv;
+
+}
+
PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFileName)
{
std::cout << "build Lagrange Multiplier View"<< std::endl;
@@ -568,6 +633,7 @@ PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFile
return pv;
}
+
PView *elasticitySolver::buildElasticEnergyView(const std::string postFileName)
{
std::cout << "build Elastic Energy View"<< std::endl;
@@ -599,7 +665,8 @@ PView *elasticitySolver::buildElasticEnergyView(const std::string postFileName)
}
PView *elasticitySolver::buildVonMisesView(const std::string postFileName)
-{std::cout << "build elastic view"<< std::endl;
+{
+ std::cout << "build elastic view"<< std::endl;
std::map<int, std::vector<double> > data;
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); ++i)
@@ -672,6 +739,10 @@ void elasticitySolverRegisterBindings(binding *b)
cm->setDescription ("assembles the problem");
cm->setArgNames ("linearSystem",NULL);
+ cm = cb->addMethod("solve", &elasticitySolver::solve);
+ cm->setDescription ("solve the problem");
+ cm->setArgNames (NULL);
+
cm = cb->addMethod("addDirichletBC", &elasticitySolver::addDirichletBC);
cm->setDescription ("add a Dirichlet (displacement) boundary condition on a given entity");
cm->setArgNames ("dim", "entityId", "component", "value", NULL);
@@ -702,6 +773,10 @@ void elasticitySolverRegisterBindings(binding *b)
cm->setDescription ("create a new view.");
cm->setArgNames ("fileName", NULL);
+ cm = cb->addMethod ("buildStressesView", &elasticitySolver::buildStressesView);
+ cm->setDescription ("create a new view.");
+ cm->setArgNames ("fileName", NULL);
+
cm = cb->addMethod ("buildLagrangeMultiplierView", &elasticitySolver::buildLagrangeMultiplierView);
cm->setDescription ("create a new view.");
cm->setArgNames ("fileName", NULL);
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 6e1b59747119d1c4449c4956c5bffb891efe2238..2b5c839ff9a7cc247e1963bd77125715906e607f 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -12,6 +12,7 @@
#include "dofManager.h"
#include "simpleFunction.h"
#include "functionSpace.h"
+#include "function.h"
class GModel;
class PView;
@@ -86,6 +87,7 @@ class elasticitySolver
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);
+ void addNeumannBCFct (int dim, int entityId, const function* luaFunction, lua_State *L);
#endif
@@ -103,6 +105,7 @@ class elasticitySolver
void postSolve();
void getSolutionOnElement(MElement *el, fullMatrix<double> &sol);
virtual PView *buildDisplacementView(const std::string postFileName);
+ virtual PView *buildStressesView(const std::string postFileName);
virtual PView *buildLagrangeMultiplierView(const std::string posFileName);
virtual PView *buildElasticEnergyView(const std::string postFileName);
virtual PView *buildVonMisesView(const std::string postFileName);