pp

contrib/cm3/CMakeLists.txt~

-# pour boris:
-## boris cela ne marche que si le fichier est dans SVN ...
-#include(/home/boris/MyGmshProjects/CMakeLists.txt)
-# tests pour eric et christophe:
-list(APPEND EXTERNAL_INCLUDES contrib/cm3)
-add_executable(dgsolver EXCLUDE_FROM_ALL contrib/cm3/mainDG.cpp contrib/cm3/DgC0PlateSolver.cpp contrib/cm3/MInterfaceElement.cpp contrib/cm3/GModelDg.cpp contrib/cm3/C0DgPlateTerms.h ${GMSH_SRC})
-target_link_libraries(dgsolver ${LINK_LIBRARIES})
-

contrib/cm3/DGterms.h~

-//
-// C++ Interface: terms
-//
-// Description: 
-//
-//
-// Author:  <Eric Bechet>, (C) 2009
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-
-#ifndef _TERMS_H_
-#define _TERMS_H_
-
-#include "SVector3.h"
-#include <vector>
-#include <iterator>
-#include "Numeric.h"
-#include "functionSpace.h"
-#include "groupOfElements.h"
-#include "materialLaw.h"
-
-class  BilinearTermBase
-{
- public :
-  virtual ~BilinearTermBase() {}
-  virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
-};
-
-template<class T1,class T2> class BilinearTerm : public BilinearTermBase
-{
- protected :
-  FunctionSpace<T1>& space1;
-  FunctionSpace<T2>& space2;
- public :
-  BilinearTerm(FunctionSpace<T1>& space1_,FunctionSpace<T1>& space2_) : space1(space1_),space2(space2_) {}
-  virtual ~BilinearTerm() {}
-};
-
-class  LinearTermBase
-{
-  public:
-  virtual ~LinearTermBase() {}
-  virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0;
-  virtual void get(MVertex *ver,fullVector<double> &m) =0;
-};
-
-template<class T1> class LinearTerm : public LinearTermBase
-{
- protected :
-  FunctionSpace<T1>& space1;
- 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
-{
-  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);
-  }
-};
-
-
-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 : 
-  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)
-  {
-    Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
-  }
-  virtual void get(MVertex *ver,fullMatrix<double> &m)
-  {
-    Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
-  }
-}; // class
-
-
-
-template<class T1> class LaplaceTerm<T1,T1> : public BilinearTerm<T1,T1> // symmetric
-{
- public : 
-  LaplaceTerm(FunctionSpace<T1>& space1_) : BilinearTerm<T1,T1>(space1_,space1_)
-  {}
-  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]);
-          m(j,k)+=contrib;
-          if (j!=k) m(k,j)+=contrib;
-        }
-      }
-    }
-//    m.print("");
-//    exit(0);
-  }
-}; // class
-
-
-
-
-class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3> 
-{
- protected : 
-  double E,nu;
-  bool sym;
-  fullMatrix<double> H;/* =
-    { {C11, C12, C12,    0,   0,   0},
-      {C12, C11, C12,    0,   0,   0},
-      {C12, C12, C11,    0,   0,   0},
-      {  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;
-    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;
-  }
-  virtual ~IsotropicElasticTerm() {}
-  virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
-  {
-    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.);
-      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.);
-      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);
-        m.gemm(BTH, B, weight * detJ, 1.);
-      }
-    }
-  }
-}; // class
-
-
-inline double dot(const double &a, const double &b)
-{ return a*b; }
-
-
-template<class T1> class LoadTerm : public LinearTerm<T1>
-{
-  simpleFunction<typename TensorialTraits<T1>::ValType> &Load;
- public : 
-  LoadTerm(FunctionSpace<T1>& space1_,simpleFunction<typename TensorialTraits<T1>::ValType> &Load_) :LinearTerm<T1>(space1_),Load(Load_) {}
-  virtual ~LoadTerm() {}
-
-  virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
-  {
-    double 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 void get(MVertex *ver,fullVector<double> &m)
-  {
-    double nbFF=LinearTerm<T1>::space1.getNumKeys(ver);
-    double jac[3][3];
-    m.resize(nbFF);
-    std::vector<typename TensorialTraits<T1>::ValType> Vals;
-    LinearTerm<T1>::space1.f(ver, Vals);
-    typename TensorialTraits<T1>::ValType load=Load(ver->x(),ver->y(),ver->z());
-    for (int j = 0; j < nbFF ; ++j)
-    {
-      m(j)=dot(Vals[j],load);
-    }
-  }
-};
-
-
-
-
-#endif// _TERMS_H_

contrib/cm3/DgC0PlateElementaryTerms1ddl.h~

-// Compute the component (j,k) of the elementary stiffness matrix
-inline double BulkC0PlateDGStiffnessBendingTerms(TensorialTraits<double>::HessType &hessj, TensorialTraits<double>::HessType &hessk, LinearElasticShellHookeTensor *H, const LocalBasis *lb){
-  double val = 0.;
-  for(int alpha=0;alpha<2;alpha++)
-    for(int beta=0;beta<2;beta++)
-      for(int gamma=0;gamma<2;gamma++)
-        for(int delta=0;delta<2;delta++)
-          val += hessj(alpha,beta)*hessk(gamma,delta)*H->get(alpha,beta,gamma,delta);
-  return lb->gett0(2)*lb->gett0(2)*val;
-}
-
-inline double consAndCompC0PlateStiffnessTerms(LinearElasticShellHookeTensor *Hhat,const fullMatrix<double> &Bhat,const SVector3 &dt, const LocalBasis *lb){
-  double val=0.;
-  for(int alpha=0;alpha<2;alpha++)
-    for(int beta=0;beta<2;beta++)
-      for(int gamma=0;gamma<2;gamma++)
-        for(int delta=0;delta<2;delta++)
-          val += Hhat->get(alpha,beta,gamma,delta)*Bhat(gamma,delta)*dot(dt,lb->getphi0(alpha))*(-lb->getphi0(1,beta));
-  return 0.5*val;
-}
-
-inline double stabilityC0PlateStiffnessTerms(LinearElasticShellHookeTensor *Hhat, const SVector3 &dta, const SVector3 &dtb, const LocalBasis *lb){
-  double val=0.;
-  for(int alpha=0;alpha<2;alpha++)
-    for(int beta=0;beta<2;beta++)
-      for(int gamma=0;gamma<2;gamma++)
-        for(int delta=0;delta<2;delta++)
-          val += dot(dta,lb->getphi0(gamma))*dot(dtb,lb->getphi0(alpha))*Hhat->get(alpha,beta,gamma,delta)*(-lb->getphi0(1,delta))*(-lb->getphi0(1,beta));
-
-  return val;
-}
-
-inline double BulkC0PlateDGForceTerms(const TensorialTraits<double>::HessType &hessj,const std::vector<TensorialTraits<double>::HessType> &Hess,const LinearElasticShellHookeTensor *H,const LocalBasis *lb,const fullMatrix<double> &disp){
-  const int n = Hess.size();
-  double sum,val;
-  val=0.;
-  for(int a=0;a<2;a++)
-    for(int b=0;b<2;b++){
-      sum=0.;
-      for(int j=0;j<n;j++)
-        sum+=-lb->gett0(2)*Hess[j](a,b)*disp(j,0);
-      for(int c=0;c<2;c++)
-        for(int d=0;d<2;d++)
-          val += -lb->gett0(2)*hessj(c,d)*H->get(a,b,c,d)*sum;
-    }
-  return val;
-}
-
-inline double consC0PlateForceTerms(const LinearElasticShellHookeTensor *Hhat, const fullMatrix<double> &Bhat, const std::vector<SVector3> &Dt_m, const std::vector<SVector3> &Dt_p,const LocalBasis *lb, const fullMatrix<double> &disp){
-  const int n_m = Dt_m.size();
-  const int n_p = Dt_p.size();
-  double sum,val;
-  val=0.;
-  for(int a=0;a<2;a++){
-    sum=0.;
-    for(int j=0;j<n_m;j++)
-      sum -= dot(Dt_m[j],lb->getphi0(a))*disp(j,0);
-    for(int j=0;j<n_p;j++)
-      sum += dot(Dt_p[j],lb->getphi0(a))*disp(j+n_m,0);
-
-    for(int b=0;b<2;b++)
-      for(int c=0;c<2;c++)
-        for(int d=0;d<2;d++)
-          val+= Hhat->get(a,b,c,d)*Bhat(c,d)*sum*(-lb->getphi0(1,b));
-  }
-  return 0.5*val;
-}
-
-

contrib/cm3/mainElasticity.cpp~

-#include "Gmsh.h"
-#include "elasticitySolver.h"
-#include "PView.h"
-#include "PViewData.h"
-#include "highlevel.h"/usr/bin/gmake -f "/home/gauthier/distgmsh/trunk/Makefile" elastic -o ../../bin/elastic
-#include "groupOfElements.h"
-#include <iterator>
-
-int main (int argc, char* argv[])
-{
-
-  if (argc != 2){
-    printf("usage : elasticity input_file_name\n");
-    return -1;
-  }
-
-
-  GmshInitialize(argc, argv);
-  // globals are still present in Gmsh
-
-  // instanciate a solver
-  elasticitySolver mySolver (1000);
-
-  // read some input file
-  mySolver.readInputFile(argv[1]);
-
-  // solve the problem
-  mySolver.solve();
-
-  PView *pv = mySolver.buildDisplacementView("displacement");
-  pv->getData()->writeMSH("disp.msh", false);
-  delete pv;
-  pv = mySolver.buildElasticEnergyView("elastic energy");
-  pv->getData()->writeMSH("energ.msh", false);
-  delete pv;
-
-
-
-  // stop gmsh
-  GmshFinalize();
-
-}
-
-
-
-/*
-  groupOfElements *g = new groupOfElements (2, 7);
-
-  MElement *e=*(g->begin());
-  std::cout << e->getNumPrimaryVertices() << "vertices" << std::endl;
-  const double uvw[3]={0.,0.,0.};
-  std::vector<Dof> dofs;
-  std::vector<double> vals;
-  std::vector<SVector3> grads;
-  std::vector<SVector3> vals2;
-  std::vector<STensor3> grads2;
-
-  std::ostream_iterator< double > output( std::cout, " " );
-
-  ScalarLagrangeFunctionSpace L(100);
-  std::cout << L.getNumKeys(e) << "fonctions de formes L" << std::endl;
-  L.getKeys(e,dofs);
-  for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl;
-  dofs.clear();
-  L.f(e,0.1,0.1,0,vals);
-  L.gradf(e,0.1,0.1,0,grads);
-  std::copy(vals.begin(),vals.end(),output); std::cout << std::endl;
-  for (std::vector<SVector3>::iterator it=grads.begin();it!=grads.end();++it) { std::cout << (*it)[0]<< " " << (*it)[1] <<" " << (*it)[2] << std::endl; }
-
-  VectorLagrangeFunctionSpace L1(100,VectorLagrangeFunctionSpace::VECTOR_X);
-  VectorLagrangeFunctionSpace L2(100,VectorLagrangeFunctionSpace::VECTOR_Y);
-  std::cout << L2.getNumKeys(e) << "fonctions de formes L2" << std::endl;
-  L2.f(e,0.1,0.1,0,vals2);
-  L2.gradf(e,0.1,0.1,0,grads2);
-  for (std::vector<SVector3>::iterator it=vals2.begin();it!=vals2.end();++it) { std::cout << (*it)[0]<< " " << (*it)[1] <<" " << (*it)[2] << std::endl; }
-  for (std::vector<STensor3>::iterator it=grads2.begin();it!=grads2.end();++it) { (*it).print(""); }
-
-  VectorLagrangeFunctionSpace L3(100,VectorLagrangeFunctionSpace::VECTOR_Z);
-
-  VectorLagrangeFunctionSpace P123(100);
-  std::cout << P123.getNumKeys(e) << "fonctions de formes P123" << std::endl;
-  P123.getKeys(e,dofs);
-  std::cout << dofs.size() << std::endl;
-  for (int i=0;i<dofs.size();++i) std::cout << "entity: " << dofs[i].getEntity() << " id: " << dofs[i].getType() << std::endl;
-
-  vals2.clear();
-  grads2.clear();
-  P123.f(e,0.1,0.1,0,vals2);
-  P123.gradf(e,0.1,0.1,0,grads2);
-  for (std::vector<SVector3>::iterator it=vals2.begin();it!=vals2.end();++it) { std::cout << (*it)[0]<< " " << (*it)[1] <<" " << (*it)[2] << std::endl; }
-  for (std::vector<STensor3>::iterator it=grads2.begin();it!=grads2.end();++it) { (*it).print(""); }
-
-
-
-  FormBilinear<TermBilinearMeca,ScalarLagrangeFunctionSpace,ScalarLagrangeFunctionSpace > f(L,L);
-  f.func();
-  f.Accumulate(e,uvw);
-
-
-  FormBilinear<TermBilinearMecaNL,ScalarLagrangeFunctionSpace,ScalarLagrangeFunctionSpace > fnl(L,L);
-  fnl.func();
-
-*/
-
-
-
-
-