diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp index a4a70c44bcd8d52ec4af757ec3298416230d9b46..eb86ec0d4b697e34d1f9bf5123d79940085a7089 100644 --- a/Solver/elasticitySolver.cpp +++ b/Solver/elasticitySolver.cpp @@ -198,7 +198,7 @@ void elasticitySolver::solve() // LaplaceTerm<SVector3,SVector3> Eterm(*LagSpace); Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler); } - + printf("nDofs=%d\n",pAssembler->sizeOfR()); printf("-- done assembling!\n"); lsys->systemSolve(); printf("-- done solving!\n"); diff --git a/Solver/terms.h b/Solver/terms.h index 837ed4f126515826b2e5c63b5631f4e248000424..e7d22cb2937a098a17372921de89709ba1c417cf 100644 --- a/Solver/terms.h +++ b/Solver/terms.h @@ -208,7 +208,8 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3> virtual ~IsotropicElasticTerm() {} virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) { - if (sym) + if (ele->getParent()) ele=ele->getParent(); + if (sym) { int nbFF = BilinearTerm<SVector3,SVector3>::space1.getNumKeys(ele); double jac[3][3]; @@ -222,7 +223,7 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3> 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]; - if (ele->getParent()) ele=ele->getParent(); + const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac); std::vector<TensorialTraits<SVector3>::GradType> Grads; BilinearTerm<SVector3,SVector3>::space1.gradf(ele,u, v, w, Grads); // a optimiser ?? @@ -300,7 +301,8 @@ template<class T1> class LoadTerm : public LinearTerm<T1> virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m) { - int nbFF=LinearTerm<T1>::space1.getNumKeys(ele); + if (ele->getParent()) ele=ele->getParent(); + int nbFF=LinearTerm<T1>::space1.getNumKeys(ele); double jac[3][3]; m.resize(nbFF); m.scale(0.);