diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index be544939fc1b15d55efc410c1b388cd398127fe9..ceb4c2d9df1645b0b7db2543f0027e7b502c3151 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -25,31 +25,29 @@
#include "PViewData.h"
#endif
+
/*
-static void printLinearSystem(linearSystemCSRTaucs<double> *lsys)
-{
- int *startIndex;
- int *columns;
- double *values;
- lsys->getMatrix(startIndex, columns, values);
- for(int i = 0; i < lsys->getNbUnk(); i++){
- std::cout<<"a ";
- for(int j = 0; j < lsys->getNbUnk(); j++){
- double val = 0.;
- for(int k = startIndex[i]; k < startIndex[i+1]; k++){
- if(columns[k] == j) {val = values[k]; break;}
- }
- std::cout<< val <<" ";
+static void printLinearSystem(linearSystem<double> *lsys, int size)
+{
+ if(!lsys->isAllocated()){printf("Linear system not allocated\n"); return;}
+ double valeur;
+ printf("K\n");
+ for (int i = 0; i < size; i++){
+ for (int j = 0; j < size; j++ ){
+ lsys->getFromMatrix (i, j, valeur);
+ printf ("%g ", valeur) ;
}
- std::cout<<std::endl;
- }std::cout<<std::endl;
- for(int i = 0; i < lsys->getNbUnk(); i++) {
- double val; lsys->getFromRightHandSide(i,val);
- std::cout<<"b "<<val<<std::endl;
- }std::cout<<std::endl;
- for(int i = 0; i < lsys->getNbUnk(); i++) {
- double val; lsys->getFromSolution(i,val);
- std::cout<<"x "<<val<<std::endl;
+ printf("\n");
+ }
+ printf("b\n");
+ for (int i = 0; i < size; i++){
+ lsys->getFromRightHandSide(i, valeur);
+ printf("%g\n", valeur);
+ }
+ printf("x\n");
+ for (int i = 0; i < size; i++){
+ lsys->getFromSolution(i, valeur);
+ printf("%g\n", valeur);
}
}
*/
@@ -67,7 +65,7 @@ void elasticitySolver::setMesh(const std::string &meshFileName)
VectorLagrangeFunctionSpace::VECTOR_Y);
if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace;
- LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpace(_tag+1);
+ LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpaceOfElement(_tag+1);
}
void elasticitySolver::exportKb()
@@ -97,7 +95,6 @@ void elasticitySolver::exportKb()
void elasticitySolver::solve()
{
-
//linearSystemFull<double> *lsys = new linearSystemFull<double>;
#if defined(HAVE_TAUCS)
@@ -110,13 +107,12 @@ void elasticitySolver::solve()
#endif
assemble(lsys);
+ //printLinearSystem(lsys,pAssembler->sizeOfR());
lsys->systemSolve();
printf("-- done solving!\n");
- GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
-
-// printLinearSystem(lsys);
double energ=0;
+ GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); i++){
SolverField<SVector3> Field(pAssembler, LagSpace);
IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
@@ -157,7 +153,8 @@ void elasticitySolver::readInputFile(const std::string &fn)
}
if(what[0]=='#'){
char buffer[1024];
- fgets(buffer, sizeof(buffer), f);
+ if(fgets(buffer, sizeof(buffer), f) == NULL)
+ Msg::Error("Cannot read line.");
}
else if (!strcmp(what, "ElasticDomain")){
elasticField field;
@@ -348,6 +345,46 @@ void elasticitySolver::readInputFile(const std::string &fn)
}
+void elasticitySolver::cutMesh(gLevelset *ls)
+{
+ pModel = pModel->buildCutGModel(ls);
+ pModel->writeMSH("cutMesh.msh");
+}
+
+void elasticitySolver::setElasticDomain(int phys, double E, double nu)
+{
+ elasticField field;
+ field._E = E;
+ field._nu = nu;
+ field._tag = _tag;
+ field.g = new groupOfElements (_dim, phys);
+ elasticFields.push_back(field);
+}
+
+void elasticitySolver::setLagrangeMultipliers(int phys, double tau, SVector3 d,
+ int tag, simpleFunction<double> *f)
+{
+ LagrangeMultiplierField field;
+ field._tau = tau;
+ field._tag = tag;
+ field._f = f;
+ field._d = d;
+ field.g = new groupOfElements (_dim - 1, phys);
+ LagrangeMultiplierFields.push_back(field);
+}
+
+void elasticitySolver::setEdgeDisp(int edge, int comp, simpleFunction<double> *f)
+{
+ dirichletBC diri;
+ diri.g = new groupOfElements (1, edge);
+ diri._f = f;
+ diri._comp = comp;
+ diri._tag = edge;
+ diri.onWhat = BoundaryCondition::ON_EDGE;
+ allDirichlet.push_back(diri);
+}
+
+
void elasticitySolver::addDirichletBC(int dim, int entityId, int component, double value)
{
dirichletBC diri;
@@ -418,7 +455,7 @@ elasticitySolver::elasticitySolver(GModel *model, int tag)
if (_dim==2) LagSpace=new VectorLagrangeFunctionSpace(_tag,
VectorLagrangeFunctionSpace::VECTOR_X,
VectorLagrangeFunctionSpace::VECTOR_Y);
- LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpace(_tag+1);
+ LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpaceOfElement(_tag+1);
}
void elasticitySolver::assemble(linearSystem<double> *lsys)
@@ -465,18 +502,21 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
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<SVector3> LagTerm(*LagSpace, *LagrangeMultiplierSpace,
LagrangeMultiplierFields[i]._d);
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);
- LoadTerm<double> Lterm(*LagrangeMultiplierSpace, LagrangeMultiplierFields[i]._f);
+ 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);
}
@@ -488,15 +528,6 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
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());
@@ -606,6 +637,53 @@ void elasticitySolver::computeEffectiveStrain(std::vector<double> strain)
strain[i] = st[i] / volTot;
}
+double elasticitySolver::computeDisplacementError(int comp, simpleFunction<double> *sol) {
+ std::cout << "compute displacement error"<< std::endl;
+ double err = 0.;
+ std::set<MVertex*> v;
+ 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;
+ 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));
+ }
+ }
+ }
+ SolverField<SVector3> Field(pAssembler, LagSpace);
+ for (std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it){
+ SVector3 val;
+ MPoint p(*it);
+ Field.f(&p, 0, 0, 0, val);
+ double diff = (*sol)((*it)->x(), (*it)->y(), (*it)->z()) - val(comp);
+ err += diff * diff;
+ }
+ 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);
+ double diff = (*sol)(xyz[0], xyz[1], xyz[2]) - val(comp);
+ err += diff * diff;
+ }
+printf("Displacement Error (comp=%d) = %g\n",comp,sqrt(err));
+ return sqrt(err);
+}
+
+double elasticitySolver::computeLagNorm(int tag, simpleFunction<double> *sol) {
+ return 0;
+}
+
#if defined(HAVE_POST)
PView* elasticitySolver::buildDisplacementView (const std::string postFileName)
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 16d69d3aaa5727987431bb63717544ff54ec12c0..7419c637f1456f32b57d85b1e2ec841518bf7d11 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -16,6 +16,7 @@ template <class scalar> class simpleFunction;
class GModel;
class PView;
class groupOfElements;
+class gLevelset;
struct LagrangeMultiplierField {
int _tag;
@@ -95,6 +96,10 @@ class elasticitySolver
void readInputFile(const std::string &meshFileName);
void read(const std::string s) {readInputFile(s.c_str());}
virtual void setMesh(const std::string &meshFileName);
+ void cutMesh(gLevelset *ls);
+ void setElasticDomain(int phys, double E, double nu);
+ void setLagrangeMultipliers(int phys, double tau, SVector3 d, int tag, simpleFunction<double> *f);
+ void setEdgeDisp(int edge, int comp, simpleFunction<double> *f);
void solve();
void postSolve();
void exportKb();
@@ -107,6 +112,8 @@ class elasticitySolver
virtual PView *buildElasticEnergyView(const std::string postFileName);
virtual PView *buildVonMisesView(const std::string postFileName);
virtual PView *buildVolumeView(const std::string postFileName);
+ double computeDisplacementError(int comp, 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
diff --git a/Solver/thermicSolver.cpp b/Solver/thermicSolver.cpp
index d4240179187fba0cfa1c390522910ef7154cb06a..d2b8672cec6d8412de8eb6863d764bbc8542df99 100644
--- a/Solver/thermicSolver.cpp
+++ b/Solver/thermicSolver.cpp
@@ -68,6 +68,15 @@ void thermicSolver::setThermicDomain(int phys, double k)
thermicFields.push_back(field);
}
+void thermicSolver::changeLMTau(int tag, double tau)
+{
+ for(unsigned int i = 0; i < LagrangeMultiplierFields.size(); i++){
+ if(LagrangeMultiplierFields[i]._tag == tag){
+ LagrangeMultiplierFields[i]._tau = tau;
+ }
+ }
+}
+
void thermicSolver::setLagrangeMultipliers(int phys, double tau, int tag, simpleFunction<double> *f)
{
LagrangeMultiplierFieldT field;
@@ -88,6 +97,16 @@ void thermicSolver::setEdgeTemp(int edge, simpleFunction<double> *f)
allDirichlet.push_back(diri);
}
+void thermicSolver::setFaceTemp(int face, simpleFunction<double> *f)
+{
+ dirichletBCT diri;
+ diri.g = new groupOfElements (2, face);
+ diri._f = f;
+ diri._tag = face;
+ diri.onWhat = BoundaryConditionT::ON_FACE;
+ allDirichlet.push_back(diri);
+}
+
void thermicSolver::assemble(linearSystem<double> *lsys)
{
if (pAssembler) delete pAssembler;
@@ -201,7 +220,7 @@ double thermicSolver::computeLagNorm(int tag, simpleFunction<double> *sol) {
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());
+//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];
@@ -220,7 +239,7 @@ printf("element (%g,%g) (%g,%g)\n",e->getVertex(0)->x(),e->getVertex(0)->y(),e->
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("(%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);
}
}
}
@@ -298,6 +317,45 @@ PView* thermicSolver::buildLagrangeMultiplierView (const std::string postFileNam
return pv;
}
+PView* thermicSolver::buildErrorEstimateView(const std::string errorFileName,
+ simpleFunction<double> *sol)
+{
+ std::cout << "build Error View"<< std::endl;
+ std::map<int, std::vector<double> > data;
+
+ 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;
+ int npts;
+ IntPt *GP;
+ double jac[3][3];
+ int integrationOrder = 2 * (e->getPolynomialOrder()+5);
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ double val = 0.0;
+ 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;
+ }
+ std::vector<double> vec;
+ vec.push_back(sqrt(val));
+ data[e->getNum()] = vec;
+ }
+ }
+
+ PView *pv = new PView (errorFileName, "ElementData", pModel, data, 0.0, 1);
+ return pv;
+}
+
#else
PView* thermicSolver::buildTemperatureView(const std::string postFileName)
{
@@ -310,6 +368,12 @@ PView* thermicSolver::buildLagrangeMultiplierView (const std::string postFileNam
Msg::Error("Post-pro module not available");
return 0;
}
+
+PView* thermicSolver::buildErrorEstimateView(const std::string errorFileName)
+{
+ Msg::Error("Post-pro module not available");
+ return 0;
+}
#endif
diff --git a/Solver/thermicSolver.h b/Solver/thermicSolver.h
index 2f715c9149f085b3e520fa33e983ce82a5fceed3..d884c6f5140334f01db0a316c583eb1070403776 100644
--- a/Solver/thermicSolver.h
+++ b/Solver/thermicSolver.h
@@ -85,14 +85,16 @@ class thermicSolver
void cutMesh(gLevelset *ls);
void setThermicDomain(int phys, double k);
void setLagrangeMultipliers(int phys, double tau, int tag, simpleFunction<double> *f);
+ void changeLMTau(int tag, double tau);
void setEdgeTemp(int edge, simpleFunction<double> *f);
+ void setFaceTemp(int face, simpleFunction<double> *f);
void solve();
virtual PView *buildTemperatureView(const std::string postFileName);
- virtual PView *buildLagrangeMultiplierView(const std::string posFileName);
+ virtual PView *buildLagrangeMultiplierView(const std::string postFileName);
double computeL2Norm(simpleFunction<double> *f);
double computeLagNorm(int tag, simpleFunction<double> *f);
- // std::pair<PView *, PView*> buildErrorEstimateView
- // (const std::string &errorFileName, double, int);
+ PView* buildErrorEstimateView(const std::string errorFileName,
+ simpleFunction<double> *sol);
// std::pair<PView *, PView*> buildErrorEstimateView
// (const std::string &errorFileName, const elasticityData &ref, double, int);
};