diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index b936f23343264acdb8672d0fb083da46fea26e89..467bfbb270a2866a0014bf08b96ee3757a098291 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -19,6 +19,8 @@ set(SRC
functionSpace.cpp
filters.cpp
sparsityPattern.cpp
+ STensor43.cpp
+ terms.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
diff --git a/Solver/STensor43.cpp b/Solver/STensor43.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cb46e4ebc1bd4ed8533c438565cd3eb14083a48b
--- /dev/null
+++ b/Solver/STensor43.cpp
@@ -0,0 +1,11 @@
+
+#include "STensor43.h"
+
+void STensor43::print (const char *s) const
+{
+ printf("STensor43::print to be implemented \n");
+/* printf(" tensor4 %s : \n %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n %12.5E %12.5E %12.5E \n",s,
+ (*this)(0,0),(*this)(0,1),(*this)(0,2),
+ (*this)(1,0),(*this)(1,1),(*this)(1,2),
+ (*this)(2,0),(*this)(2,1),(*this)(2,2));*/
+}
diff --git a/Solver/STensor43.h b/Solver/STensor43.h
new file mode 100644
index 0000000000000000000000000000000000000000..dd461ed8b671c0dd35a60d9a2e1bb32abfb3e907
--- /dev/null
+++ b/Solver/STensor43.h
@@ -0,0 +1,130 @@
+#ifndef _STENSOR43_H_
+#define _STENSOR43_H_
+
+#include "STensor3.h"
+#include "fullMatrix.h"
+#include "Numeric.h"
+
+
+// concrete class for general 3x3 matrix
+
+class STensor43 {
+ protected:
+ // 0000 0001 0002 0010 ... 2211 2212 2220 2221 2222
+ double _val[81];
+ public:
+ inline int getIndex(int i, int j, int k, int l) const
+ {
+ static int _index[3][3][3][3] = {{{{0,1,2},{3,4,5},{6,7,8}},{{9,10,11},{12,13,14},{15,16,17}},{{18,19,20},{21,22,23},{24,25,26}}},
+ {{{27,28,29},{30,31,32},{33,34,35}},{{36,37,38},{39,40,41},{42,43,44}},{{45,46,47},{48,49,50},{51,52,53}}},
+ {{{54,55,56},{57,58,59},{60,61,62}},{{63,64,65},{66,67,68},{69,70,71}},{{72,73,74},{75,76,77},{77,79,80}}}};
+ return _index[i][j][k][l];
+ }
+ STensor43(const STensor43& other)
+ {
+ for (int i = 0; i < 81; i++) _val[i] = other._val[i];
+ }
+ // default constructor, null tensor
+ STensor43(const double v = 0.0)
+ {
+ for (int i = 0; i < 3; i++)
+ for (int j = 0; j < 3; j++)
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ if ((i==k)&&(j==l))
+ _val[getIndex(i, j, k, l)]=v;
+ else
+ _val[getIndex(i, j, k, l)]=0.0;
+ }
+ inline double &operator()(int i, int j,int k, int l)
+ {
+ return _val[getIndex(i, j, k, l)];
+ }
+ inline double operator()(int i, int j, int k, int l) const
+ {
+ return _val[getIndex(i, j, k ,l)];
+ }
+ STensor43 operator + (const STensor43 &other) const
+ {
+ STensor43 res(*this);
+ for (int i = 0; i < 81; i++) res._val[i] += other._val[i];
+ return res;
+ }
+ STensor43& operator += (const STensor43 &other)
+ {
+ for (int i = 0; i < 81; i++) _val[i] += other._val[i];
+ return *this;
+ }
+ STensor43& operator *= (const double &other)
+ {
+ for (int i = 0; i < 81; i++) _val[i] *= other;
+ return *this;
+ }
+/* STensor43& operator *= (const STensor43 &other)
+ {
+// to be implemented
+ return *this;
+ }*/
+ void print(const char *) const;
+};
+
+// tensor product
+inline void tensprod(const STensor3 &a, const STensor3 &b, STensor43 &c)
+{
+ for (int i = 0; i < 3; i++)
+ for (int j = 0; j < 3; j++)
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ c(i,j,k,l)=a(i,j)*b(k,l);
+}
+
+inline double dot(const STensor43 &a, const STensor43 &b)
+{
+ double prod=0;
+ for (int i = 0; i < 3; i++)
+ for (int j = 0; j < 3; j++)
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ prod+=a(i,j,k,l)*b(i,j,k,l);
+ return prod;
+}
+
+inline STensor43 operator*(const STensor43 &t, double m)
+{
+ STensor43 val(t);
+ val *= m;
+ return val;
+}
+
+inline STensor43 operator*(double m,const STensor43 &t)
+{
+ STensor43 val(t);
+ val *= m;
+ return val;
+}
+
+inline STensor3 operator*(const STensor43 &t, const STensor3 &m)
+{
+ STensor3 val(0.);
+ for (int i = 0; i < 3; i++)
+ for (int j = 0; j < 3; j++)
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ val(i,j)+=t(i,j,k,l)*m(k,l);
+ return val;
+}
+
+inline STensor3 operator*( const STensor3 &m , const STensor43 &t)
+{
+ STensor3 val(0.);
+ for (int i = 0; i < 3; i++)
+ for (int j = 0; j < 3; j++)
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ val(k,l)+=t(i,j,k,l)*m(i,j);
+ return val;
+}
+
+
+
+#endif
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 31cbea5fe5956271037b87a6fa97f74a408afc24..1fdd6715bd3d0bdb4e7673acecac62737a246e80 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -92,7 +92,7 @@ void elasticitySolver::solve()
{
SolverField<SVector3> Field(pAssembler, LagSpace);
IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
- BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+ BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
}
printf("elastic energy=%f\n",energ);
@@ -107,7 +107,7 @@ void elasticitySolver::postSolve()
{
SolverField<SVector3> Field(pAssembler, LagSpace);
IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
- BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+ BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
}
printf("elastic energy=%f\n",energ);
@@ -551,7 +551,7 @@ PView* elasticitySolver::buildStressesView (const std::string postFileName)
double nu = elasticFields[i]._nu;
SolverField<SVector3> Field(pAssembler, LagSpace);
IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
- BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+ BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
{
MElement *e=*it;
@@ -644,8 +644,8 @@ PView *elasticitySolver::buildElasticEnergyView(const std::string postFileName)
if(elasticFields[i]._E == 0.) continue;
SolverField<SVector3> Field(pAssembler, LagSpace);
IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
- BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
- ScalarTermConstant One(1.0);
+ BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
+ ScalarTermConstant<double> One(1.0);
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
{
MElement *e = *it;
@@ -673,7 +673,7 @@ PView *elasticitySolver::buildVonMisesView(const std::string postFileName)
{
SolverField<SVector3> Field(pAssembler, LagSpace);
IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
- BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+ BilinearTermToScalarTerm Elastic_Energy_Term(Eterm);
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
{
MElement *e=*it;
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index ec0da235b6a2e8d73a8dde19cb4859ed06dc51d0..68c30abedfc9b98895e79734d27dcd176b485641 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -3,6 +3,7 @@
#include "SVector3.h"
#include "STensor3.h"
+#include "STensor43.h"
#include <vector>
#include <iterator>
#include <iostream>
@@ -16,8 +17,8 @@ template<class T> struct TensorialTraits
{
typedef T ValType;
typedef T GradType[3];
-/* typedef T HessType[3][3];
- typedef SVoid DivType;
+ typedef T HessType[3][3];
+/* typedef SVoid DivType;
typedef SVoid CurlType;*/
};
@@ -26,6 +27,7 @@ template<> struct TensorialTraits<double>
typedef double ValType;
typedef SVector3 GradType;
typedef STensor3 HessType;
+ typedef double TensProdType;
/* typedef SVoid DivType;
typedef SVoid CurlType;*/
};
@@ -35,10 +37,23 @@ template<> struct TensorialTraits<SVector3>
typedef SVector3 ValType;
typedef STensor3 GradType;
typedef STensor3 HessType;
- typedef double DivType;
- typedef SVector3 CurlType;
+ typedef STensor3 TensProdType;
+// typedef double DivType;
+// typedef SVector3 CurlType;
};
+template<> struct TensorialTraits<STensor3>
+{
+ typedef STensor3 ValType;
+// typedef STensor3 GradType;
+// typedef STensor3 HessType;
+// typedef STensor3 TensProdType;
+ typedef STensor43 TensProdType;
+// typedef double DivType;
+// typedef SVector3 CurlType;
+};
+
+
class FunctionSpaceBase
{
public:
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index 3b3f7bbfac21dc4eebd1968addfd32fd441df2a2..2126746254711822d93536f31cc7ce80e22edaac 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -97,7 +97,7 @@ template<class Iterator, class Assembler> void Assemble(BilinearTermBase &term,
-template<class Iterator, class Assembler> void Assemble(LinearTermBase &term, FunctionSpaceBase &space,
+template<class Iterator, class Assembler> void Assemble(LinearTermBase<double> &term, FunctionSpaceBase &space,
Iterator itbegin, Iterator itend,
QuadratureBase &integrator, Assembler &assembler)
{
@@ -114,7 +114,7 @@ template<class Iterator, class Assembler> void Assemble(LinearTermBase &term, Fu
}
}
-template<class Assembler> void Assemble(LinearTermBase &term, FunctionSpaceBase &space, MElement *e,
+template<class Assembler> void Assemble(LinearTermBase<double> &term, FunctionSpaceBase &space, MElement *e,
QuadratureBase &integrator, Assembler &assembler)
{
fullVector<typename Assembler::dataMat> localVector;
@@ -126,7 +126,7 @@ template<class Assembler> void Assemble(LinearTermBase &term, FunctionSpaceBase
assembler.assemble(R, localVector);
}
-template<class Iterator, class dataMat> void Assemble(ScalarTermBase &term,
+template<class Iterator, class dataMat> void Assemble(ScalarTermBase<double> &term,
Iterator itbegin, Iterator itend,
QuadratureBase &integrator, dataMat & val)
{
@@ -140,7 +140,7 @@ template<class Iterator, class dataMat> void Assemble(ScalarTermBase &term,
}
}
-template<class Iterator, class dataMat> void Assemble(ScalarTermBase &term, MElement *e,
+template<class Iterator, class dataMat> void Assemble(ScalarTermBase<double> &term, MElement *e,
QuadratureBase &integrator, dataMat & val)
{
dataMat localval;
diff --git a/Solver/terms.cpp b/Solver/terms.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..08622668b3dcae207854fbb68607b07358c31ce3
--- /dev/null
+++ b/Solver/terms.cpp
@@ -0,0 +1,198 @@
+//
+// C++ Implementation: terms
+//
+// Description:
+//
+//
+// Author: <Eric Bechet>, (C) 2011
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "terms.h"
+
+
+void BilinearTermToScalarTerm::get(MElement *ele, int npts, IntPt *GP, double &val) const
+{
+ fullMatrix<double> localMatrix;
+ bilterm.get(ele, npts, GP, localMatrix);
+ val = localMatrix(0, 0);
+}
+
+void BilinearTermBase::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+ std::vector<fullMatrix<double> > mv(npts);
+ get(ele,npts,GP,mv);
+ m.resize(mv[0].size1(), mv[0].size2());
+ m.setAll(0.);
+ double jac[3][3];
+ for (int k=0;k<npts;k++)
+ {
+ const double u = GP[k].pt[0]; const double v = GP[k].pt[1]; const double w = GP[k].pt[2];
+ const double weight = GP[k].weight; const double detJ = ele->getJacobian(u, v, w, jac);
+ const double coeff=weight*detJ;
+ for (int i=0;i<mv[k].size1();++i)
+ for (int j=0;j<mv[k].size2();++j)
+ m(i,j)+=mv[k](i,j)*coeff;
+ }
+}
+
+IsotropicElasticTerm::IsotropicElasticTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<SVector3>& space2_, double E_, double nu_) :
+ BilinearTerm<SVector3, SVector3>(space1_, space2_), E(E_), nu(nu_), H(6, 6)
+{
+ double FACT = E / (1 + nu);
+ double C11 = FACT * (1 - nu) / (1 - 2 * nu);
+ double C12 = FACT * nu / (1 - 2 * nu);
+ double C44 = (C11 - C12) / 2;
+ H.scale(0.);
+ for(int i = 0; i < 3; ++i) { H(i, i) = C11; H(i + 3, i + 3) = C44; }
+ H(1, 0) = H(0, 1) = H(2, 0) = H(0, 2) = H(1, 2) = H(2, 1) = C12;
+ sym = (&space1_ == &space2_);
+}
+IsotropicElasticTerm::IsotropicElasticTerm(FunctionSpace<SVector3>& space1_, double E_, double nu_) :
+ BilinearTerm<SVector3, SVector3>(space1_, space1_), E(E_), nu(nu_), H(6, 6)
+{
+ double FACT = E / (1 + nu);
+ double C11 = FACT * (1 - nu) / (1 - 2 * nu);
+ double C12 = FACT * nu / (1 - 2 * nu);
+ double C44 = (C11 - C12) / 2;
+/* FACT = E / (1 - nu * nu); // plane stress (plates)
+ C11 = FACT;
+ C12 = nu * FACT;
+ C44 = (1. - nu) * .5 * FACT;*/
+ H.scale(0.);
+ for(int i = 0; i < 3; ++i) { H(i, i) = C11; H(i + 3, i + 3) = C44; }
+ H(1, 0) = H(0, 1) = H(2, 0) = H(0, 2) = H(1, 2) = H(2, 1) = C12;
+ sym = true;
+}
+
+void IsotropicElasticTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+ if(ele->getParent()) ele = ele->getParent();
+ if(sym)
+ {
+ int nbFF = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele);
+ double jac[3][3];
+ fullMatrix<double> B(6, nbFF);
+ fullMatrix<double> BTH(nbFF, 6);
+ fullMatrix<double> BT(nbFF, 6);
+ m.resize(nbFF, nbFF);
+ m.setAll(0.);
+ //std::cout << m.size1() << " " << m.size2() << std::endl;
+ 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<TensorialTraits<SVector3>::GradType> Grads;
+ BilinearTerm<SVector3, SVector3>::space1.gradf(ele, u, v, w, Grads); // a optimiser ??
+ for(int j = 0; j < nbFF; j++)
+ {
+ BT(j, 0) = B(0, j) = Grads[j](0, 0);
+ BT(j, 1) = B(1, j) = Grads[j](1, 1);
+ BT(j, 2) = B(2, j) = Grads[j](2, 2);
+ BT(j, 3) = B(3, j) = Grads[j](0, 1) + Grads[j](1, 0);
+ BT(j, 4) = B(4, j) = Grads[j](1, 2) + Grads[j](2, 1);
+ BT(j, 5) = B(5, j) = Grads[j](0, 2) + Grads[j](2, 0);
+ }
+ BTH.setAll(0.);
+ BTH.gemm(BT, H);
+ m.gemm(BTH, B, weight * detJ, 1.);
+ }
+ }
+ else
+ {
+ int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele);
+ int nbFF2 = BilinearTerm<SVector3, SVector3>::space2.getNumKeys(ele);
+ double jac[3][3];
+ fullMatrix<double> B(6, nbFF2);
+ fullMatrix<double> BTH(nbFF2, 6);
+ fullMatrix<double> BT(nbFF1, 6);
+ m.resize(nbFF1, nbFF2);
+ m.setAll(0.);
+ // Sum on Gauss Points i
+ 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<TensorialTraits<SVector3>::GradType> Grads;// tableau de matrices...
+ std::vector<TensorialTraits<SVector3>::GradType> GradsT;// tableau de matrices...
+ BilinearTerm<SVector3, SVector3>::space1.gradf(ele, u, v, w, Grads);
+ BilinearTerm<SVector3, SVector3>::space2.gradf(ele, u, v, w, GradsT);
+ for(int j = 0; j < nbFF1; j++)
+ {
+ BT(j, 0) = Grads[j](0, 0);
+ BT(j, 1) = Grads[j](1, 1);
+ BT(j, 2) = Grads[j](2, 2);
+ BT(j, 3) = Grads[j](0, 1) + Grads[j](1, 0);
+ BT(j, 4) = Grads[j](1, 2) + Grads[j](2, 1);
+ BT(j, 5) = Grads[j](0, 2) + Grads[j](2, 0);
+ }
+ for(int j = 0; j < nbFF2; j++)
+ {
+ B(0, j) = GradsT[j](0, 0);
+ B(1, j) = GradsT[j](1, 1);
+ B(2, j) = GradsT[j](2, 2);
+ B(3, j) = GradsT[j](0, 1) + GradsT[j](1, 0);
+ B(4, j) = GradsT[j](1, 2) + GradsT[j](2, 1);
+ B(5, j) = GradsT[j](0, 2) + GradsT[j](2, 0);
+ }
+ BTH.setAll(0.);
+ BTH.gemm(BT, H);
+ // gemm add the product to m so there is a sum on gauss' points here
+ m.gemm(BTH, B, weight * detJ, 1.);
+ }
+ }
+}
+
+void LagrangeMultiplierTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+ int nbFF1 = BilinearTerm<SVector3, double>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
+ int nbFF2 = BilinearTerm<SVector3, double>::space2.getNumKeys(ele); //nbVertices of boundary
+ double jac[3][3];
+ m.resize(nbFF1, nbFF2);
+ m.setAll(0.);
+ for(int i = 0; i < npts; i++)
+ {
+ double u = GP[i].pt[0]; double v = GP[i].pt[1]; double w = GP[i].pt[2];
+ const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
+ std::vector<TensorialTraits<SVector3>::ValType> Vals;
+ std::vector<TensorialTraits<double>::ValType> ValsT;
+ BilinearTerm<SVector3,double>::space1.f(ele, u, v, w, Vals);
+ BilinearTerm<SVector3,double>::space2.f(ele, u, v, w, ValsT);
+ for(int j = 0; j < nbFF1; j++)
+ {
+ for(int k = 0; k < nbFF2; k++)
+ {
+ m(j, k) += dot(Vals[j], _d) * ValsT[k] * weight * detJ;
+ }
+ }
+ }
+}
+
+void LagMultTerm::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+ int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
+ int nbFF2 = BilinearTerm<SVector3, SVector3>::space2.getNumKeys(ele); //nbVertices of boundary
+ double jac[3][3];
+ m.resize(nbFF1, nbFF2);
+ m.setAll(0.);
+ for(int i = 0; i < npts; i++)
+ {
+ double u = GP[i].pt[0]; double v = GP[i].pt[1]; double w = GP[i].pt[2];
+ const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
+ std::vector<TensorialTraits<SVector3>::ValType> Vals;
+ std::vector<TensorialTraits<SVector3>::ValType> ValsT;
+ BilinearTerm<SVector3,SVector3>::space1.f(ele, u, v, w, Vals);
+ BilinearTerm<SVector3,SVector3>::space2.f(ele, u, v, w, ValsT);
+ for(int j = 0; j < nbFF1; j++)
+ {
+ for(int k = 0; k < nbFF2; k++)
+ {
+ m(j, k) += _eqfac * dot(Vals[j], ValsT[k]) * weight * detJ;
+ }
+ }
+ }
+}
+
+
diff --git a/Solver/terms.h b/Solver/terms.h
index a336c2490a9af709e9d24ae7e36d1eafee2c6178..e3540ec8ee36df169806b3c2895e32e2bb8d5390 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -4,7 +4,7 @@
// Description:
//
//
-// Author: <Eric Bechet>, (C) 2009
+// Author: <Eric Bechet>, (C) 2011
//
// Copyright: See COPYING file that comes with this distribution
//
@@ -15,102 +15,181 @@
#define _TERMS_H_
#include "SVector3.h"
-#include <vector>
-#include <iterator>
+#include "STensor3.h"
+#include "STensor43.h"
#include "Numeric.h"
#include "functionSpace.h"
#include "groupOfElements.h"
#include "materialLaw.h"
+#include <vector>
+#include <iterator>
+
+template<class T2> class ScalarTermBase;
+
+template<class T2> class LinearTermBase;
+template<class T2> class PlusTerm;
+
+class BilinearTermBase;
+
+inline double dot(const double &a, const double &b)
+{
+ return a * b;
+}
+
+inline int delta(int i,int j) {if (i==j) return 1; else return 0;}
+
+
+
+
+template<class T2=double> class ScalarTermBase
+{
+public :
+ virtual ~ScalarTermBase() {}
+ virtual void get(MElement *ele, int npts, IntPt *GP, T2 &val) const = 0 ;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<T2> &vval) const = 0 ;
+ virtual ScalarTermBase<T2>* clone () const =0;
+};
+
+template<class T2=double> class ScalarTerm : public ScalarTermBase<T2>
+{
+public :
+ virtual ~ScalarTerm() {}
+};
+
+
+template<class T2=double> class ScalarTermConstant : public ScalarTerm<T2>
+{
+protected:
+ T2 cst;
+public :
+ ScalarTermConstant(T2 val_ = T2()) : cst(val_) {}
+ virtual ~ScalarTermConstant() {}
+ virtual void get(MElement *ele, int npts, IntPt *GP, T2 &val) const;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<T2> &vval) const;
+ virtual void get(MVertex *ver, T2 &val) const;
+ virtual ScalarTermBase<T2>* clone () const {return new ScalarTermConstant<T2>(cst);}
+};
+
+class BilinearTermToScalarTerm : public ScalarTerm<double>
+{
+protected:
+ BilinearTermBase &bilterm;
+public :
+ BilinearTermToScalarTerm(BilinearTermBase &bilterm_): bilterm(bilterm_){}
+ virtual ~BilinearTermToScalarTerm() {}
+ virtual void get(MElement *ele, int npts, IntPt *GP, double &val) const;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<double> &vval) const {};
+ virtual ScalarTermBase<double>* clone () const {return new BilinearTermToScalarTerm(bilterm);}
+};
+
+
class BilinearTermBase
{
- public :
+public :
virtual ~BilinearTermBase() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) = 0 ;
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const ;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullMatrix<double> > &mv) const = 0 ;
+ virtual BilinearTermBase* clone () const =0;
+};
+
+
+template<class T2> class BilinearTermContract : public BilinearTermBase
+{
+protected:
+ const LinearTermBase<T2> *a;
+ const LinearTermBase<T2> *b;
+public:
+ BilinearTermContract(const LinearTermBase<T2> &a_,const LinearTermBase<T2> &b_) : a(a_.clone()),b(b_.clone()) {}
+ virtual ~BilinearTermContract() {delete a;delete b;}
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const ;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullMatrix<double> > &mv) const {};
+ virtual BilinearTermBase* clone () const {return new BilinearTermContract<T2>(*a,*b);}
};
+template<class T2> class BilinearTermContractWithLaw : public BilinearTermContract<T2>
+{
+public:
+protected:
+ const ScalarTermBase< typename TensorialTraits<T2>::TensProdType > *c;
+public:
+ BilinearTermContractWithLaw(const LinearTermBase<T2> &a_,const ScalarTermBase<typename TensorialTraits<T2>::TensProdType> &c_, const LinearTermBase<T2> &b_) : BilinearTermContract<T2>(a_,b_), c(c_.clone()) {}
+ virtual ~BilinearTermContractWithLaw() {delete c;}
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const ;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullMatrix<double> > &mv) const ;
+ virtual BilinearTermBase* clone () const {return new BilinearTermContractWithLaw<T2>(*BilinearTermContract<T2>::a,*c,*BilinearTermContract<T2>::b);}
+};
+
+template<class T2> BilinearTermContract<T2> operator |(const LinearTermBase<T2>& L1,const LinearTermBase<T2>& L2);
+
+
template<class T1,class T2> class BilinearTerm : public BilinearTermBase
{
- protected :
+protected :
FunctionSpace<T1>& space1;
FunctionSpace<T2>& space2;
- public :
+public :
BilinearTerm(FunctionSpace<T1>& space1_, FunctionSpace<T2>& space2_) : space1(space1_), space2(space2_) {}
virtual ~BilinearTerm() {}
};
-class LinearTermBase
+template<class T2=double> class LinearTermBase
{
- public:
+public:
virtual ~LinearTermBase() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<double> &v) = 0;
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &v) const ;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<T2> > &vv) const =0;
+ virtual LinearTermBase<T2>* clone () const = 0;
+ PlusTerm<T2> operator +(const LinearTermBase<T2>& other);
+};
+
+
+template<class T2> class PlusTerm:public LinearTermBase<T2>
+{
+ const LinearTermBase<T2> *a;
+ const LinearTermBase<T2> *b;
+public:
+ PlusTerm(const LinearTermBase<T2> &a_,const LinearTermBase<T2> &b_) : a(a_.clone()),b(b_.clone()) {}
+ virtual ~PlusTerm() {delete a;delete b;}
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &v) const ;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<T2> > &vv) const {};
+ virtual LinearTermBase<T2>* clone () const { return new PlusTerm<T2>(*a,*b);}
};
-template<class T1> class LinearTerm : public LinearTermBase
+
+template<class T1, class T2=double> class LinearTerm : public LinearTermBase<T2>
{
- protected :
+protected :
FunctionSpace<T1>& space1;
- public :
+public :
LinearTerm(FunctionSpace<T1>& space1_) : space1(space1_) {}
virtual ~LinearTerm() {}
};
-class ScalarTermBase
-{
- public :
- virtual ~ScalarTermBase() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, double &val) = 0;
-};
-class ScalarTerm : public ScalarTermBase
-{
- public :
- virtual ~ScalarTerm() {}
-};
-template<class T1, class T2> class BilinearTermToScalarTerm : public ScalarTerm
+template<class T1> class GradTerm : public LinearTerm<T1, typename TensorialTraits<T1>::GradType >
{
- BilinearTerm<T1, T2> &bilterm;
- public :
- BilinearTermToScalarTerm(BilinearTerm<T1, T2> &bilterm_): bilterm(bilterm_){}
- virtual ~BilinearTermToScalarTerm() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, double &val)
- {
- fullMatrix<double> localMatrix;
- bilterm.get(ele, npts, GP, localMatrix);
- val = localMatrix(0, 0);
- }
+public :
+ GradTerm(FunctionSpace<T1>& space1_) : LinearTerm<T1,typename TensorialTraits<T1>::GradType >(space1_) {}
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<typename TensorialTraits<T1>::GradType > &vec) const {LinearTermBase<typename TensorialTraits<T1>::GradType>::get(ele,npts,GP,vec);}
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<typename TensorialTraits<T1>::GradType > > &vvec) const;
+ virtual LinearTermBase<typename TensorialTraits<T1>::GradType>* clone () const { return new GradTerm<T1>(LinearTerm<T1,typename TensorialTraits<T1>::GradType>::space1);}
+ virtual ~GradTerm() {}
};
-class ScalarTermConstant : public ScalarTerm
-{
- double val;
- public :
- ScalarTermConstant(double val_ = 1.0) : val(val_) {}
- virtual ~ScalarTermConstant() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, double &val)
- {
- double jac[3][3];
- val = 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);
- val += weight * detJ;
- }
- }
- virtual void get(MVertex *ver, double &val)
- {
- val = 1;
- }
-};
template<class T1, class T2> class LaplaceTerm : public BilinearTerm<T1, T2>
{
- public :
+public :
LaplaceTerm(FunctionSpace<T1>& space1_, FunctionSpace<T2>& space2_) : BilinearTerm<T1, T2>(space1_, space2_)
{}
virtual ~LaplaceTerm() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+ {
+ Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
+ }
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< fullMatrix<double> > &vm) const
{
Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
}
@@ -118,40 +197,25 @@ template<class T1, class T2> class LaplaceTerm : public BilinearTerm<T1, T2>
{
Msg::Error("LaplaceTerm<S1, S2> w/ S1 != S2 not implemented");
}
+ virtual BilinearTermBase* clone () const {return new LaplaceTerm< T1, T2 >(BilinearTerm<T1, T2>::space1,BilinearTerm<T1, T2>::space2);}
}; // class
template<class T1> class LaplaceTerm<T1, T1> : public BilinearTerm<T1, T1> // symmetric
{
+protected:
double diffusivity;
- public :
+public :
LaplaceTerm(FunctionSpace<T1>& space1_, double diff = 1) :
BilinearTerm<T1, T1>(space1_, space1_), diffusivity(diff) {}
virtual ~LaplaceTerm() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
- {
- int nbFF = BilinearTerm<T1, T1>::space1.getNumKeys(ele);
- double jac[3][3];
- m.resize(nbFF, nbFF);
- m.setAll(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 TensorialTraits<T1>::GradType> Grads;
- BilinearTerm<T1, T1>::space1.gradf(ele, u, v, w, Grads);
- for(int j = 0; j < nbFF; j++){
- for(int k = j; k < nbFF; k++){
- double contrib = weight * detJ * dot(Grads[j], Grads[k]) * diffusivity;
- m(j, k) += contrib;
- if(j != k) m(k, j) += contrib;
- }
- }
- }
- }
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< fullMatrix<double> > &vm) const{};
+ virtual BilinearTermBase* clone () const {return new LaplaceTerm< T1, T1 >(BilinearTerm<T1, T1>::space1,diffusivity);}
}; // class
class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
{
- protected :
+protected :
double E, nu;
bool sym;
fullMatrix<double> H;/* =
@@ -161,237 +225,71 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
{ 0, 0, 0, C44, 0, 0},
{ 0, 0, 0, 0, C44, 0},
{ 0, 0, 0, 0, 0, C44} };*/
- public :
- IsotropicElasticTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<SVector3>& space2_, double E_, double nu_) :
- BilinearTerm<SVector3, SVector3>(space1_, space2_), E(E_), nu(nu_), H(6, 6)
- {
- double FACT = E / (1 + nu);
- double C11 = FACT * (1 - nu) / (1 - 2 * nu);
- double C12 = FACT * nu / (1 - 2 * nu);
- double C44 = (C11 - C12) / 2;
- H.scale(0.);
- for(int i = 0; i < 3; ++i) { H(i, i) = C11; H(i + 3, i + 3) = C44; }
- H(1, 0) = H(0, 1) = H(2, 0) = H(0, 2) = H(1, 2) = H(2, 1) = C12;
- sym = (&space1_ == &space2_);
- }
- IsotropicElasticTerm(FunctionSpace<SVector3>& space1_, double E_, double nu_) :
- BilinearTerm<SVector3, SVector3>(space1_, space1_), E(E_), nu(nu_), H(6, 6)
- {
- double FACT = E / (1 + nu);
- double C11 = FACT * (1 - nu) / (1 - 2 * nu);
- double C12 = FACT * nu / (1 - 2 * nu);
- double C44 = (C11 - C12) / 2;
- FACT = E / (1 - nu * nu);
- C11 = FACT;
- C12 = nu * FACT;
- C44 = (1. - nu) * .5 * FACT;
- H.scale(0.);
- for(int i = 0; i < 3; ++i) { H(i, i) = C11; H(i + 3, i + 3) = C44; }
- H(1, 0) = H(0, 1) = H(2, 0) = H(0, 2) = H(1, 2) = H(2, 1) = C12;
- sym = true;
- }
+public :
+ IsotropicElasticTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<SVector3>& space2_, double E_, double nu_);
+ IsotropicElasticTerm(FunctionSpace<SVector3>& space1_, double E_, double nu_);
virtual ~IsotropicElasticTerm() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
- {
- if(ele->getParent()) ele = ele->getParent();
- if(sym){
- int nbFF = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele);
- double jac[3][3];
- fullMatrix<double> B(6, nbFF);
- fullMatrix<double> BTH(nbFF, 6);
- fullMatrix<double> BT(nbFF, 6);
- m.resize(nbFF, nbFF);
- m.setAll(0.);
- //std::cout << m.size1() << " " << m.size2() << std::endl;
- 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<TensorialTraits<SVector3>::GradType> Grads;
- BilinearTerm<SVector3, SVector3>::space1.gradf(ele, u, v, w, Grads); // a optimiser ??
- for(int j = 0; j < nbFF; j++){
- BT(j, 0) = B(0, j) = Grads[j](0, 0);
- BT(j, 1) = B(1, j) = Grads[j](1, 1);
- BT(j, 2) = B(2, j) = Grads[j](2, 2);
- BT(j, 3) = B(3, j) = Grads[j](0, 1) + Grads[j](1, 0);
- BT(j, 4) = B(4, j) = Grads[j](1, 2) + Grads[j](2, 1);
- BT(j, 5) = B(5, j) = Grads[j](0, 2) + Grads[j](2, 0);
- }
- BTH.setAll(0.);
- BTH.gemm(BT, H);
- m.gemm(BTH, B, weight * detJ, 1.);
- }
- }
- else{
- int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele);
- int nbFF2 = BilinearTerm<SVector3, SVector3>::space2.getNumKeys(ele);
- double jac[3][3];
- fullMatrix<double> B(6, nbFF2);
- fullMatrix<double> BTH(nbFF2, 6);
- fullMatrix<double> BT(nbFF1, 6);
- m.resize(nbFF1, nbFF2);
- m.setAll(0.);
- // Sum on Gauss Points i
- 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<TensorialTraits<SVector3>::GradType> Grads;// tableau de matrices...
- std::vector<TensorialTraits<SVector3>::GradType> GradsT;// tableau de matrices...
- BilinearTerm<SVector3, SVector3>::space1.gradf(ele, u, v, w, Grads);
- BilinearTerm<SVector3, SVector3>::space2.gradf(ele, u, v, w, GradsT);
- for(int j = 0; j < nbFF1; j++){
- BT(j, 0) = Grads[j](0, 0);
- BT(j, 1) = Grads[j](1, 1);
- BT(j, 2) = Grads[j](2, 2);
- BT(j, 3) = Grads[j](0, 1) + Grads[j](1, 0);
- BT(j, 4) = Grads[j](1, 2) + Grads[j](2, 1);
- BT(j, 5) = Grads[j](0, 2) + Grads[j](2, 0);
- }
- for(int j = 0; j < nbFF2; j++){
- B(0, j) = GradsT[j](0, 0);
- B(1, j) = GradsT[j](1, 1);
- B(2, j) = GradsT[j](2, 2);
- B(3, j) = GradsT[j](0, 1) + GradsT[j](1, 0);
- B(4, j) = GradsT[j](1, 2) + GradsT[j](2, 1);
- B(5, j) = GradsT[j](0, 2) + GradsT[j](2, 0);
- }
- BTH.setAll(0.);
- BTH.gemm(BT, H);
- // gemm add the product to m so there is a sum on gauss' points here
- m.gemm(BTH, B, weight * detJ, 1.);
- }
- }
- }
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< fullMatrix<double> > &vm) const{};
+ virtual BilinearTermBase* clone () const {return new IsotropicElasticTerm(BilinearTerm<SVector3, SVector3>::space1,BilinearTerm<SVector3, SVector3>::space2,E,nu);}
}; // class
-inline double dot(const double &a, const double &b)
-{
- return a * b;
-}
template<class T1> class LoadTerm : public LinearTerm<T1>
{
+protected:
simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
- public :
+public :
LoadTerm(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> &Load_) :
- LinearTerm<T1>(space1_), Load(Load_) {}
+ LinearTerm<T1>(space1_), Load(Load_) {}
virtual ~LoadTerm() {}
-
- virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m)
- {
- if(ele->getParent()) ele = ele->getParent();
- int nbFF = LinearTerm<T1>::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 TensorialTraits<T1>::ValType> Vals;
- LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
- SPoint3 p;
- ele->pnt(u, v, w, p);
- typename TensorialTraits<T1>::ValType load = Load(p.x(), p.y(), p.z());
- for(int j = 0; j < nbFF ; ++j){
- m(j) += dot(Vals[j], load) * weight * detJ;
- }
- }
- }
+ virtual LinearTermBase<double>* clone () const { return new LoadTerm<T1>(LinearTerm<T1>::space1,Load);}
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m) const ;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<double> > &vv) const {};
};
class LagrangeMultiplierTerm : public BilinearTerm<SVector3, double>
{
SVector3 _d;
- public :
+public :
LagrangeMultiplierTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<double>& space2_, const SVector3 &d) :
- BilinearTerm<SVector3, double>(space1_, space2_)
+ BilinearTerm<SVector3, double>(space1_, space2_)
{
- for(int i = 0; i < 3; i++) _d(i) = d(i);
+ for (int i = 0; i < 3; i++) _d(i) = d(i);
}
virtual ~LagrangeMultiplierTerm() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
- {
- int nbFF1 = BilinearTerm<SVector3, double>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
- int nbFF2 = BilinearTerm<SVector3, double>::space2.getNumKeys(ele); //nbVertices of boundary
- double jac[3][3];
- m.resize(nbFF1, nbFF2);
- m.setAll(0.);
- for(int i = 0; i < npts; i++){
- double u = GP[i].pt[0]; double v = GP[i].pt[1]; double w = GP[i].pt[2];
- const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
- std::vector<TensorialTraits<SVector3>::ValType> Vals;
- std::vector<TensorialTraits<double>::ValType> ValsT;
- BilinearTerm<SVector3,double>::space1.f(ele, u, v, w, Vals);
- BilinearTerm<SVector3,double>::space2.f(ele, u, v, w, ValsT);
- for(int j = 0; j < nbFF1; j++){
- for(int k = 0; k < nbFF2; k++){
- m(j, k) += dot(Vals[j], _d) * ValsT[k] * weight * detJ;
- }
- }
- }
- }
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< fullMatrix<double> > &vm) const{};
+ virtual BilinearTermBase* clone () const {return new LagrangeMultiplierTerm(BilinearTerm<SVector3, double>::space1,BilinearTerm<SVector3, double>::space2,_d);}
};
class LagMultTerm : public BilinearTerm<SVector3, SVector3>
{
- private :
+private :
double _eqfac;
- public :
+public :
LagMultTerm(FunctionSpace<SVector3>& space1_, FunctionSpace<SVector3>& space2_, double eqfac = 1.0) :
- BilinearTerm<SVector3, SVector3>(space1_, space2_), _eqfac(eqfac) {}
+ BilinearTerm<SVector3, SVector3>(space1_, space2_), _eqfac(eqfac) {}
virtual ~LagMultTerm() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m)
- {
- int nbFF1 = BilinearTerm<SVector3, SVector3>::space1.getNumKeys(ele); //nbVertices*nbcomp of parent
- int nbFF2 = BilinearTerm<SVector3, SVector3>::space2.getNumKeys(ele); //nbVertices of boundary
- double jac[3][3];
- m.resize(nbFF1, nbFF2);
- m.setAll(0.);
- for(int i = 0; i < npts; i++) {
- double u = GP[i].pt[0]; double v = GP[i].pt[1]; double w = GP[i].pt[2];
- const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
- std::vector<TensorialTraits<SVector3>::ValType> Vals;
- std::vector<TensorialTraits<SVector3>::ValType> ValsT;
- BilinearTerm<SVector3,SVector3>::space1.f(ele, u, v, w, Vals);
- BilinearTerm<SVector3,SVector3>::space2.f(ele, u, v, w, ValsT);
- for(int j = 0; j < nbFF1; j++){
- for(int k = 0; k < nbFF2; k++){
- m(j, k) += _eqfac * dot(Vals[j], ValsT[k]) * weight * detJ;
- }
- }
- }
- }
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector< fullMatrix<double> > &vm) const{};
+ virtual BilinearTermBase* clone () const {return new LagMultTerm(BilinearTerm<SVector3, SVector3>::space1,BilinearTerm<SVector3, SVector3>::space2,_eqfac);}
};
template<class T1> class LoadTermOnBorder : public LinearTerm<T1>
{
- private :
+private :
double _eqfac;
simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
- public :
+public :
LoadTermOnBorder(FunctionSpace<T1>& space1_, simpleFunction<typename TensorialTraits<T1>::ValType> &Load_, double eqfac = 1.0) :
- LinearTerm<T1>(space1_), _eqfac(eqfac), Load(Load_) {}
+ LinearTerm<T1>(space1_), _eqfac(eqfac), Load(Load_) {}
virtual ~LoadTermOnBorder() {}
- virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m)
- {
- MElement *elep;
- if (ele->getParent()) elep = ele->getParent();
- int nbFF = LinearTerm<T1>::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 TensorialTraits<T1>::ValType> Vals;
- LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
- SPoint3 p;
- ele->pnt(u, v, w, p);
- typename TensorialTraits<T1>::ValType load = Load(p.x(), p.y(), p.z());
- for(int j = 0; j < nbFF ; ++j){
- m(j) += _eqfac * dot(Vals[j], load) * weight * detJ;
- }
- }
- }
+ virtual LinearTermBase<double>* clone () const { return new LoadTermOnBorder<T1>(LinearTerm<T1>::space1,Load,_eqfac);}
+ virtual void get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m) const ;
+ virtual void get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<double> > &vv) const {};
};
+#include "terms.hpp"
+
#endif// _TERMS_H_
diff --git a/Solver/terms.hpp b/Solver/terms.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..844efa9d84761e7b32d35a13606498454ac1898e
--- /dev/null
+++ b/Solver/terms.hpp
@@ -0,0 +1,220 @@
+//
+// C++ Template Implementations: terms
+//
+// Description:
+//
+//
+// Author: <Eric Bechet>, (C) 2011
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "terms.h"
+
+template<class T2> void LinearTermBase<T2>::get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &vec) const
+{
+ std::vector<fullVector<T2> > vv;
+ vv.resize(npts);
+ get(ele,npts,GP,vv);
+ int nbFF=vv[0].size();
+ vec.resize(nbFF);
+ vec.setAll(T2());
+ double jac[3][3];
+ 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);
+ for(int j = 0; j < nbFF; j++)
+ {
+ double contrib = weight * detJ;
+ vec(j) += contrib*vv[i](j);
+ }
+ }
+}
+
+template<class T2> void PlusTerm<T2>::get(MElement *ele, int npts, IntPt *GP, fullVector<T2> &v) const
+{
+ fullVector<T2> v2;
+ a->get(ele,npts,GP,v);
+ b->get(ele,npts,GP,v2);
+ for (int i=0;i<v2.size();++i) v(i)+=v2(i);
+}
+
+template<class T2> void ScalarTermConstant<T2>::get(MElement *ele, int npts, IntPt *GP, T2 &val) const
+{
+ double jac[3][3];
+ double eval = 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);
+ eval += weight * detJ;
+ }
+ val=cst * eval;
+}
+
+template<class T2> void ScalarTermConstant<T2>::get(MElement *ele, int npts, IntPt *GP, std::vector<T2> &vval) const
+{
+ for(int i = 0; i < npts; i++)
+ {
+ vval[i] = cst;
+ }
+}
+
+template<class T2> void ScalarTermConstant<T2>::get(MVertex *ver, T2 &val) const
+{
+ val = cst;
+}
+
+template<class T1> void LaplaceTerm<T1, T1>::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+ int nbFF = BilinearTerm<T1, T1>::space1.getNumKeys(ele);
+ double jac[3][3];
+ m.resize(nbFF, nbFF);
+ m.setAll(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 TensorialTraits<T1>::GradType> Grads;
+ BilinearTerm<T1, T1>::space1.gradf(ele, u, v, w, Grads);
+ for(int j = 0; j < nbFF; j++)
+ {
+ for(int k = j; k < nbFF; k++)
+ {
+ double contrib = weight * detJ * dot(Grads[j], Grads[k]) * diffusivity;
+ m(j, k) += contrib;
+ if(j != k) m(k, j) += contrib;
+ }
+ }
+ }
+}
+
+template<class T1> void LoadTerm<T1>::get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m) const
+{
+ if(ele->getParent()) ele = ele->getParent();
+ int nbFF = LinearTerm<T1>::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 TensorialTraits<T1>::ValType> Vals;
+ LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
+ SPoint3 p;
+ ele->pnt(u, v, w, p);
+ typename TensorialTraits<T1>::ValType load = Load(p.x(), p.y(), p.z());
+ for(int j = 0; j < nbFF ; ++j)
+ {
+ m(j) += dot(Vals[j], load) * weight * detJ;
+ }
+ }
+}
+
+template<class T2> BilinearTermContract<T2> operator |(const LinearTermBase<T2>& L1,const LinearTermBase<T2>& L2)
+{
+ return BilinearTermContract<T2>(L1,L2);
+}
+
+template<class T2> void BilinearTermContract<T2>::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+ fullVector<T2> va;
+ fullVector<T2> vb;
+ a->get(ele,npts,GP,va);
+ b->get(ele,npts,GP,vb);
+ m.resize(va.size(), vb.size());
+ m.setAll(0.);
+ for (int i=0;i<va.size();++i)
+ for (int j=0;j<vb.size();++j)
+ m(i,j)=dot(va(i),vb(j));
+}
+
+
+template<class T2> void BilinearTermContractWithLaw<T2>::get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const
+{
+ BilinearTermBase::get(ele,npts,GP,m);
+}
+
+template<class T2> void BilinearTermContractWithLaw<T2>::get(MElement *ele, int npts, IntPt *GP, std::vector<fullMatrix<double> > &mv) const
+{
+ std::vector<fullVector<T2> > va(npts);
+ std::vector<fullVector<T2> > vb(npts);
+ std::vector<typename TensorialTraits<T2>::TensProdType> tens(npts);
+ BilinearTermContract<T2>::a->get(ele,npts,GP,va);
+ BilinearTermContract<T2>::b->get(ele,npts,GP,vb);
+ c->get(ele,npts,GP,tens);
+ for (int k=0;k<npts;k++)
+ {
+ mv[k].resize(va[k].size(), vb[k].size());
+ for (int i=0;i<va[k].size();++i)
+ for (int j=0;j<vb[k].size();++j)
+ mv[k](i,j)=dot(va[k](i),tens[k]*vb[k](j));
+ }
+}
+
+template<class T2> PlusTerm<T2> LinearTermBase<T2>::operator +(const LinearTermBase<T2>& other)
+{
+ return PlusTerm<T2>(*this,other);
+}
+/*
+template<class T1> void GradTerm<T1>::get(MElement *ele, int npts, IntPt *GP, fullVector<typename TensorialTraits<T1>::GradType > &vec) const
+{
+ int nbFF = LinearTerm<T1,typename TensorialTraits<T1>::GradType>::space1.getNumKeys(ele);
+ double jac[3][3];
+ vec.resize(nbFF);
+ vec.setAll(typename TensorialTraits<T1>::GradType());
+ for(int i = 0; i < npts; i++)
+ {
+ std::vector<typename TensorialTraits<T1>::GradType> Grads;
+ 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);
+ LinearTerm<T1,typename TensorialTraits<T1>::GradType>::space1.gradf(ele, u, v, w, Grads);
+ for(int j = 0; j < nbFF; j++)
+ {
+ double contrib = weight * detJ;
+ vec(j) += contrib*Grads[j];
+ }
+ }
+}
+*/
+template<class T1> void GradTerm<T1>::get(MElement *ele, int npts, IntPt *GP, std::vector<fullVector<typename TensorialTraits<T1>::GradType > > &vvec) const
+{
+ int nbFF = LinearTerm<T1,typename TensorialTraits<T1>::GradType>::space1.getNumKeys(ele);
+ for(int i = 0; i < npts; i++)
+ {
+ vvec[i].resize(nbFF);
+ std::vector<typename TensorialTraits<T1>::GradType> Grads;
+ const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
+ LinearTerm<T1,typename TensorialTraits<T1>::GradType>::space1.gradf(ele, u, v, w, Grads);
+ for(int j = 0; j < nbFF; j++)
+ {
+ vvec[i](j) = Grads[j];
+ }
+ }
+}
+
+
+
+template<class T1> void LoadTermOnBorder<T1>::get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m) const
+{
+ MElement *elep;
+ if (ele->getParent()) elep = ele->getParent();
+ int nbFF = LinearTerm<T1>::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 TensorialTraits<T1>::ValType> Vals;
+ LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
+ SPoint3 p;
+ ele->pnt(u, v, w, p);
+ typename TensorialTraits<T1>::ValType load = Load(p.x(), p.y(), p.z());
+ for(int j = 0; j < nbFF ; ++j){
+ m(j) += _eqfac * dot(Vals[j], load) * weight * detJ;
+ }
+ }
+}