diff --git a/CMakeLists.txt b/CMakeLists.txt
index ddcc217f08ae130239b9740285a101aa1b625db4..6f491a28dbe9be5cef6507ad4739d19ef81684fb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -109,6 +109,7 @@ set(GMSH_API
Numeric/Numeric.h Numeric/GaussIntegration.h Numeric/polynomialBasis.h
Numeric/JacobianBasis.h Numeric/MetricBasis.h Numeric/bezierBasis.h Numeric/fullMatrix.h
Numeric/simpleFunction.h Numeric/cartesian.h Numeric/ElementType.h
+ Numeric/BasisFactory.h
Geo/GModel.h Geo/GEntity.h Geo/GPoint.h Geo/GVertex.h Geo/GEdge.h
Geo/GFace.h Geo/GRegion.h Geo/GEdgeLoop.h Geo/GEdgeCompound.h
Geo/GFaceCompound.h Geo/GRegionCompound.h Geo/GRbf.h Geo/MVertex.h
@@ -131,6 +132,7 @@ set(GMSH_API
Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h
Solver/linearSystemFull.h Solver/elasticitySolver.h Solver/sparsityPattern.h
Solver/groupOfElements.h Solver/linearSystemPETSc.h Solver/linearSystemMUMPS.h
+ Solver/thermicSolver.h
Post/PView.h Post/PViewData.h Plugin/PluginManager.h Post/OctreePost.h
Post/PViewDataList.h Post/PViewDataGModel.h Post/PViewOptions.h Post/ColorTable.h
Numeric/nodalBasis.h
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index a4d1657167df2d0bdbabf1e03d90118e35128a66..4b14fb8f1ed1c990e27cd4fd46d8d29c0460eebc 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -12,6 +12,7 @@ set(SRC
groupOfElements.cpp
elasticityTerm.cpp
elasticitySolver.cpp
+ thermicSolver.cpp
SElement.cpp
eigenSolver.cpp
multiscaleLaplace.cpp
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 58d1aa06968da333c12de9b687d668ee6a869e69..610785068bb2be0beb8ea243d57a98fc401b40ea 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -119,12 +119,11 @@ void elasticitySolver::solve()
// printLinearSystem(lsys);
double energ=0;
- for (unsigned int i = 0; i < elasticFields.size(); i++)
- {
+ for (unsigned int i = 0; i < elasticFields.size(); i++){
SolverField<SVector3> Field(pAssembler, LagSpace);
IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
- Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),
+ Assemble(Elastic_Energy_Term, elasticFields[i].g->begin(), elasticFields[i].g->end(),
Integ_Bulk,energ);
}
printf("elastic energy=%f\n",energ);
@@ -135,8 +134,7 @@ void elasticitySolver::postSolve()
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
double energ=0;
- for (unsigned int i = 0; i < elasticFields.size(); i++)
- {
+ for (unsigned int i = 0; i < elasticFields.size(); i++){
SolverField<SVector3> Field(pAssembler, LagSpace);
IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
@@ -156,12 +154,12 @@ void elasticitySolver::readInputFile(const std::string &fn)
return;
}
if(what[0]=='#'){
- /*
+
char *line=NULL;
size_t l = 0;
- int r = getline(&line,&l,f);
+ getline(&line,&l,f);
free(line);
- */
+
}
else if (!strcmp(what, "ElasticDomain")){
elasticField field;
@@ -185,7 +183,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
}
field._tag = _tag;
field._d = SVector3(d1, d2, d3);
- field._f = simpleFunction<double>(val);
+ field._f = new simpleFunction<double>(val);
field.g = new groupOfElements (_dim - 1, physical);
LagrangeMultiplierFields.push_back(field);
}
@@ -319,6 +317,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
int tag=1;
gLevelsetPlane ls(a,b,c,d,tag);
pModel = pModel->buildCutGModel(&ls);
+ pModel->writeMSH("cutMesh.msh");
}
else if (!strcmp(what, "CutMeshSphere")){
double x, y, z, r;
@@ -329,6 +328,19 @@ void elasticitySolver::readInputFile(const std::string &fn)
int tag=1;
gLevelsetSphere ls(x,y,z,r,tag);
pModel = pModel->buildCutGModel(&ls);
+ pModel->writeMSH("cutMesh.msh");
+ }
+ else if (!strcmp(what, "CutMeshMathExpr")){
+ char expr[256];
+ if(fscanf(f, "%s", expr) != 1){
+ fclose(f);
+ return;
+ }
+ std::string exprS(expr);
+ int tag=1;
+ gLevelsetMathEval ls(exprS,tag);
+ pModel = pModel->buildCutGModel(&ls);
+ pModel->writeMSH("cutMesh.msh");
}
else {
Msg::Error("Invalid input : '%s'", what);
@@ -338,7 +350,8 @@ void elasticitySolver::readInputFile(const std::string &fn)
}
-void elasticitySolver::addDirichletBC (int dim, int entityId, int component, double value) {
+void elasticitySolver::addDirichletBC(int dim, int entityId, int component, double value)
+{
dirichletBC diri;
diri.g = new groupOfElements (dim, entityId);
diri._f= new simpleFunction<double>(value);
@@ -353,8 +366,14 @@ void elasticitySolver::addDirichletBC (int dim, int entityId, int component, dou
allDirichlet.push_back(diri);
}
+void elasticitySolver::addDirichletBC(int dim, std::string phys, int component, double value)
+{
+ int entityId = pModel->getPhysicalNumber(dim, phys);
+ addDirichletBC(dim, entityId, component, value);
+}
-void elasticitySolver::addNeumannBC (int dim, int entityId, const std::vector<double> value) {
+void elasticitySolver::addNeumannBC(int dim, int entityId, const std::vector<double> value)
+{
if(value.size()!=3) return;
neumannBC neu;
neu.g = new groupOfElements (dim, entityId);
@@ -369,9 +388,14 @@ void elasticitySolver::addNeumannBC (int dim, int entityId, const std::vector<do
allNeumann.push_back(neu);
}
+void elasticitySolver::addNeumannBC(int dim, std::string phys, const std::vector<double> value)
+{
+ int entityId = pModel->getPhysicalNumber(dim, phys);
+ addNeumannBC(dim, entityId, value);
+}
-
-void elasticitySolver::addElasticDomain (int physical, double e, double nu){
+void elasticitySolver::addElasticDomain(int physical, double e, double nu)
+{
elasticField field;
field._tag = _tag;
field._E = e;
@@ -380,6 +404,12 @@ void elasticitySolver::addElasticDomain (int physical, double e, double nu){
elasticFields.push_back(field);
}
+void elasticitySolver::addElasticDomain(std::string phys, double e, double nu)
+{
+ int entityId = pModel->getPhysicalNumber(_dim, phys);
+ addElasticDomain(entityId, e, nu);
+}
+
elasticitySolver::elasticitySolver(GModel *model, int tag)
{
pModel = model;
@@ -387,7 +417,9 @@ elasticitySolver::elasticitySolver(GModel *model, int tag)
_tag = tag;
pAssembler = NULL;
if (_dim==3) LagSpace=new VectorLagrangeFunctionSpace(_tag);
- if (_dim==2) LagSpace=new VectorLagrangeFunctionSpace(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+ if (_dim==2) LagSpace=new VectorLagrangeFunctionSpace(_tag,
+ VectorLagrangeFunctionSpace::VECTOR_X,
+ VectorLagrangeFunctionSpace::VECTOR_Y);
LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpace(_tag+1);
}
@@ -401,28 +433,24 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
// numbered afterwards
// Dirichlet conditions
- for (unsigned int i = 0; i < allDirichlet.size(); i++)
- {
+ for (unsigned int i = 0; i < allDirichlet.size(); i++){
FilterDofComponent filter(allDirichlet[i]._comp);
- FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),
- *pAssembler,*allDirichlet[i]._f,filter);
+ FixNodalDofs(*LagSpace, allDirichlet[i].g->begin(), allDirichlet[i].g->end(),
+ *pAssembler, *allDirichlet[i]._f, filter);
}
// LagrangeMultipliers
- for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i)
- {
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i){
NumberDofs(*LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(),
LagrangeMultiplierFields[i].g->end(), *pAssembler);
}
// Elastic Fields
- for (unsigned int i = 0; i < elasticFields.size(); ++i)
- {
+ for (unsigned int i = 0; i < elasticFields.size(); ++i){
if(elasticFields[i]._E != 0.)
NumberDofs(*LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),
*pAssembler);
}
// Voids
- for (unsigned int i = 0; i < elasticFields.size(); ++i)
- {
+ for (unsigned int i = 0; i < elasticFields.size(); ++i){
if(elasticFields[i]._E == 0.)
FixVoidNodalDofs(*LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),
*pAssembler);
@@ -430,19 +458,17 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
// Neumann conditions
GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
- for (unsigned int i = 0; i < allNeumann.size(); i++)
- {
- LoadTerm<SVector3> Lterm(*LagSpace,*allNeumann[i]._f);
- Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),
+ for (unsigned int i = 0; i < allNeumann.size(); i++){
+ LoadTerm<SVector3> Lterm(*LagSpace, allNeumann[i]._f);
+ Assemble(Lterm, *LagSpace, allNeumann[i].g->begin(), allNeumann[i].g->end(),
Integ_Boundary,*pAssembler);
}
// Assemble cross term, laplace term and rhs term for LM
GaussQuadrature Integ_LagrangeMult(GaussQuadrature::ValVal);
GaussQuadrature Integ_Laplace(GaussQuadrature::GradGrad);
- for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); i++)
- {
- LagrangeMultiplierTerm LagTerm(*LagSpace, *LagrangeMultiplierSpace,
- LagrangeMultiplierFields[i]._d);
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); i++){
+ LagrangeMultiplierTerm<SVector3> LagTerm(*LagSpace, *LagrangeMultiplierSpace,
+ LagrangeMultiplierFields[i]._d);
Assemble(LagTerm, *LagSpace, *LagrangeMultiplierSpace,
LagrangeMultiplierFields[i].g->begin(),
LagrangeMultiplierFields[i].g->end(), Integ_LagrangeMult, *pAssembler);
@@ -458,11 +484,10 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
}
// Assemble elastic term for
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
- for (unsigned int i = 0; i < elasticFields.size(); i++)
- {
+ for (unsigned int i = 0; i < elasticFields.size(); i++){
IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
- Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),
- Integ_Bulk,*pAssembler);
+ Assemble(Eterm, *LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),
+ Integ_Bulk, *pAssembler);
}
/*for (int i=0;i<pAssembler->sizeOfR();i++){
@@ -479,14 +504,9 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
}
-void elasticitySolver::getSolutionOnElement (MElement *el, fullMatrix<double> &sol) {
-
-}
-
-#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];
@@ -512,27 +532,10 @@ 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 valx[256];
- double valy[256];
- double valz[256];
- for (int k = 0; k < e->getNumVertices(); k++){
- a->getDofValue(e->getVertex(k), 0, _tag, valx[k]);
- a->getDofValue(e->getVertex(k), 1, _tag, valy[k]);
- a->getDofValue(e->getVertex(k), 2, _tag, valz[k]);
- }
- double gradux[3];
- double graduy[3];
- double graduz[3];
- 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 E, double nu, int _tag)
+{
+ double eps[6];
+ deformation(a, e, u, v, w, _tag, eps);
double A = E / (1. + nu);
double B = A * (nu / (1. - 2 * nu));
double trace = eps[0] + eps[1] + eps[2] ;
@@ -547,19 +550,123 @@ static double vonMises(dofManager<double> *a, MElement *e,
return ComputeVonMises(s);
}
-*/
+
+void elasticitySolver::computeEffectiveStiffness(std::vector<double> stiff)
+{
+ double st[6] = {0., 0., 0., 0., 0., 0.};
+ double volTot = 0.;
+ for (unsigned int i = 0; i < elasticFields.size(); ++i){
+ double E = elasticFields[i]._E;
+ double nu = elasticFields[i]._nu;
+ SolverField<SVector3> Field(pAssembler, LagSpace);
+ for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
+ it != elasticFields[i].g->end(); ++it){
+ MElement *e = *it;
+ double vol = e->getVolume() * e->getVolumeSign();
+ 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];
+ SPoint3 center = e->barycenterUVW();
+ double u = center.x(), v = center.y(), w = center.z();
+ 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] ;
+ st[0] += (A * eps[0] + B * trace) * vol;
+ st[1] += (A * eps[1] + B * trace) * vol;
+ st[2] += (A * eps[2] + B * trace) * vol;
+ st[3] += (A * eps[3]) * vol;
+ st[4] += (A * eps[4]) * vol;
+ st[5] += (A * eps[5]) * vol;
+ volTot += vol;
+ }
+ }
+ for(int i = 0; i < 6;i++)
+ stiff[i] = st[i] / volTot;
+}
+
+void elasticitySolver::computeEffectiveStrain(std::vector<double> strain)
+{
+ double st[6] = {0., 0., 0., 0., 0., 0.};
+ double volTot = 0.;
+ for (unsigned int i = 0; i < elasticFields.size(); ++i){
+ SolverField<SVector3> Field(pAssembler, LagSpace);
+ for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
+ it != elasticFields[i].g->end(); ++it){
+ MElement *e = *it;
+ double vol = e->getVolume() * e->getVolumeSign();
+ 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];
+ SPoint3 center = e->barycenterUVW();
+ double u = center.x(), v = center.y(), w = center.z();
+ e->interpolateGrad(valx, u, v, w, gradux);
+ e->interpolateGrad(valy, u, v, w, graduy);
+ e->interpolateGrad(valz, u, v, w, graduz);
+
+ st[0] += gradux[0] * vol;
+ st[1] += graduy[1] * vol;
+ st[2] += graduz[2] * vol;
+ st[3] += 0.5 * (gradux[1] + graduy[0]) * vol;
+ st[4] += 0.5 * (gradux[2] + graduz[0]) * vol;
+ st[5] += 0.5 * (graduy[2] + graduz[1]) * vol;
+ volTot += vol;
+ }
+ }
+ for(int i = 0; i < 6;i++)
+ strain[i] = st[i] / volTot;
+}
+
+#if defined(HAVE_POST)
PView* elasticitySolver::buildDisplacementView (const std::string postFileName)
{
std::cout << "build Displacement View"<< std::endl;
std::set<MVertex*> v;
- std::map<MVertex*,MElement*> vCut;
- for (unsigned int i = 0; i < elasticFields.size(); ++i)
- {
+ std::map<MVertex*, MElement*> vCut;
+ for (unsigned int i = 0; i < elasticFields.size(); ++i){
if(elasticFields[i]._E == 0.) continue;
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
it != elasticFields[i].g->end(); ++it){
- MElement *e=*it;
+ MElement *e = *it;
if(e->getParent()) {
for (int j = 0; j < e->getNumVertices(); ++j) {
if(vCut.find(e->getVertex(j)) == vCut.end())
@@ -577,18 +684,20 @@ PView* elasticitySolver::buildDisplacementView (const std::string postFileName)
for (std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it){
SVector3 val;
MPoint p(*it);
- Field.f(&p,0,0,0,val);
- std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2);
- data[(*it)->getNum()]=vec;
+ Field.f(&p, 0, 0, 0, val);
+ std::vector<double> vec(3);
+ vec[0] = val(0); vec[1] = val(1); vec[2] = val(2);
+ data[(*it)->getNum()] = vec;
}
for (std::map<MVertex*,MElement*>::iterator it = vCut.begin(); it != vCut.end(); ++it){
SVector3 val;
double uvw[3];
double xyz[3] = {it->first->x(), it->first->y(), it->first->z()};
it->second->xyz2uvw(xyz, uvw);
- Field.f(it->second,uvw[0],uvw[1],uvw[2],val);
- std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2);
- data[it->first->getNum()]=vec;
+ Field.f(it->second, uvw[0], uvw[1], uvw[2], val);
+ std::vector<double> vec(3);
+ vec[0] = val(0); vec[1] = val(1); vec[2] = val(2);
+ data[it->first->getNum()] = vec;
}
PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
return pv;
@@ -596,20 +705,19 @@ PView* elasticitySolver::buildDisplacementView (const std::string postFileName)
PView* elasticitySolver::buildStressesView (const std::string postFileName)
{
+ double sti[6] = {0., 0., 0., 0., 0., 0.};
+ double str[6] = {0., 0., 0., 0., 0., 0.};
+ double volTot = 0.;
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)
- {
+ 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 Elastic_Energy_Term(Eterm);
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
- it != elasticFields[i].g->end(); ++it)
- {
- MElement *e=*it;
+ it != elasticFields[i].g->end(); ++it){
+ MElement *e = *it;
+ double vol = e->getVolume() * e->getVolumeSign();
int nbVertex = e->getNumVertices();
std::vector<SVector3> val(nbVertex);
@@ -619,16 +727,17 @@ PView* elasticitySolver::buildStressesView (const std::string postFileName)
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);
+ 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;
+ SPoint3 center = e->barycenterUVW();
+ double u = center.x(), v = center.y(), w = center.z();
e->interpolateGrad(valx, u, v, w, gradux);
e->interpolateGrad(valy, u, v, w, graduy);
e->interpolateGrad(valz, u, v, w, graduz);
@@ -638,7 +747,6 @@ PView* elasticitySolver::buildStressesView (const std::string postFileName)
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] ;
@@ -650,15 +758,79 @@ PView* elasticitySolver::buildStressesView (const std::string postFileName)
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;
+ 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;
+
+ for(int k = 0; k < 6; k++)
+ str[k] += eps[k] * vol;
+ sti[0] += sxx * vol;
+ sti[1] += syy * vol;
+ sti[2] += szz * vol;
+ sti[3] += sxy * vol;
+ sti[4] += sxz * vol;
+ sti[5] += syz * vol;
+ volTot += vol;
}
}
+ for(int i = 0; i < 6; i++){
+ str[i] = str[i] / volTot;
+ sti[i] = sti[i] / volTot;
+ }
+ printf("effective stiffn = ");for(int i = 0; i < 6; i++) printf("%g ",sti[i]);printf("\n");
+ printf("effective strain = ");for(int i = 0; i < 6; i++) printf("%g ",str[i]);printf("\n");
+
PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
return pv;
+}
+
+PView* elasticitySolver::buildStrainView (const std::string postFileName)
+{
+ std::cout << "build strain view"<< std::endl;
+ std::map<int, std::vector<double> > data;
+ for (unsigned int i = 0; i < elasticFields.size(); ++i){
+ SolverField<SVector3> Field(pAssembler, LagSpace);
+ 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);
+ std::vector<double> vec(9);
+ vec[0] = gradux[0];
+ vec[4] = graduy[1];
+ vec[8] = graduy[2];
+ vec[1] = vec[3] = 0.5 * (gradux[0] + graduy[1]);
+ vec[2] = vec[6] = 0.5 * (gradux[0] + graduz[2]);
+ vec[5] = vec[7] = 0.5 * (gradux[1] + graduz[2]);
+
+ data[e->getNum()] = vec;
+ }
+ }
+ PView *pv = new PView(postFileName, "ElementData", pModel, data, 0.0);
+ return pv;
}
PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFileName)
@@ -666,20 +838,17 @@ PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFile
std::cout << "build Lagrange Multiplier View"<< std::endl;
if(!LagrangeMultiplierSpace) return new PView();
std::set<MVertex*> v;
- for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i)
- {
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i){
for(groupOfElements::elementContainer::const_iterator it =
LagrangeMultiplierFields[i].g->begin();
- it != LagrangeMultiplierFields[i].g->end(); ++it)
- {
+ it != LagrangeMultiplierFields[i].g->end(); ++it){
MElement *e = *it;
for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
}
}
std::map<int, std::vector<double> > data;
SolverField<double> Field(pAssembler, LagrangeMultiplierSpace);
- for(std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
- {
+ for(std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it){
double val;
MPoint p(*it);
Field.f(&p, 0, 0, 0, val);
@@ -697,60 +866,95 @@ PView *elasticitySolver::buildElasticEnergyView(const std::string postFileName)
std::cout << "build Elastic Energy View"<< std::endl;
std::map<int, std::vector<double> > data;
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
- for (unsigned int i = 0; i < elasticFields.size(); ++i)
- {
+ for (unsigned int i = 0; i < elasticFields.size(); ++i){
if(elasticFields[i]._E == 0.) continue;
SolverField<SVector3> Field(pAssembler, LagSpace);
- IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+ IsotropicElasticTerm Eterm(Field, elasticFields[i]._E, elasticFields[i]._nu);
BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
ScalarTermConstant<double> One(1.0);
for (groupOfElements::elementContainer::const_iterator it =
- elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
- {
+ elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it){
MElement *e = *it;
double energ;
double vol;
IntPt *GP;
- int npts=Integ_Bulk.getIntPoints(e,&GP);
- Elastic_Energy_Term.get(e,npts,GP,energ);
- One.get(e,npts,GP,vol);
+ int npts = Integ_Bulk.getIntPoints(e, &GP);
+ Elastic_Energy_Term.get(e, npts, GP, energ);
+ One.get(e, npts, GP, vol);
std::vector<double> vec;
- vec.push_back(energ/vol);
- data[e->getNum()]=vec;
+ vec.push_back(energ / vol);
+ data[e->getNum()] = vec;
}
}
PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
return pv;
}
+PView *elasticitySolver::buildVolumeView(const std::string postFileName)
+{
+ std::cout << "build Volume View";
+ std::map<int, std::vector<double> > data;
+ double voltot = 0; double length = 0;
+ GaussQuadrature Integ_Vol(GaussQuadrature::Val);
+ for (unsigned int i = 0; i < elasticFields.size(); ++i){
+ ScalarTermConstant<double> One(1.0);
+ for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
+ it != elasticFields[i].g->end(); ++it){
+ MElement *e = *it;
+ double vol;
+ IntPt *GP;
+ int npts = Integ_Vol.getIntPoints(e, &GP);
+ One.get(e, npts, GP, vol);
+ voltot += vol;
+ std::vector<double> vec;
+ vec.push_back(vol);
+ data[e->getNum()] = vec;
+ }
+ }
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i){
+ ScalarTermConstant<double> One(1.0);
+ for (groupOfElements::elementContainer::const_iterator it = LagrangeMultiplierFields[i].g->begin();
+ it != LagrangeMultiplierFields[i].g->end(); ++it){
+ MElement *e = *it;
+ double l;
+ IntPt *GP;
+ int npts = Integ_Vol.getIntPoints(e, &GP);
+ One.get(e, npts, GP, l);
+ length += l;
+ }
+ std::cout << " : length " << LagrangeMultiplierFields[i]._tag << " = " << length;
+ length = 0;
+ }
+ PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0, 1);
+ std::cout << " / total vol = " << voltot << std::endl;
+ return pv;
+}
+
PView *elasticitySolver::buildVonMisesView(const std::string postFileName)
{
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)
- {
+ for (unsigned int i = 0; i < elasticFields.size(); ++i){
SolverField<SVector3> Field(pAssembler, LagSpace);
- IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+ IsotropicElasticTerm Eterm(Field, elasticFields[i]._E, elasticFields[i]._nu);
BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
- it != elasticFields[i].g->end(); ++it)
- {
- MElement *e=*it;
+ it != elasticFields[i].g->end(); ++it){
+ MElement *e = *it;
double energ;
IntPt *GP;
- int npts=Integ_Bulk.getIntPoints(e,&GP);
- Elastic_Energy_Term.get(e,npts,GP,energ);
+ int npts = Integ_Bulk.getIntPoints(e, &GP);
+ Elastic_Energy_Term.get(e, npts, GP, energ);
std::vector<double> vec;
vec.push_back(energ);
- data[(*it)->getNum()]=vec;
+ data[(*it)->getNum()] = vec;
}
}
PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
return pv;
}
-
#else
PView* elasticitySolver::buildDisplacementView(const std::string postFileName)
@@ -759,6 +963,12 @@ PView* elasticitySolver::buildDisplacementView(const std::string postFileName)
return 0;
}
+PView* elasticitySolver::buildStrainView(const std::string postFileName)
+{
+ Msg::Error("Post-pro module not available");
+ return 0;
+}
+
PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFileName)
{
Msg::Error("Post-pro module not available");
@@ -782,5 +992,10 @@ PView* elasticitySolver::buildStressesView (const std::string postFileName)
Msg::Error("Post-pro module not available");
return 0;
}
+PView *elasticitySolver::buildVolumeView(const std::string postFileName)
+{
+ Msg::Error("Post-pro module not available");
+ return 0;
+}
#endif
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index cfb879139b276126a7fd3ed190af93b807f84b85..423e5ff400dcc5cae0f1cadcadc809004209a991 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -22,7 +22,7 @@ struct LagrangeMultiplierField {
groupOfElements *g;
double _tau;
SVector3 _d;
- simpleFunction<double> _f;
+ simpleFunction<double> *_f;
LagrangeMultiplierField() : _tag(0), g(0){}
};
@@ -78,9 +78,12 @@ class elasticitySolver
elasticitySolver(GModel *model, int tag);
- void addDirichletBC (int dim, int entityId, int component, double value);
- void addNeumannBC (int dim, int entityId, const std::vector<double> value);
- void addElasticDomain (int tag, double e, double nu);
+ void addDirichletBC(int dim, int entityId, int component, double value);
+ void addDirichletBC(int dim, std::string phys, int component, double value);
+ void addNeumannBC(int dim, int entityId, const std::vector<double> value);
+ void addNeumannBC(int dim, std::string phys, const std::vector<double> value);
+ void addElasticDomain(int tag, double e, double nu);
+ void addElasticDomain(std::string phys, double e, double nu);
virtual ~elasticitySolver()
{
@@ -95,12 +98,15 @@ class elasticitySolver
void solve();
void postSolve();
void exportKb();
- void getSolutionOnElement(MElement *el, fullMatrix<double> &sol);
+ void computeEffectiveStiffness(std::vector<double> stiff);
+ void computeEffectiveStrain(std::vector<double> strain);
virtual PView *buildDisplacementView(const std::string postFileName);
+ virtual PView *buildStrainView(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);
+ virtual PView *buildVolumeView(const std::string postFileName);
// std::pair<PView *, PView*> buildErrorEstimateView
// (const std::string &errorFileName, double, int);
// std::pair<PView *, PView*> buildErrorEstimateView
diff --git a/Solver/femTerm.h b/Solver/femTerm.h
index 85d8ccd94abd9cee987367a7113746bdc92f312e..584b5be3d7fa6f40b5c9984859d4ba39ca83b861 100644
--- a/Solver/femTerm.h
+++ b/Solver/femTerm.h
@@ -54,7 +54,7 @@ class femTerm {
groupOfElements &C) const
{
groupOfElements::elementContainer::const_iterator it = L.begin();
- for ( ; it != L.end() ; ++it){
+ for ( ; it != L.end(); ++it){
MElement *eL = *it;
if (&C == &L || C.find(eL)){
SElement se(eL);
@@ -78,7 +78,7 @@ class femTerm {
{
const int nbR = localMatrix.size1();
const int nbC = localMatrix.size2();
- std::vector<Dof> R,C; // better use default consdtructors and reserve the right amount of space to avoid reallocation
+ std::vector<Dof> R, C; // better use default consdtructors and reserve the right amount of space to avoid reallocation
R.reserve(nbR);
C.reserve(nbC);
bool sym=true;
@@ -90,12 +90,12 @@ class femTerm {
Dof c(getLocalDofC(se, j));
R.push_back(r);
C.push_back(c);
- if (!(r==c)) sym=false;
+ if (!(r == c)) sym = false;
}
}
else
{
- sym=false;
+ sym = false;
for (int j = 0; j < nbR; j++)
R.push_back(getLocalDofR(se, j));
for (int k = 0; k < nbC; k++)
@@ -118,6 +118,31 @@ class femTerm {
dm.fixVertex(v[i], comp, field, e(v[i]->x(), v[i]->y(), v[i]->z()));
}
+ void neumannNodalBC(MElement *e, int comp, int field,
+ const simpleFunction<dataVec> &fct,
+ dofManager<dataVec> &dm)
+ {
+ double jac[3][3];
+ double sf[256];
+ int integrationOrder = 2 * e->getPolynomialOrder();
+ int npts;
+ IntPt *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];
+ const double w = GP[ip].pt[2];
+ const double weight = GP[ip].weight;
+ const double detJ = e->getJacobian(u, v, w, jac);
+ SPoint3 p; e->pnt(u, v, w, p);
+ e->getShapeFunctions(u, v, w, sf);
+ const dataVec FCT = fct(p.x(), p.y(), p.z());
+ for (int k = 0; k < e->getNumShapeFunctions(); k++){
+ dm.assemble(e->getShapeFunctionNode(k), comp, field, detJ * weight * sf[k] * FCT);
+ }
+ }
+ }
+
void neumannNodalBC(int physical, int dim, int comp, int field,
const simpleFunction<dataVec> &fct,
dofManager<dataVec> &dm)
@@ -127,30 +152,31 @@ class femTerm {
m->getPhysicalGroups(groups);
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];
for (unsigned int i = 0; i < it->second.size(); ++i){
GEntity *ge = it->second[i];
for (unsigned int j = 0; j < ge->getNumMeshElements(); j++){
MElement *e = ge->getMeshElement(j);
- int integrationOrder = 2 * e->getPolynomialOrder();
- int nbNodes = e->getNumVertices();
- int npts;
- IntPt *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];
- const double w = GP[ip].pt[2];
- const double weight = GP[ip].weight;
- const double detJ = e->getJacobian(u, v, w, jac);
- SPoint3 p; e->pnt(u, v, w, p);
- e->getShapeFunctions(u, v, w, sf);
- const dataVec FCT = fct(p.x(), p.y(), p.z());
- for (int k = 0; k < nbNodes; k++){
- dm.assemble(e->getVertex(k), comp, field, detJ * weight * sf[k] * FCT);
- }
- }
+ neumannNodalBC(e, comp, field, fct, dm);
+ }
+ }
+ }
+ void neumannNormalNodalBC(int physical, int dim, int field,
+ const simpleFunction<dataVec> &fct,
+ dofManager<dataVec> &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);
+ if (it == groups[dim].end()) return;
+ for (unsigned int i = 0; i < it->second.size(); ++i){
+ GEntity *ge = it->second[i];
+ for (unsigned int j = 0; j < ge->getNumMeshElements(); j++){
+ MElement *e = ge->getMeshElement(j);
+
+ neumannNodalBC(e, 0, field, fct, dm);
+ neumannNodalBC(e, 1, field, fct, dm);
+ neumannNodalBC(e, 2, field, fct, dm);
}
}
}
@@ -182,8 +208,8 @@ class DummyfemTerm : public femTerm<double>
private : // i dont want to mess with this anymore
virtual int sizeOfC(SElement *se) const {return 0;}
virtual int sizeOfR(SElement *se) const {return 0;}
- virtual Dof getLocalDofR(SElement *se, int iRow) const {return Dof(0,0);}
- virtual Dof getLocalDofC(SElement *se, int iCol) const {return Dof(0,0);}
+ virtual Dof getLocalDofR(SElement *se, int iRow) const {return Dof(0, 0);}
+ virtual Dof getLocalDofC(SElement *se, int iCol) const {return Dof(0, 0);}
virtual void elementMatrix(SElement *se, fullMatrix<dataMat> &m) const {m.scale(0.);}
virtual void elementVector(SElement *se, fullVector<dataVec> &m) const {m.scale(0.);}
};
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index d890cd09f9dfbdee3524b6b5ccfd8c6664bb76cf..a3f356ae9015e0f08c5f878ffa532d259103c898 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -115,7 +115,8 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
if(ele->getParent()) {
- if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+ if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B ||
+ ele->getTypeForMSH() == MSH_POLYG_B) { //FIXME MPolygonBorders...
ele->movePointFromParentSpaceToElementSpace(u, v, w);
}
}
@@ -128,7 +129,8 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
{
if(ele->getParent()) {
- if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+ if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B ||
+ ele->getTypeForMSH() == MSH_POLYG_B) { //FIXME MPolygonBorders...
ele->movePointFromParentSpaceToElementSpace(u, v, w);
}
}
@@ -150,7 +152,8 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
virtual void hessfuvw(MElement *ele, double u, double v, double w, std::vector<HessType> &hess)
{
if(ele->getParent()) {
- if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+ if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B ||
+ ele->getTypeForMSH() == MSH_POLYG_B) { //FIXME MPolygonBorders...
ele->movePointFromParentSpaceToElementSpace(u, v, w);
}
}
@@ -169,7 +172,8 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
virtual void gradfuvw(MElement *ele, double u, double v, double w, std::vector<GradType> &grads)
{
if(ele->getParent()) {
- if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B) { //FIXME MPolygonBorders...
+ if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B ||
+ ele->getTypeForMSH() == MSH_POLYG_B) { //FIXME MPolygonBorders...
ele->movePointFromParentSpaceToElementSpace(u, v, w);
}
}
diff --git a/Solver/terms.cpp b/Solver/terms.cpp
index 97f2bf5fc73f54f8810ed88dcc102e31a434c144..cff7fbb8aa25edf1dcdbcf8ab46ebcab5e7db875 100644
--- a/Solver/terms.cpp
+++ b/Solver/terms.cpp
@@ -142,31 +142,6 @@ void IsotropicElasticTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<do
}
}
-void LagrangeMultiplierTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
-{
- int nbFF1 = BilinearTerm<SVector3, double>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
- int nbFF2 = BilinearTerm<SVector3, double>::space2.getNumKeys(ele); //nbVertices of boundary
- double jac[3][3];
- m.resize(nbFF1, nbFF2);
- m.setAll(0.);
- for(int i = 0; i < npts; i++)
- {
- double u = GP[i].pt[0]; double v = GP[i].pt[1]; double w = GP[i].pt[2];
- const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
- std::vector<TensorialTraits<SVector3>::ValType> Vals;
- std::vector<TensorialTraits<double>::ValType> ValsT;
- BilinearTerm<SVector3,double>::space1.f(ele, u, v, w, Vals);
- BilinearTerm<SVector3,double>::space2.f(ele, u, v, w, ValsT);
- for(int j = 0; j < nbFF1; j++)
- {
- for(int k = 0; k < nbFF2; k++)
- {
- m(j, k) += dot(Vals[j], _d) * ValsT[k] * weight * detJ;
- }
- }
- }
-}
-
void LagMultTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
{
int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
diff --git a/Solver/terms.h b/Solver/terms.h
index 2e0041137d500f6b5411750a283be7bc45abcc2d..b49c75f54dbbebaa56ebe0ee65e84fad711fe1bd 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -34,8 +34,14 @@ inline double dot(const double &a, const double &b)
inline int delta(int i,int j) {if (i==j) return 1; else return 0;}
-
-
+inline void setDirection(double &a, const double &b)
+{
+ a = b;
+}
+inline void setDirection(SVector3 &a, const SVector3 &b)
+{
+ for(int i = 0; i < 3; i++) a(i) = b(i);
+}
template<class T2=double> class ScalarTermBase
{
@@ -234,9 +240,9 @@ public :
template<class T1> class LoadTerm : public LinearTerm<T1>
{
protected:
- simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
+ simpleFunction<typename TensorialTraits<T1>::ValType> *Load;
public :
- LoadTerm(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> &Load_) :
+ LoadTerm(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> *Load_) :
LinearTerm<T1>(space1_), Load(Load_) {}
virtual ~LoadTerm() {}
virtual LinearTermBase<double>* clone () const { return new LoadTerm<T1>(LinearTerm<T1>::space1,Load);}
@@ -244,19 +250,19 @@ public :
virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<double> > &vv) const {};
};
-class LagrangeMultiplierTerm : public BilinearTerm<SVector3, double>
+template<class T1> class LagrangeMultiplierTerm : public BilinearTerm<T1, double>
{
- SVector3 _d;
+ T1 _d;
public :
- LagrangeMultiplierTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<double>& space2_, const SVector3 &d) :
- BilinearTerm<SVector3, double>(space1_, space2_)
+ LagrangeMultiplierTerm(FunctionSpace<T1>& space1_, FunctionSpace<double>& space2_, const T1 &d) :
+ BilinearTerm<T1, double>(space1_, space2_)
{
- for (int i = 0; i < 3; i++) _d(i) = d(i);
+ setDirection(_d, d);
}
virtual ~LagrangeMultiplierTerm() {}
virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const;
virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< fullMatrix<double> > &vm) const{};
- virtual BilinearTermBase* clone () const {return new LagrangeMultiplierTerm(BilinearTerm<SVector3, double>::space1,BilinearTerm<SVector3, double>::space2,_d);}
+ virtual BilinearTermBase* clone () const {return new LagrangeMultiplierTerm(BilinearTerm<T1, double>::space1, BilinearTerm<T1, double>::space2, _d);}
};
class LagMultTerm : public BilinearTerm<SVector3, SVector3>
@@ -276,9 +282,9 @@ template<class T1> class LoadTermOnBorder : public LinearTerm<T1>
{
private :
double _eqfac;
- simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
+ simpleFunction<typename TensorialTraits<T1>::ValType> *Load;
public :
- LoadTermOnBorder(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> &Load_, double eqfac = 1.0) :
+ LoadTermOnBorder(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> *Load_, double eqfac = 1.0) :
LinearTerm<T1>(space1_), _eqfac(eqfac), Load(Load_) {}
virtual ~LoadTermOnBorder() {}
virtual LinearTermBase<double>* clone () const { return new LoadTermOnBorder<T1>(LinearTerm<T1>::space1,Load,_eqfac);}
diff --git a/Solver/terms.hpp b/Solver/terms.hpp
index fef362b6eaf5920154b3fded57e8bee4bfab591d..a6b0e3f43affd598b45117ad57886235acf90c4f 100644
--- a/Solver/terms.hpp
+++ b/Solver/terms.hpp
@@ -7,7 +7,7 @@
// Eric Bechet
//
-#include "terms.h"
+//#include "terms.h"
template<class T2> void LinearTermBase<T2>::get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &vec) const
{
@@ -101,7 +101,7 @@ template<class T1> void LoadTerm<T1>::get(MElement *ele, int npts, IntPt *GP, fu
LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
SPoint3 p;
ele->pnt(u, v, w, p);
- typename TensorialTraits<T1>::ValType load = Load(p.x(), p.y(), p.z());
+ typename TensorialTraits<T1>::ValType load = (*Load)(p.x(), p.y(), p.z());
for(int j = 0; j < nbFF ; ++j)
{
m(j) += dot(Vals[j], load) * weight * detJ;
@@ -191,25 +191,49 @@ template<class T1> void GradTerm<T1>::get(MElement *ele, int npts, IntPt *GP, st
}
}
-
+template<class T1> void LagrangeMultiplierTerm<T1>::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+ int nbFF1 = BilinearTerm<T1, double>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
+ int nbFF2 = BilinearTerm<T1, double>::space2.getNumKeys(ele); //nbVertices of boundary
+ double jac[3][3];
+ m.resize(nbFF1, nbFF2);
+ m.setAll(0.);
+ for(int i = 0; i < npts; i++)
+ {
+ double u = GP[i].pt[0]; double v = GP[i].pt[1]; double w = GP[i].pt[2];
+ const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
+ std::vector<typename TensorialTraits<T1>::ValType> Vals;
+ std::vector<TensorialTraits<double>::ValType> ValsT;
+ BilinearTerm<T1,double>::space1.f(ele, u, v, w, Vals);
+ BilinearTerm<T1,double>::space2.f(ele, u, v, w, ValsT);
+ for(int j = 0; j < nbFF1; j++)
+ {
+ for(int k = 0; k < nbFF2; k++)
+ {
+ m(j, k) += dot(Vals[j], _d) * ValsT[k] * weight * detJ;
+ }
+ }
+ }
+}
template<class T1> void LoadTermOnBorder<T1>::get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m) const
{
- MElement *elep;
- if (ele->getParent()) elep = ele->getParent();
int nbFF = LinearTerm<T1>::space1.getNumKeys(ele);
double jac[3][3];
m.resize(nbFF);
m.scale(0.);
for(int i = 0; i < npts; i++)
{
- const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
+ double u = GP[i].pt[0]; double v = GP[i].pt[1]; double w = GP[i].pt[2];
const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
std::vector<typename TensorialTraits<T1>::ValType> Vals;
LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
+ if(ele->getTypeForMSH() == MSH_LIN_B || ele->getTypeForMSH() == MSH_TRI_B ||
+ ele->getTypeForMSH() == MSH_POLYG_B)
+ ele->movePointFromParentSpaceToElementSpace(u, v, w);
SPoint3 p;
ele->pnt(u, v, w, p);
- typename TensorialTraits<T1>::ValType load = Load(p.x(), p.y(), p.z());
+ typename TensorialTraits<T1>::ValType load = (*Load)(p.x(), p.y(), p.z());
for(int j = 0; j < nbFF ; ++j){
m(j) += _eqfac * dot(Vals[j], load) * weight * detJ;
}
diff --git a/Solver/thermicSolver.cpp b/Solver/thermicSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a3ef7885ff69d36bc04fa13f0311821f29ba871c
--- /dev/null
+++ b/Solver/thermicSolver.cpp
@@ -0,0 +1,318 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#include <string.h>
+#include "GmshConfig.h"
+#include "thermicSolver.h"
+#include "linearSystemCSR.h"
+#include "linearSystemPETSc.h"
+#include "linearSystemGMM.h"
+#include "linearSystemFull.h"
+#include "Numeric.h"
+#include "GModel.h"
+#include "functionSpace.h"
+#include "terms.h"
+#include "solverAlgorithms.h"
+#include "quadratureRules.h"
+#include "solverField.h"
+#include "MPoint.h"
+#include "gmshLevelset.h"
+#if defined(HAVE_POST)
+#include "PView.h"
+#include "PViewData.h"
+#endif
+
+void thermicSolver::setMesh(const std::string &meshFileName)
+{
+ pModel = new GModel();
+ pModel->readMSH(meshFileName.c_str());
+ _dim = pModel->getNumRegions() ? 3 : 2;
+
+ if (LagSpace) delete LagSpace;
+ LagSpace = new ScalarLagrangeFunctionSpace(_tag);
+
+ if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace;
+ LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpaceOfElement(_tag + 1);
+}
+
+void thermicSolver::solve()
+{
+ linearSystemFull<double> *lsys = new linearSystemFull<double>;
+/*#if defined(HAVE_TAUCS)
+ linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#elif defined(HAVE_PETSC)
+ linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
+#else
+ linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
+ lsys->setNoisy(2);
+#endif*/
+ assemble(lsys);
+ lsys->systemSolve();
+ printf("-- done solving!\n");
+}
+
+void thermicSolver::cutMesh(gLevelset *ls)
+{
+ pModel = pModel->buildCutGModel(ls);
+ pModel->writeMSH("cutMesh.msh");
+}
+
+void thermicSolver::setThermicDomain(int phys, double k)
+{
+ thermicField field;
+ field._k = k;
+ field._tag = _tag;
+ field.g = new groupOfElements (_dim, phys);
+ thermicFields.push_back(field);
+}
+
+void thermicSolver::setLagrangeMultipliers(int phys, double tau, int tag, simpleFunction<double> *f)
+{
+ LagrangeMultiplierFieldT field;
+ field._tau = tau;
+ field._tag = tag;
+ field._f = f;
+ field.g = new groupOfElements (_dim - 1, phys);
+ LagrangeMultiplierFields.push_back(field);
+}
+
+void thermicSolver::setEdgeTemp(int edge, simpleFunction<double> *f)
+{
+ dirichletBCT diri;
+ diri.g = new groupOfElements (1, edge);
+ diri._f = f;
+ diri._tag = edge;
+ diri.onWhat = BoundaryConditionT::ON_EDGE;
+ allDirichlet.push_back(diri);
+}
+
+void thermicSolver::assemble(linearSystem<double> *lsys)
+{
+ if (pAssembler) delete pAssembler;
+ pAssembler = new dofManager<double>(lsys);
+
+ // 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
+
+ // Dirichlet conditions
+ for (unsigned int i = 0; i < allDirichlet.size(); i++){
+ FilterDofTrivial filter;
+ FixNodalDofs(*LagSpace, allDirichlet[i].g->begin(), allDirichlet[i].g->end(),
+ *pAssembler, *allDirichlet[i]._f, filter);
+ }
+ // LagrangeMultipliers
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i){
+ NumberDofs(*LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(),
+ LagrangeMultiplierFields[i].g->end(), *pAssembler);
+ }
+ // Thermic Fields
+ for (unsigned int i = 0; i < thermicFields.size(); ++i){
+ NumberDofs(*LagSpace, thermicFields[i].g->begin(), thermicFields[i].g->end(), *pAssembler);
+ }
+ // Neumann conditions
+ GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
+ for (unsigned int i = 0; i < allNeumann.size(); i++){
+ std::cout << "Neumann BC" << std::endl;
+ LoadTerm<double> Lterm(*LagSpace, allNeumann[i]._f);
+ Assemble(Lterm, *LagSpace, allNeumann[i].g->begin(), allNeumann[i].g->end(),
+ Integ_Boundary, *pAssembler);
+ }
+ // Assemble cross term, laplace term and rhs term for LM
+ GaussQuadrature Integ_LagrangeMult(GaussQuadrature::ValVal);
+ GaussQuadrature Integ_Laplace(GaussQuadrature::GradGrad);
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); i++){
+ printf("Lagrange Mult Lag\n");
+ LagrangeMultiplierTerm<double> LagTerm(*LagSpace, *LagrangeMultiplierSpace, 1.);
+ Assemble(LagTerm, *LagSpace, *LagrangeMultiplierSpace,
+ LagrangeMultiplierFields[i].g->begin(), LagrangeMultiplierFields[i].g->end(),
+ Integ_LagrangeMult, *pAssembler);
+ printf("Lagrange Mult Lap\n");
+ LaplaceTerm<double, double> LapTerm(*LagrangeMultiplierSpace,
+ -LagrangeMultiplierFields[i]._tau);
+ Assemble(LapTerm, *LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(),
+ LagrangeMultiplierFields[i].g->end(), Integ_Laplace, *pAssembler);
+ printf("Lagrange Mult Load\n");
+ LoadTermOnBorder<double> Lterm(*LagrangeMultiplierSpace, LagrangeMultiplierFields[i]._f);
+ Assemble(Lterm, *LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(),
+ LagrangeMultiplierFields[i].g->end(), Integ_Boundary, *pAssembler);
+ }
+ // Assemble thermic term
+ GaussQuadrature Integ_Bulk(GaussQuadrature::ValVal);
+ for (unsigned int i = 0; i < thermicFields.size(); i++){
+ printf("Thermic Term\n");
+ LaplaceTerm<double, double> Tterm(*LagSpace, thermicFields[i]._k);
+ Assemble(Tterm, *LagSpace, thermicFields[i].g->begin(), thermicFields[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("%g ", d);
+ }
+ double d; lsys->getFromRightHandSide(i, d);
+ printf(" | %g\n", d);
+ }*/
+
+ printf("nDofs=%d\n", pAssembler->sizeOfR());
+ printf("nFixed=%d\n", pAssembler->sizeOfF());
+}
+
+double thermicSolver::computeL2Norm(simpleFunction<double> *sol) {
+ double val = 0.0;
+ SolverField<double> solField(pAssembler, LagSpace);
+ for (unsigned int i = 0; i < thermicFields.size(); ++i){
+ for (groupOfElements::elementContainer::const_iterator it = thermicFields[i].g->begin(); it != thermicFields[i].g->end(); ++it){
+ MElement *e = *it;
+//printf("element (%g,%g) (%g,%g) (%g,%g)\n",e->getVertex(0)->x(),e->getVertex(0)->y(),e->getVertex(1)->x(),e->getVertex(1)->y(),e->getVertex(2)->x(),e->getVertex(2)->y());
+ int npts;
+ IntPt *GP;
+ double jac[3][3];
+ int integrationOrder = 2 * (e->getPolynomialOrder()+5);
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ for (int j = 0; j < npts; j++){
+ double u = GP[j].pt[0];
+ double v = GP[j].pt[1];
+ double w = GP[j].pt[2];
+ double weight = GP[j].weight;
+ double detJ = fabs(e->getJacobian (u, v, w, jac));
+ SPoint3 p;
+ e->pnt(u, v, w, p);
+ double FEMVALUE;
+ solField.f(e, u, v, w, FEMVALUE);
+ double diff = (*sol)(p.x(), p.y(), p.z()) - FEMVALUE;
+ val += diff * diff * detJ * weight;
+//printf("(%g %g) : detJ=%g we=%g FV=%g sol=%g diff=%g\n",p.x(),p.y(),detJ,weight,FEMVALUE,(*sol)(p.x(), p.y(), p.z()),diff);
+ }
+ }
+ }
+printf("L2Norm = %g\n",sqrt(val));
+ return sqrt(val);
+}
+
+double thermicSolver::computeLagNorm(int tag, simpleFunction<double> *sol) {
+ double val = 0.0, val2 = 0.0;
+ SolverField<double> solField(pAssembler, LagrangeMultiplierSpace);
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i){
+ if(tag != LagrangeMultiplierFields[i]._tag) continue;
+ for(groupOfElements::elementContainer::const_iterator it = LagrangeMultiplierFields[i].g->begin();
+ it != LagrangeMultiplierFields[i].g->end(); ++it){
+ MElement *e = *it;
+printf("element (%g,%g) (%g,%g)\n",e->getVertex(0)->x(),e->getVertex(0)->y(),e->getVertex(1)->x(),e->getVertex(1)->y());
+ int npts;
+ IntPt *GP;
+ double jac[3][3];
+ int integrationOrder = 2 * (e->getPolynomialOrder()+1);
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ for (int j = 0; j < npts; j++){
+ double u = GP[j].pt[0];
+ double v = GP[j].pt[1];
+ double w = GP[j].pt[2];
+ double weight = GP[j].weight;
+ double detJ = fabs(e->getJacobian (u, v, w, jac));
+ SPoint3 p;
+ e->getParent()->pnt(u, v, w, p);
+ double FEMVALUE;
+ solField.f(e, u, v, w, FEMVALUE);
+ double diff = (*sol)(p.x(), p.y(), p.z()) - FEMVALUE;
+ val += diff * diff * detJ * weight;
+ val2 += (*sol)(p.x(), p.y(), p.z()) * (*sol)(p.x(), p.y(), p.z()) * detJ * weight;
+printf("(%g %g) : u,v=(%g,%g) detJ=%g we=%g FV=%g sol=%g diff=%g\n",p.x(),p.y(),u,v,detJ,weight,FEMVALUE,(*sol)(p.x(), p.y(), p.z()),diff);
+ }
+ }
+ }
+printf("LagNorm = %g\n",sqrt(val/val2));
+ return sqrt(val/val2);
+}
+
+#if defined(HAVE_POST)
+
+PView* thermicSolver::buildTemperatureView (const std::string postFileName)
+{
+ std::cout << "build Temperature View"<< std::endl;
+ std::set<MVertex*> v;
+ std::map<MVertex*, MElement*> vCut;
+ for (unsigned int i = 0; i < thermicFields.size(); ++i){
+ for (groupOfElements::elementContainer::const_iterator it = thermicFields[i].g->begin(); it != thermicFields[i].g->end(); ++it){
+ MElement *e = *it;
+ if(e->getParent()) {
+ for (int j = 0; j < e->getNumVertices(); ++j) {
+ if(vCut.find(e->getVertex(j)) == vCut.end())
+ vCut[e->getVertex(j)] = e->getParent();
+ }
+ }
+ else {
+ for (int j = 0; j < e->getNumVertices(); ++j)
+ v.insert(e->getVertex(j));
+ }
+ }
+ }
+ std::map<int, std::vector<double> > data;
+ SolverField<double> Field(pAssembler, LagSpace);
+ for (std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it){
+ double val;
+ MPoint p(*it);
+ Field.f(&p, 0, 0, 0, val); //printf("valv=%g\n",val);
+ std::vector<double> vec; vec.push_back(val);
+ data[(*it)->getNum()] = vec;
+ }
+ for (std::map<MVertex*, MElement*>::iterator it = vCut.begin(); it != vCut.end(); ++it){
+ double val;
+ double uvw[3];
+ double xyz[3] = {it->first->x(), it->first->y(), it->first->z()};
+ it->second->xyz2uvw(xyz, uvw);
+ Field.f(it->second, uvw[0], uvw[1], uvw[2], val); //printf("valvc=%g\n",val);
+ std::vector<double> vec; vec.push_back(val);
+ data[it->first->getNum()] = vec;
+ }
+ PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0, 1);
+ return pv;
+}
+
+PView* thermicSolver::buildLagrangeMultiplierView (const std::string postFileName)
+{
+ std::cout << "build Lagrange Multiplier View"<< std::endl;
+ if(!LagrangeMultiplierSpace) return new PView();
+ std::set<MVertex*> v;
+ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i){
+ for(groupOfElements::elementContainer::const_iterator it = LagrangeMultiplierFields[i].g->begin();
+ it != LagrangeMultiplierFields[i].g->end(); ++it){
+ MElement *e = *it;
+ for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
+ }
+ }
+ std::map<int, std::vector<double> > data;
+ SolverField<double> Field(pAssembler, LagrangeMultiplierSpace);
+ for(std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it){
+ double val;
+ MPoint p(*it);
+ Field.f(&p, 0, 0, 0, val);
+ std::vector<double> vec;
+ vec.push_back(val);
+ data[(*it)->getNum()] = vec;
+ }
+ PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0, 1);
+ return pv;
+}
+
+#else
+PView* thermicSolver::buildTemperatureView(const std::string postFileName)
+{
+ Msg::Error("Post-pro module not available");
+ return 0;
+}
+
+PView* thermicSolver::buildLagrangeMultiplierView (const std::string postFileName)
+{
+ Msg::Error("Post-pro module not available");
+ return 0;
+}
+#endif
+
+
+
+
+
diff --git a/Solver/thermicSolver.h b/Solver/thermicSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..59c42456f0b5432edeef422b923024771c56d055
--- /dev/null
+++ b/Solver/thermicSolver.h
@@ -0,0 +1,100 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _THERMIC_SOLVER_H_
+#define _THERMIC_SOLVER_H_
+
+#include <map>
+#include <string>
+#include "SVector3.h"
+#include "dofManager.h"
+#include "simpleFunction.h"
+#include "functionSpace.h"
+
+class GModel;
+class PView;
+class groupOfElements;
+class gLevelset;
+
+struct LagrangeMultiplierFieldT {
+ int _tag;
+ groupOfElements *g;
+ double _tau;
+ simpleFunction<double> *_f;
+ LagrangeMultiplierFieldT() : _tag(0), g(0){}
+};
+
+struct thermicField {
+ int _tag; // tag for the dofManager
+ groupOfElements *g; // support for this field
+ double _k; // diffusivity
+ thermicField () : _tag(0), g(0){}
+};
+
+struct BoundaryConditionT
+{
+ 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
+ BoundaryConditionT() : _tag(0), onWhat(UNDEF), g(0) {}
+};
+
+struct dirichletBCT : public BoundaryConditionT
+{
+ simpleFunction<double> *_f;
+ dirichletBCT ():BoundaryConditionT(), _f(0){}
+};
+
+struct neumannBCT : public BoundaryConditionT
+{
+ simpleFunction<double> *_f;
+ neumannBCT () : BoundaryConditionT(), _f(0){}
+};
+// a thermic solver ...
+class thermicSolver
+{
+ protected:
+ GModel *pModel;
+ int _dim, _tag;
+ dofManager<double> *pAssembler;
+ FunctionSpace<double> *LagSpace;
+ FunctionSpace<double> *LagrangeMultiplierSpace;
+
+ // young modulus and poisson coefficient per physical
+ std::vector<thermicField> thermicFields;
+ std::vector<LagrangeMultiplierFieldT> LagrangeMultiplierFields;
+ // neumann BC
+ std::vector<neumannBCT> allNeumann;
+ // dirichlet BC
+ std::vector<dirichletBCT> allDirichlet;
+
+ public:
+ thermicSolver(int tag) : _tag(tag), pAssembler(0), LagSpace(0), LagrangeMultiplierSpace(0) {}
+
+ virtual ~thermicSolver()
+ {
+ if (LagSpace) delete LagSpace;
+ if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace;
+ if (pAssembler) delete pAssembler;
+ }
+ void assemble (linearSystem<double> *lsys);
+ virtual void setMesh(const std::string &meshFileName);
+ void cutMesh(gLevelset *ls);
+ void setThermicDomain(int phys, double k);
+ void setLagrangeMultipliers(int phys, double tau, int tag, simpleFunction<double> *f);
+ void setEdgeTemp(int edge, simpleFunction<double> *f);
+ void solve();
+ virtual PView *buildTemperatureView(const std::string postFileName);
+ virtual PView *buildLagrangeMultiplierView(const std::string posFileName);
+ double computeL2Norm(simpleFunction<double> *f);
+ double computeLagNorm(int tag, simpleFunction<double> *f);
+ // std::pair<PView *, PView*> buildErrorEstimateView
+ // (const std::string &errorFileName, double, int);
+ // std::pair<PView *, PView*> buildErrorEstimateView
+ // (const std::string &errorFileName, const elasticityData &ref, double, int);
+};
+
+#endif
diff --git a/wrappers/gmshpy/gmshPost.i b/wrappers/gmshpy/gmshPost.i
index d30cc19dba658416f19e463a0fc9416965cd6aae..6b68a2caa2915b865c5fe2c6d6109a5afdd07c7a 100644
--- a/wrappers/gmshpy/gmshPost.i
+++ b/wrappers/gmshpy/gmshPost.i
@@ -16,6 +16,7 @@
#include "PViewFactory.h"
#include "PViewData.h"
#include "PViewAsSimpleFunction.h"
+ #include "PViewDataGModel.h"
#endif
%}
@@ -34,6 +35,7 @@ namespace std {
%include "simpleFunction.h"
%template(simpleFunctionDouble) simpleFunction<double>;
%include "PViewAsSimpleFunction.h"
+%include "PViewDataGModel.h"
%include "Plugin.h"
%include "PluginManager.h"
#endif
diff --git a/wrappers/gmshpy/gmshSolver.i b/wrappers/gmshpy/gmshSolver.i
index 881d8ad93859331e7ff5367fec2375c36ad83a09..6c22c518dd5088610ba635137cdb65074af16b63 100644
--- a/wrappers/gmshpy/gmshSolver.i
+++ b/wrappers/gmshpy/gmshSolver.i
@@ -10,6 +10,7 @@
#if defined(HAVE_SOLVER)
#include "dofManager.h"
#include "elasticitySolver.h"
+ #include "thermicSolver.h"
#include "frameSolver.h"
#include "linearSystem.h"
#include "linearSystemCSR.h"
@@ -24,6 +25,7 @@
%include "dofManager.h"
%template(dofManagerDouble) dofManager<double>;
%include "elasticitySolver.h"
+%include "thermicSolver.h"
%include "frameSolver.h"
%include "linearSystem.h"
%template(linearSystemDouble) linearSystem<double>;