Skip to content
Snippets Groups Projects
Commit a2729cc1 authored by Gaetan Bricteux's avatar Gaetan Bricteux
Browse files

enable multiple LM fields + postpro

parent 86c615dd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Solver/elasticitySolver.cpp
+ 155
− 41
View file @ a2729cc1
...@@ -25,7 +25,6 @@ ...@@ -25,7 +25,6 @@
#include "PViewData.h" #include "PViewData.h"
#endif #endif
/* /*
static void printLinearSystem(linearSystem<double> *lsys, int size) static void printLinearSystem(linearSystem<double> *lsys, int size)
{ {
...@@ -52,25 +51,41 @@ static void printLinearSystem(linearSystem<double> *lsys, int size) ...@@ -52,25 +51,41 @@ static void printLinearSystem(linearSystem<double> *lsys, int size)
} }
*/ */
void elasticitySolver::setMesh(const std::string &meshFileName) elasticitySolver::elasticitySolver(GModel *model, int tag)
{
pModel = model;
_dim = pModel->getNumRegions() ? 3 : 2;
_tag = tag;
pAssembler = NULL;
if (_dim==3) LagSpace=new VectorLagrangeFunctionSpace(_tag);
if (_dim==2) LagSpace=new VectorLagrangeFunctionSpace(_tag,
VectorLagrangeFunctionSpace::VECTOR_X,
VectorLagrangeFunctionSpace::VECTOR_Y);
}
void elasticitySolver::setMesh(const std::string &meshFileName, int dim)
{ {
pModel = new GModel(); pModel = new GModel();
pModel->readMSH(meshFileName.c_str()); pModel->readMSH(meshFileName.c_str());
_dim = pModel->getNumRegions() ? 3 : 2; _dim = pModel->getNumRegions() ? 3 : 2;
if (LagSpace) delete LagSpace; if (LagSpace) delete LagSpace;
if (_dim==3) LagSpace = new VectorLagrangeFunctionSpace(_tag); if (dim == 3 || _dim==3) LagSpace = new VectorLagrangeFunctionSpace(_tag);
if (_dim==2) LagSpace = new VectorLagrangeFunctionSpace else if (dim == 2 || _dim==2) LagSpace = new VectorLagrangeFunctionSpace
(_tag,VectorLagrangeFunctionSpace::VECTOR_X, (_tag,VectorLagrangeFunctionSpace::VECTOR_X,
VectorLagrangeFunctionSpace::VECTOR_Y); VectorLagrangeFunctionSpace::VECTOR_Y);
if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace; for (unsigned int i = 0; i < LagrangeMultiplierSpaces.size(); i++)
LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpaceOfElement(_tag+1); if (LagrangeMultiplierSpaces[i]) delete LagrangeMultiplierSpaces[i];
LagrangeMultiplierSpaces.clear();
} }
void elasticitySolver::exportKb() void elasticitySolver::exportKb()
{ {
std::string sysname = "A"; std::string sysname = "A";
if(!pAssembler ||
!pAssembler->getLinearSystem(sysname) ||
!pAssembler->getLinearSystem(sysname)->isAllocated()) return;
double valeur; double valeur;
FILE *f = Fopen("K.txt", "w"); FILE *f = Fopen("K.txt", "w");
if(f){ if(f){
...@@ -95,8 +110,11 @@ void elasticitySolver::exportKb() ...@@ -95,8 +110,11 @@ void elasticitySolver::exportKb()
void elasticitySolver::solve() void elasticitySolver::solve()
{ {
//linearSystemFull<double> *lsys = new linearSystemFull<double>; std::string sysname = "A";
if(pAssembler && pAssembler->getLinearSystem(sysname))
delete pAssembler->getLinearSystem(sysname);
//linearSystemFull<double> *lsys = new linearSystemFull<double>;
#if defined(HAVE_TAUCS) #if defined(HAVE_TAUCS)
linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>; linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#elif defined(HAVE_PETSC) #elif defined(HAVE_PETSC)
...@@ -176,11 +194,12 @@ void elasticitySolver::readInputFile(const std::string &fn) ...@@ -176,11 +194,12 @@ void elasticitySolver::readInputFile(const std::string &fn)
fclose(f); fclose(f);
return; return;
} }
field._tag = _tag; SVector3 sv(d1, d2, d3);
field._d = SVector3(d1, d2, d3); field._d = sv.unit();
field._f = new simpleFunction<double>(val); field._f = new simpleFunction<double>(val);
field.g = new groupOfElements (_dim - 1, physical); field.g = new groupOfElements (_dim - 1, physical);
LagrangeMultiplierFields.push_back(field); LagrangeMultiplierFields.push_back(field);
LagrangeMultiplierSpaces.push_back(new ScalarLagrangeFunctionSpaceOfElement(field._tag));
} }
else if (!strcmp(what, "Void")){ else if (!strcmp(what, "Void")){
elasticField field; elasticField field;
...@@ -368,9 +387,20 @@ void elasticitySolver::setLagrangeMultipliers(int phys, double tau, SVector3 d, ...@@ -368,9 +387,20 @@ void elasticitySolver::setLagrangeMultipliers(int phys, double tau, SVector3 d,
field._tau = tau; field._tau = tau;
field._tag = tag; field._tag = tag;
field._f = f; field._f = f;
field._d = d; field._d = d.unit();
field.g = new groupOfElements (_dim - 1, phys); field.g = new groupOfElements (_dim - 1, phys); // LM applied on entity of dimension = dim-1!
LagrangeMultiplierFields.push_back(field); LagrangeMultiplierFields.push_back(field);
LagrangeMultiplierSpaces.push_back(new ScalarLagrangeFunctionSpaceOfElement(tag));
}
void elasticitySolver::changeLMTau(int tag, double tau)
{
for(unsigned int i = 0; i < LagrangeMultiplierFields.size(); i++){
if(LagrangeMultiplierFields[i]._tag == tag){
LagrangeMultiplierFields[i]._tau = tau;
}
}
} }
void elasticitySolver::setEdgeDisp(int edge, int comp, simpleFunction<double> *f) void elasticitySolver::setEdgeDisp(int edge, int comp, simpleFunction<double> *f)
...@@ -445,19 +475,6 @@ void elasticitySolver::addElasticDomain(std::string phys, double e, double nu) ...@@ -445,19 +475,6 @@ void elasticitySolver::addElasticDomain(std::string phys, double e, double nu)
addElasticDomain(entityId, e, nu); addElasticDomain(entityId, e, nu);
} }
elasticitySolver::elasticitySolver(GModel *model, int tag)
{
pModel = model;
_dim = pModel->getNumRegions() ? 3 : 2;
_tag = tag;
pAssembler = NULL;
if (_dim==3) LagSpace=new VectorLagrangeFunctionSpace(_tag);
if (_dim==2) LagSpace=new VectorLagrangeFunctionSpace(_tag,
VectorLagrangeFunctionSpace::VECTOR_X,
VectorLagrangeFunctionSpace::VECTOR_Y);
LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpaceOfElement(_tag+1);
}
void elasticitySolver::assemble(linearSystem<double> *lsys) void elasticitySolver::assemble(linearSystem<double> *lsys)
{ {
if (pAssembler) delete pAssembler; if (pAssembler) delete pAssembler;
...@@ -475,7 +492,10 @@ void elasticitySolver::assemble(linearSystem<double> *lsys) ...@@ -475,7 +492,10 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
} }
// LagrangeMultipliers // 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(), unsigned int j = 0;
for ( ; j < LagrangeMultiplierSpaces.size(); j++)
if (LagrangeMultiplierSpaces[j]->getId() == LagrangeMultiplierFields[i]._tag) break;
NumberDofs(*(LagrangeMultiplierSpaces[j]), LagrangeMultiplierFields[i].g->begin(),
LagrangeMultiplierFields[i].g->end(), *pAssembler); LagrangeMultiplierFields[i].g->end(), *pAssembler);
} }
// Elastic Fields // Elastic Fields
...@@ -502,27 +522,31 @@ void elasticitySolver::assemble(linearSystem<double> *lsys) ...@@ -502,27 +522,31 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
GaussQuadrature Integ_LagrangeMult(GaussQuadrature::ValVal); GaussQuadrature Integ_LagrangeMult(GaussQuadrature::ValVal);
GaussQuadrature Integ_Laplace(GaussQuadrature::GradGrad); GaussQuadrature Integ_Laplace(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); i++){ for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); i++){
unsigned int j = 0;
for ( ; j < LagrangeMultiplierSpaces.size(); j++)
if (LagrangeMultiplierSpaces[j]->getId() == LagrangeMultiplierFields[i]._tag) break;
printf("Lagrange Mult Lag\n"); printf("Lagrange Mult Lag\n");
LagrangeMultiplierTerm<SVector3> LagTerm(*LagSpace, *LagrangeMultiplierSpace, LagrangeMultiplierTerm<SVector3> LagTerm(*LagSpace, *(LagrangeMultiplierSpaces[j]),
LagrangeMultiplierFields[i]._d); LagrangeMultiplierFields[i]._d);
Assemble(LagTerm, *LagSpace, *LagrangeMultiplierSpace, Assemble(LagTerm, *LagSpace, *(LagrangeMultiplierSpaces[j]),
LagrangeMultiplierFields[i].g->begin(), LagrangeMultiplierFields[i].g->begin(),
LagrangeMultiplierFields[i].g->end(), Integ_LagrangeMult, *pAssembler); LagrangeMultiplierFields[i].g->end(), Integ_LagrangeMult, *pAssembler);
printf("Lagrange Mult Lap\n"); printf("Lagrange Mult Lap\n");
LaplaceTerm<double,double> LapTerm(*LagrangeMultiplierSpace, LaplaceTerm<double,double> LapTerm(*(LagrangeMultiplierSpaces[j]),
LagrangeMultiplierFields[i]._tau); -LagrangeMultiplierFields[i]._tau);
Assemble(LapTerm, *LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(), Assemble(LapTerm, *(LagrangeMultiplierSpaces[j]), LagrangeMultiplierFields[i].g->begin(),
LagrangeMultiplierFields[i].g->end(), Integ_Laplace, *pAssembler); LagrangeMultiplierFields[i].g->end(), Integ_Laplace, *pAssembler);
printf("Lagrange Mult Load\n"); printf("Lagrange Mult Load\n");
LoadTermOnBorder<double> Lterm(*LagrangeMultiplierSpace, LagrangeMultiplierFields[i]._f); LoadTermOnBorder<double> Lterm(*(LagrangeMultiplierSpaces[j]), LagrangeMultiplierFields[i]._f);
Assemble(Lterm, *LagrangeMultiplierSpace, LagrangeMultiplierFields[i].g->begin(), Assemble(Lterm, *(LagrangeMultiplierSpaces[j]), LagrangeMultiplierFields[i].g->begin(),
LagrangeMultiplierFields[i].g->end(), Integ_Boundary, *pAssembler); LagrangeMultiplierFields[i].g->end(), Integ_Boundary, *pAssembler);
} }
// Assemble elastic term for // Assemble elastic term for
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad); GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); i++){ for (unsigned int i = 0; i < elasticFields.size(); i++){
printf("Elastic\n");
IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu); IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
Assemble(Eterm, *LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(), Assemble(Eterm, *LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),
Integ_Bulk, *pAssembler); Integ_Bulk, *pAssembler);
...@@ -637,7 +661,8 @@ void elasticitySolver::computeEffectiveStrain(std::vector<double> strain) ...@@ -637,7 +661,8 @@ void elasticitySolver::computeEffectiveStrain(std::vector<double> strain)
strain[i] = st[i] / volTot; strain[i] = st[i] / volTot;
} }
double elasticitySolver::computeDisplacementError(int comp, simpleFunction<double> *sol) { double elasticitySolver::computeDisplacementError(simpleFunction<double> *f0, simpleFunction<double> *f1,
simpleFunction<double> *f2) {
std::cout << "compute displacement error"<< std::endl; std::cout << "compute displacement error"<< std::endl;
double err = 0.; double err = 0.;
std::set<MVertex*> v; std::set<MVertex*> v;
...@@ -664,8 +689,10 @@ double elasticitySolver::computeDisplacementError(int comp, simpleFunction<doubl ...@@ -664,8 +689,10 @@ double elasticitySolver::computeDisplacementError(int comp, simpleFunction<doubl
SVector3 val; SVector3 val;
MPoint p(*it); MPoint p(*it);
Field.f(&p, 0, 0, 0, val); Field.f(&p, 0, 0, 0, val);
double diff = (*sol)((*it)->x(), (*it)->y(), (*it)->z()) - val(comp); SVector3 sol((*f0)((*it)->x(), (*it)->y(), (*it)->z()), (*f1)((*it)->x(), (*it)->y(), (*it)->z()),
err += diff * diff; (*f2)((*it)->x(), (*it)->y(), (*it)->z()));
double diff = normSq(sol - val);
err += diff;
} }
for (std::map<MVertex*,MElement*>::iterator it = vCut.begin(); it != vCut.end(); ++it){ for (std::map<MVertex*,MElement*>::iterator it = vCut.begin(); it != vCut.end(); ++it){
SVector3 val; SVector3 val;
...@@ -673,19 +700,95 @@ double elasticitySolver::computeDisplacementError(int comp, simpleFunction<doubl ...@@ -673,19 +700,95 @@ double elasticitySolver::computeDisplacementError(int comp, simpleFunction<doubl
double xyz[3] = {it->first->x(), it->first->y(), it->first->z()}; double xyz[3] = {it->first->x(), it->first->y(), it->first->z()};
it->second->xyz2uvw(xyz, uvw); it->second->xyz2uvw(xyz, uvw);
Field.f(it->second, uvw[0], uvw[1], uvw[2], val); Field.f(it->second, uvw[0], uvw[1], uvw[2], val);
double diff = (*sol)(xyz[0], xyz[1], xyz[2]) - val(comp); SVector3 sol((*f0)(xyz[0], xyz[1], xyz[2]), (*f1)(xyz[0], xyz[1], xyz[2]),
err += diff * diff; (*f2)(xyz[0], xyz[1], xyz[2]));
double diff = normSq(sol - val);
err += diff;
} }
printf("Displacement Error (comp=%d) = %g\n",comp,sqrt(err)); printf("Displacement Error = %g\n",sqrt(err));
return sqrt(err); return sqrt(err);
} }
double elasticitySolver::computeL2Norm(simpleFunction<double> *f0, simpleFunction<double> *f1,
simpleFunction<double> *f2) {
double val = 0.0;
SolverField<SVector3> solField(pAssembler, LagSpace);
for (unsigned int i = 0; i < elasticFields.size(); ++i){
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[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);
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);
SVector3 FEMVALUE;
solField.f(e, u, v, w, FEMVALUE);
SVector3 sol((*f0)(p.x(), p.y(), p.z()), (*f1)(p.x(), p.y(), p.z()),
(*f2)(p.x(), p.y(), p.z()));
double diff = normSq(sol - FEMVALUE);
val += diff * detJ * weight;
}
}
}
printf("L2Norm = %g\n",sqrt(val));
return sqrt(val);
}
double elasticitySolver::computeLagNorm(int tag, simpleFunction<double> *sol) { double elasticitySolver::computeLagNorm(int tag, simpleFunction<double> *sol) {
return 0; return 0;
} }
#if defined(HAVE_POST) #if defined(HAVE_POST)
PView *elasticitySolver::buildErrorView(const std::string postFileName, simpleFunction<double> *f0,
simpleFunction<double> *f1, simpleFunction<double> *f2)
{
std::cout << "build Error View"<< std::endl;
std::map<int, std::vector<double> > data;
SolverField<SVector3> solField(pAssembler, LagSpace);
for (unsigned int i = 0; i < elasticFields.size(); ++i){
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[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);
SVector3 FEMVALUE;
solField.f(e, u, v, w, FEMVALUE);
SVector3 sol((*f0)(p.x(), p.y(), p.z()), (*f1)(p.x(), p.y(), p.z()),
(*f2)(p.x(), p.y(), p.z()));
double diff = normSq(sol - FEMVALUE);
val += diff * detJ * weight;
}
std::vector<double> vec;
vec.push_back(sqrt(val));
data[e->getNum()] = vec;
}
}
PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0, 1);
return pv;
}
PView* elasticitySolver::buildDisplacementView (const std::string postFileName) PView* elasticitySolver::buildDisplacementView (const std::string postFileName)
{ {
std::cout << "build Displacement View"<< std::endl; std::cout << "build Displacement View"<< std::endl;
...@@ -862,10 +965,14 @@ PView* elasticitySolver::buildStrainView (const std::string postFileName) ...@@ -862,10 +965,14 @@ PView* elasticitySolver::buildStrainView (const std::string postFileName)
return pv; return pv;
} }
PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFileName) PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFileName, int tag)
{ {
std::cout << "build Lagrange Multiplier View"<< std::endl; std::cout << "build Lagrange Multiplier View"<< std::endl;
if(!LagrangeMultiplierSpace) return new PView(); unsigned int t = 0;
if(tag != -1)
for ( ; t < LagrangeMultiplierSpaces.size(); t++)
if (LagrangeMultiplierSpaces[t]->getId() == tag) break;
if(t == LagrangeMultiplierSpaces.size()) return new PView();
std::set<MVertex*> v; 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 = for(groupOfElements::elementContainer::const_iterator it =
...@@ -876,7 +983,7 @@ PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFile ...@@ -876,7 +983,7 @@ PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFile
} }
} }
std::map<int, std::vector<double> > data; std::map<int, std::vector<double> > data;
SolverField<double> Field(pAssembler, LagrangeMultiplierSpace); SolverField<double> Field(pAssembler, LagrangeMultiplierSpaces[t]);
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; double val;
MPoint p(*it); MPoint p(*it);
...@@ -986,6 +1093,13 @@ PView *elasticitySolver::buildVonMisesView(const std::string postFileName) ...@@ -986,6 +1093,13 @@ PView *elasticitySolver::buildVonMisesView(const std::string postFileName)
#else #else
PView* elasticitySolver::buildErrorView(int comp, simpleFunction<double> *sol,
const std::string postFileName)
{
Msg::Error("Post-pro module not available");
return 0;
}
PView* elasticitySolver::buildDisplacementView(const std::string postFileName) PView* elasticitySolver::buildDisplacementView(const std::string postFileName)
{ {
Msg::Error("Post-pro module not available"); Msg::Error("Post-pro module not available");
......
Solver/elasticitySolver.h
+ 13
− 6
View file @ a2729cc1
...@@ -63,7 +63,7 @@ class elasticitySolver ...@@ -63,7 +63,7 @@ class elasticitySolver
int _dim, _tag; int _dim, _tag;
dofManager<double> *pAssembler; dofManager<double> *pAssembler;
FunctionSpace<SVector3> *LagSpace; FunctionSpace<SVector3> *LagSpace;
FunctionSpace<double> *LagrangeMultiplierSpace; std::vector<FunctionSpace<double> *> LagrangeMultiplierSpaces;
// young modulus and poisson coefficient per physical // young modulus and poisson coefficient per physical
std::vector<elasticField> elasticFields; std::vector<elasticField> elasticFields;
...@@ -75,7 +75,7 @@ class elasticitySolver ...@@ -75,7 +75,7 @@ class elasticitySolver
std::vector<dirichletBC> allDirichlet; std::vector<dirichletBC> allDirichlet;
public: public:
elasticitySolver(int tag) : _tag(tag), pAssembler(0), LagSpace(0), LagrangeMultiplierSpace(0) {} elasticitySolver(int tag) : _tag(tag), pAssembler(0), LagSpace(0) {}
elasticitySolver(GModel *model, int tag); elasticitySolver(GModel *model, int tag);
...@@ -89,16 +89,19 @@ class elasticitySolver ...@@ -89,16 +89,19 @@ class elasticitySolver
virtual ~elasticitySolver() virtual ~elasticitySolver()
{ {
if (LagSpace) delete LagSpace; if (LagSpace) delete LagSpace;
if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace; for (unsigned int i = 0; i < LagrangeMultiplierSpaces.size(); i++)
if (LagrangeMultiplierSpaces[i]) delete LagrangeMultiplierSpaces[i];
LagrangeMultiplierSpaces.clear();
if (pAssembler) delete pAssembler; if (pAssembler) delete pAssembler;
} }
void assemble (linearSystem<double> *lsys); void assemble (linearSystem<double> *lsys);
void readInputFile(const std::string &meshFileName); void readInputFile(const std::string &meshFileName);
void read(const std::string s) {readInputFile(s.c_str());} void read(const std::string s) {readInputFile(s.c_str());}
virtual void setMesh(const std::string &meshFileName); virtual void setMesh(const std::string &meshFileName, int dim = 0);
void cutMesh(gLevelset *ls); void cutMesh(gLevelset *ls);
void setElasticDomain(int phys, double E, double nu); void setElasticDomain(int phys, double E, double nu);
void setLagrangeMultipliers(int phys, double tau, SVector3 d, int tag, simpleFunction<double> *f); void setLagrangeMultipliers(int phys, double tau, SVector3 d, int tag, simpleFunction<double> *f);
void changeLMTau(int tag, double tau);
void setEdgeDisp(int edge, int comp, simpleFunction<double> *f); void setEdgeDisp(int edge, int comp, simpleFunction<double> *f);
void solve(); void solve();
void postSolve(); void postSolve();
...@@ -108,11 +111,15 @@ class elasticitySolver ...@@ -108,11 +111,15 @@ class elasticitySolver
virtual PView *buildDisplacementView(const std::string postFileName); virtual PView *buildDisplacementView(const std::string postFileName);
virtual PView *buildStrainView(const std::string postFileName); virtual PView *buildStrainView(const std::string postFileName);
virtual PView *buildStressesView(const std::string postFileName); virtual PView *buildStressesView(const std::string postFileName);
virtual PView *buildLagrangeMultiplierView(const std::string posFileName); virtual PView *buildLagrangeMultiplierView(const std::string posFileName, int tag = -1);
virtual PView *buildElasticEnergyView(const std::string postFileName); virtual PView *buildElasticEnergyView(const std::string postFileName);
virtual PView *buildVonMisesView(const std::string postFileName); virtual PView *buildVonMisesView(const std::string postFileName);
virtual PView *buildVolumeView(const std::string postFileName); virtual PView *buildVolumeView(const std::string postFileName);
double computeDisplacementError(int comp, simpleFunction<double> *f); virtual PView *buildErrorView(const std::string postFileName, simpleFunction<double> *f0,
simpleFunction<double> *f1, simpleFunction<double> *f2);
double computeDisplacementError(simpleFunction<double> *f0, simpleFunction<double> *f1,
simpleFunction<double> *f2);
double computeL2Norm(simpleFunction<double> *f0, simpleFunction<double> *f1, simpleFunction<double> *f2);
double computeLagNorm(int tag, simpleFunction<double> *f); double computeLagNorm(int tag, simpleFunction<double> *f);
// std::pair<PView *, PView*> buildErrorEstimateView // std::pair<PView *, PView*> buildErrorEstimateView
// (const std::string &errorFileName, double, int); // (const std::string &errorFileName, double, int);
......
Solver/functionSpace.h
+ 5
− 10
View file @ a2729cc1
...@@ -78,11 +78,14 @@ class FunctionSpaceBase ...@@ -78,11 +78,14 @@ class FunctionSpaceBase
template<class T> template<class T>
class FunctionSpace : public FunctionSpaceBase class FunctionSpace : public FunctionSpaceBase
{ {
protected:
int _iField; // field number (used to build dof keys)
public: public:
typedef typename TensorialTraits<T>::ValType ValType; typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType; typedef typename TensorialTraits<T>::GradType GradType;
typedef typename TensorialTraits<T>::HessType HessType; typedef typename TensorialTraits<T>::HessType HessType;
typedef typename TensorialTraits<T>::ThirdDevType ThirdDevType; typedef typename TensorialTraits<T>::ThirdDevType ThirdDevType;
virtual int getId(void) const {return _iField;}
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) = 0; virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) = 0;
virtual void fuvw(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) {} // should return to pure virtual once all is done. virtual void fuvw(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) {} // should return to pure virtual once all is done.
virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads) = 0; virtual void gradf(MElement *ele, double u, double v, double w, std::vector<GradType> &grads) = 0;
...@@ -103,9 +106,6 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double> ...@@ -103,9 +106,6 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
typedef TensorialTraits<double>::GradType GradType; typedef TensorialTraits<double>::GradType GradType;
typedef TensorialTraits<double>::HessType HessType; typedef TensorialTraits<double>::HessType HessType;
protected:
int _iField; // field number (used to build dof keys)
private: private:
virtual void getKeys(MVertex *ver, std::vector<Dof> &keys) virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
{ {
...@@ -113,8 +113,7 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double> ...@@ -113,8 +113,7 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
} }
public: public:
ScalarLagrangeFunctionSpaceOfElement(int i = 0) : _iField(i) {} ScalarLagrangeFunctionSpaceOfElement(int i = 0) { _iField = i; }
virtual int getId(void) const {return _iField;}
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{ {
if(ele->getParent()) { if(ele->getParent()) {
...@@ -207,17 +206,13 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double> ...@@ -207,17 +206,13 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
typedef TensorialTraits<double>::GradType GradType; typedef TensorialTraits<double>::GradType GradType;
typedef TensorialTraits<double>::HessType HessType; typedef TensorialTraits<double>::HessType HessType;
protected:
int _iField; // field number (used to build dof keys)
private: private:
virtual void getKeys(MVertex *ver, std::vector<Dof> &keys) virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
{ {
keys.push_back(Dof(ver->getNum(), _iField)); keys.push_back(Dof(ver->getNum(), _iField));
} }
public: public:
ScalarLagrangeFunctionSpace(int i = 0) : _iField(i) {} ScalarLagrangeFunctionSpace(int i = 0) { _iField = i; }
virtual int getId(void) const {return _iField;}
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{ {
if(ele->getParent()) ele = ele->getParent(); if(ele->getParent()) ele = ele->getParent();
......
Solver/solverAlgorithms.h
+ 1
− 1
View file @ a2729cc1
...@@ -65,7 +65,7 @@ template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term, ...@@ -65,7 +65,7 @@ template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term,
IntPt *GP; IntPt *GP;
int npts = integrator.getIntPoints(e, &GP); int npts = integrator.getIntPoints(e, &GP);
term.get(e, npts, GP, localMatrix); //localMatrix.print(); term.get(e, npts, GP, localMatrix); //localMatrix.print();
printf("local matrix size = %d %d\n", localMatrix.size1(), localMatrix.size2()); //printf("local matrix size = %d %d\n", localMatrix.size1(), localMatrix.size2());
shapeFcts.getKeys(e, R); shapeFcts.getKeys(e, R);
testFcts.getKeys(e, C); testFcts.getKeys(e, C);
// std::cout << "assembling normal test function ; lagrange trial function : " << std::endl; // std::cout << "assembling normal test function ; lagrange trial function : " << std::endl;
......
Solver/terms.cpp
+ 1
− 1
View file @ a2729cc1
...@@ -94,7 +94,7 @@ void IsotropicElasticTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<do ...@@ -94,7 +94,7 @@ void IsotropicElasticTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<do
} }
BTH.setAll(0.); BTH.setAll(0.);
BTH.gemm(BT, H); BTH.gemm(BT, H);
m.gemm(BTH, B, weight * detJ, 1.); m.gemm(BTH, B, weight * detJ, 1.); //m = m + w*detJ*BT*H*B
} }
} }
else else
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment