From a2729cc1cd7f74644340e7c823ec8d5a246bda54 Mon Sep 17 00:00:00 2001
From: Gaetan Bricteux <gaetan.bricteux@uclouvain.be>
Date: Thu, 21 Apr 2016 08:51:25 +0000
Subject: [PATCH] enable multiple LM fields + postpro
---
Solver/elasticitySolver.cpp | 196 ++++++++++++++++++++++++++++--------
Solver/elasticitySolver.h | 19 ++--
Solver/functionSpace.h | 15 +--
Solver/solverAlgorithms.h | 2 +-
Solver/terms.cpp | 2 +-
5 files changed, 175 insertions(+), 59 deletions(-)
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index ceb4c2d9df..ace825b2e0 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -25,7 +25,6 @@
#include "PViewData.h"
#endif
-
/*
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->readMSH(meshFileName.c_str());
_dim = pModel->getNumRegions() ? 3 : 2;
if (LagSpace) delete LagSpace;
- if (_dim==3) LagSpace = new VectorLagrangeFunctionSpace(_tag);
- if (_dim==2) LagSpace = new VectorLagrangeFunctionSpace
+ if (dim == 3 || _dim==3) LagSpace = new VectorLagrangeFunctionSpace(_tag);
+ else if (dim == 2 || _dim==2) LagSpace = new VectorLagrangeFunctionSpace
(_tag,VectorLagrangeFunctionSpace::VECTOR_X,
VectorLagrangeFunctionSpace::VECTOR_Y);
- if (LagrangeMultiplierSpace) delete LagrangeMultiplierSpace;
- LagrangeMultiplierSpace = new ScalarLagrangeFunctionSpaceOfElement(_tag+1);
+ for (unsigned int i = 0; i < LagrangeMultiplierSpaces.size(); i++)
+ if (LagrangeMultiplierSpaces[i]) delete LagrangeMultiplierSpaces[i];
+ LagrangeMultiplierSpaces.clear();
}
void elasticitySolver::exportKb()
{
std::string sysname = "A";
+ if(!pAssembler ||
+ !pAssembler->getLinearSystem(sysname) ||
+ !pAssembler->getLinearSystem(sysname)->isAllocated()) return;
double valeur;
FILE *f = Fopen("K.txt", "w");
if(f){
@@ -95,8 +110,11 @@ void elasticitySolver::exportKb()
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)
linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#elif defined(HAVE_PETSC)
@@ -176,11 +194,12 @@ void elasticitySolver::readInputFile(const std::string &fn)
fclose(f);
return;
}
- field._tag = _tag;
- field._d = SVector3(d1, d2, d3);
+ SVector3 sv(d1, d2, d3);
+ field._d = sv.unit();
field._f = new simpleFunction<double>(val);
field.g = new groupOfElements (_dim - 1, physical);
LagrangeMultiplierFields.push_back(field);
+ LagrangeMultiplierSpaces.push_back(new ScalarLagrangeFunctionSpaceOfElement(field._tag));
}
else if (!strcmp(what, "Void")){
elasticField field;
@@ -368,9 +387,20 @@ void elasticitySolver::setLagrangeMultipliers(int phys, double tau, SVector3 d,
field._tau = tau;
field._tag = tag;
field._f = f;
- field._d = d;
- field.g = new groupOfElements (_dim - 1, phys);
+ field._d = d.unit();
+ field.g = new groupOfElements (_dim - 1, phys); // LM applied on entity of dimension = dim-1!
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)
@@ -445,19 +475,6 @@ void elasticitySolver::addElasticDomain(std::string phys, double e, double 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)
{
if (pAssembler) delete pAssembler;
@@ -475,7 +492,10 @@ void elasticitySolver::assemble(linearSystem<double> *lsys)
}
// LagrangeMultipliers
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);
}
// Elastic Fields
@@ -502,27 +522,31 @@ 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++){
+ unsigned int j = 0;
+ for ( ; j < LagrangeMultiplierSpaces.size(); j++)
+ if (LagrangeMultiplierSpaces[j]->getId() == LagrangeMultiplierFields[i]._tag) break;
printf("Lagrange Mult Lag\n");
- LagrangeMultiplierTerm<SVector3> LagTerm(*LagSpace, *LagrangeMultiplierSpace,
+ LagrangeMultiplierTerm<SVector3> LagTerm(*LagSpace, *(LagrangeMultiplierSpaces[j]),
LagrangeMultiplierFields[i]._d);
- Assemble(LagTerm, *LagSpace, *LagrangeMultiplierSpace,
+ Assemble(LagTerm, *LagSpace, *(LagrangeMultiplierSpaces[j]),
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(),
+ LaplaceTerm<double,double> LapTerm(*(LagrangeMultiplierSpaces[j]),
+ -LagrangeMultiplierFields[i]._tau);
+ Assemble(LapTerm, *(LagrangeMultiplierSpaces[j]), 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(),
+ LoadTermOnBorder<double> Lterm(*(LagrangeMultiplierSpaces[j]), LagrangeMultiplierFields[i]._f);
+ Assemble(Lterm, *(LagrangeMultiplierSpaces[j]), LagrangeMultiplierFields[i].g->begin(),
LagrangeMultiplierFields[i].g->end(), Integ_Boundary, *pAssembler);
}
// Assemble elastic term for
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); i++){
+ printf("Elastic\n");
IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
Assemble(Eterm, *LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),
Integ_Bulk, *pAssembler);
@@ -637,7 +661,8 @@ void elasticitySolver::computeEffectiveStrain(std::vector<double> strain)
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;
double err = 0.;
std::set<MVertex*> v;
@@ -664,8 +689,10 @@ double elasticitySolver::computeDisplacementError(int comp, simpleFunction<doubl
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;
+ SVector3 sol((*f0)((*it)->x(), (*it)->y(), (*it)->z()), (*f1)((*it)->x(), (*it)->y(), (*it)->z()),
+ (*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){
SVector3 val;
@@ -673,19 +700,95 @@ double elasticitySolver::computeDisplacementError(int comp, simpleFunction<doubl
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;
+ SVector3 sol((*f0)(xyz[0], xyz[1], xyz[2]), (*f1)(xyz[0], xyz[1], xyz[2]),
+ (*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);
}
+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) {
return 0;
}
#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)
{
std::cout << "build Displacement View"<< std::endl;
@@ -862,10 +965,14 @@ PView* elasticitySolver::buildStrainView (const std::string postFileName)
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;
- 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;
for (unsigned int i = 0; i < LagrangeMultiplierFields.size(); ++i){
for(groupOfElements::elementContainer::const_iterator it =
@@ -876,7 +983,7 @@ PView* elasticitySolver::buildLagrangeMultiplierView (const std::string postFile
}
}
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){
double val;
MPoint p(*it);
@@ -986,6 +1093,13 @@ PView *elasticitySolver::buildVonMisesView(const std::string postFileName)
#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)
{
Msg::Error("Post-pro module not available");
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 7419c637f1..ff3fd67512 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -63,7 +63,7 @@ class elasticitySolver
int _dim, _tag;
dofManager<double> *pAssembler;
FunctionSpace<SVector3> *LagSpace;
- FunctionSpace<double> *LagrangeMultiplierSpace;
+ std::vector<FunctionSpace<double> *> LagrangeMultiplierSpaces;
// young modulus and poisson coefficient per physical
std::vector<elasticField> elasticFields;
@@ -75,7 +75,7 @@ class elasticitySolver
std::vector<dirichletBC> allDirichlet;
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);
@@ -89,16 +89,19 @@ class elasticitySolver
virtual ~elasticitySolver()
{
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;
}
void assemble (linearSystem<double> *lsys);
void readInputFile(const std::string &meshFileName);
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 setElasticDomain(int phys, double E, double nu);
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 solve();
void postSolve();
@@ -108,11 +111,15 @@ class elasticitySolver
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 *buildLagrangeMultiplierView(const std::string posFileName, int tag = -1);
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);
+ 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);
// std::pair<PView *, PView*> buildErrorEstimateView
// (const std::string &errorFileName, double, int);
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 4d3e7a2184..5560f086ff 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -78,11 +78,14 @@ class FunctionSpaceBase
template<class T>
class FunctionSpace : public FunctionSpaceBase
{
+ protected:
+ int _iField; // field number (used to build dof keys)
public:
typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType;
typedef typename TensorialTraits<T>::HessType HessType;
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 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;
@@ -103,9 +106,6 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
typedef TensorialTraits<double>::GradType GradType;
typedef TensorialTraits<double>::HessType HessType;
- protected:
- int _iField; // field number (used to build dof keys)
-
private:
virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
{
@@ -113,8 +113,7 @@ class ScalarLagrangeFunctionSpaceOfElement : public FunctionSpace<double>
}
public:
- ScalarLagrangeFunctionSpaceOfElement(int i = 0) : _iField(i) {}
- virtual int getId(void) const {return _iField;}
+ ScalarLagrangeFunctionSpaceOfElement(int i = 0) { _iField = i; }
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
if(ele->getParent()) {
@@ -207,17 +206,13 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
typedef TensorialTraits<double>::GradType GradType;
typedef TensorialTraits<double>::HessType HessType;
- protected:
- int _iField; // field number (used to build dof keys)
-
private:
virtual void getKeys(MVertex *ver, std::vector<Dof> &keys)
{
keys.push_back(Dof(ver->getNum(), _iField));
}
public:
- ScalarLagrangeFunctionSpace(int i = 0) : _iField(i) {}
- virtual int getId(void) const {return _iField;}
+ ScalarLagrangeFunctionSpace(int i = 0) { _iField = i; }
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
if(ele->getParent()) ele = ele->getParent();
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index f7be9c8c2c..445541dedd 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -65,7 +65,7 @@ template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term,
IntPt *GP;
int npts = integrator.getIntPoints(e, &GP);
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);
testFcts.getKeys(e, C);
// std::cout << "assembling normal test function ; lagrange trial function : " << std::endl;
diff --git a/Solver/terms.cpp b/Solver/terms.cpp
index 016b6fd63a..687c09e6ce 100644
--- a/Solver/terms.cpp
+++ b/Solver/terms.cpp
@@ -94,7 +94,7 @@ void IsotropicElasticTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<do
}
BTH.setAll(0.);
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
--
GitLab