Skip to content
Snippets Groups Projects
Commit abf2c6cf authored by Éric Béchet's avatar Éric Béchet
Browse files

nix

parent f889b7c6
No related branches found
No related tags found
No related merge requests found
...@@ -429,21 +429,21 @@ void MyelasticitySolver::solve() ...@@ -429,21 +429,21 @@ void MyelasticitySolver::solve()
for (unsigned int i = 0; i < edgeNeu.size(); i++) for (unsigned int i = 0; i < edgeNeu.size(); i++)
{ {
LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,edgeNeu[i]._f); LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,edgeNeu[i]._f);
Assemble(Lterm,P123,edgeNeu[i].g->begin(),edgeNeu[i].g->end(),Integ_Boundary,*pAssembler); Assemble(Lterm,P123,edgeNeu[i].g->begin(),edgeNeu[i].g->end(),Integ_Boundary,*pAssembler);
printf("-- Force on edge %3d \n", edgeNeu[i]._tag); printf("-- Force on edge %3d \n", edgeNeu[i]._tag);
} }
for (unsigned int i = 0; i < faceNeu.size(); i++) for (unsigned int i = 0; i < faceNeu.size(); i++)
{ {
LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,faceNeu[i]._f); LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,faceNeu[i]._f);
Assemble(Lterm,P123,faceNeu[i].g->begin(),faceNeu[i].g->end(),Integ_Boundary,*pAssembler); Assemble(Lterm,P123,faceNeu[i].g->begin(),faceNeu[i].g->end(),Integ_Boundary,*pAssembler);
printf("-- Force on face %3d \n", faceNeu[i]._tag); printf("-- Force on face %3d \n", faceNeu[i]._tag);
} }
for (unsigned int i = 0; i < volumeNeu.size(); i++) for (unsigned int i = 0; i < volumeNeu.size(); i++)
{ {
LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,volumeNeu[i]._f); LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,volumeNeu[i]._f);
Assemble(Lterm,P123,volumeNeu[i].g->begin(),volumeNeu[i].g->end(),Integ_Boundary,*pAssembler); Assemble(Lterm,P123,volumeNeu[i].g->begin(),volumeNeu[i].g->end(),Integ_Boundary,*pAssembler);
printf("-- Force on volume %3d \n", volumeNeu[i]._tag); printf("-- Force on volume %3d \n", volumeNeu[i]._tag);
} }
...@@ -451,7 +451,7 @@ void MyelasticitySolver::solve() ...@@ -451,7 +451,7 @@ void MyelasticitySolver::solve()
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++)
{ {
ElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu); IsotropicElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
Assemble(Eterm,P123,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler); Assemble(Eterm,P123,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
} }
......
...@@ -19,7 +19,7 @@ ...@@ -19,7 +19,7 @@
#include "terms.h" #include "terms.h"
#include "quadratureRules.h" #include "quadratureRules.h"
#include "MVertex.h" #include "MVertex.h"
#include <iostream>
template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler) // symmetric template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler) // symmetric
{ {
......
...@@ -105,7 +105,7 @@ class ScalarTerm : public ScalarTermBase ...@@ -105,7 +105,7 @@ class ScalarTerm : public ScalarTermBase
template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // non symmetric template<class S1,class S2> class IsotropicElasticTerm : public BilinearTerm<S1,S2> // non symmetric
{ {
protected : protected :
double E,nu; double E,nu;
...@@ -118,7 +118,7 @@ template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // no ...@@ -118,7 +118,7 @@ template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // no
{ 0, 0, 0, 0, 0, C44} };*/ { 0, 0, 0, 0, 0, C44} };*/
public : public :
ElasticTerm(S1& space1_,S2& space2_,double E_,double nu_) : BilinearTerm<S1,S2>(space1_,space2_),E(E_),nu(nu_),H(6,6) IsotropicElasticTerm(S1& space1_,S2& space2_,double E_,double nu_) : BilinearTerm<S1,S2>(space1_,space2_),E(E_),nu(nu_),H(6,6)
{ {
double FACT = E / (1 + nu); double FACT = E / (1 + nu);
double C11 = FACT * (1 - nu) / (1 - 2 * nu); double C11 = FACT * (1 - nu) / (1 - 2 * nu);
...@@ -176,7 +176,7 @@ template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // no ...@@ -176,7 +176,7 @@ template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // no
template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric template<class S1> class IsotropicElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
{ {
protected : protected :
double E,nu; double E,nu;
...@@ -189,7 +189,7 @@ template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symm ...@@ -189,7 +189,7 @@ template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symm
{ 0, 0, 0, 0, 0, C44} };*/ { 0, 0, 0, 0, 0, C44} };*/
public : public :
ElasticTerm(S1& space1_,double E_,double nu_) : BilinearTerm<S1,S1>(space1_,space1_),E(E_),nu(nu_),H(6,6) IsotropicElasticTerm(S1& space1_,double E_,double nu_) : BilinearTerm<S1,S1>(space1_,space1_),E(E_),nu(nu_),H(6,6)
{ {
double FACT = E / (1 + nu); double FACT = E / (1 + nu);
double C11 = FACT * (1 - nu) / (1 - 2 * nu); double C11 = FACT * (1 - nu) / (1 - 2 * nu);
...@@ -235,36 +235,12 @@ template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symm ...@@ -235,36 +235,12 @@ template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symm
inline double dot(const double &a, const double &b) inline double dot(const double &a, const double &b)
{ return a*b; } { return a*b; }
template<class S1> class LoadTerm : public LinearTerm<S1>
{
typename S1::ValType Load;
public :
LoadTerm(S1& space1_,typename S1::ValType Load_) :LinearTerm<S1>(space1_),Load(Load_) {};
virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
{
double nbFF=LinearTerm<S1>::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];
const double weight = GP[i].weight;const double detJ = ele->getJacobian(u, v, w, jac);
std::vector<typename S1::ValType> Vals;
LinearTerm<S1>::space1.f(ele,u, v, w, Vals);
for (int j = 0; j < nbFF ; ++j)
{
m(j)+=dot(Vals[j],Load)*weight*detJ;
}
}
}
};
template<class S1> class LoadTermSF : public LinearTerm<S1> template<class S1> class LoadTerm : public LinearTerm<S1>
{ {
simpleFunction<typename S1::ValType> Load; simpleFunction<typename S1::ValType> Load;
public : public :
LoadTermSF(S1& space1_,simpleFunction<typename S1::ValType> Load_) :LinearTerm<S1>(space1_),Load(Load_) {}; LoadTerm(S1& space1_,simpleFunction<typename S1::ValType> Load_) :LinearTerm<S1>(space1_),Load(Load_) {};
virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m) virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
{ {
double nbFF=LinearTerm<S1>::space1.getNumKeys(ele); double nbFF=LinearTerm<S1>::space1.getNumKeys(ele);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment