From f03912b585871069e6b9fb8124dae34fb82b07f0 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 10 Mar 2010 17:56:28 +0000
Subject: [PATCH] pp
---
contrib/cm3/CMakeLists.txt~ | 8 -
contrib/cm3/DGterms.h~ | 331 --------------------
contrib/cm3/DgC0PlateElementaryTerms1ddl.h~ | 69 ----
contrib/cm3/mainElasticity.cpp~ | 108 -------
4 files changed, 516 deletions(-)
delete mode 100644 contrib/cm3/CMakeLists.txt~
delete mode 100644 contrib/cm3/DGterms.h~
delete mode 100644 contrib/cm3/DgC0PlateElementaryTerms1ddl.h~
delete mode 100644 contrib/cm3/mainElasticity.cpp~
diff --git a/contrib/cm3/CMakeLists.txt~ b/contrib/cm3/CMakeLists.txt~
deleted file mode 100644
index ecca2e400f..0000000000
--- a/contrib/cm3/CMakeLists.txt~
+++ /dev/null
@@ -1,8 +0,0 @@
-# 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})
-
diff --git a/contrib/cm3/DGterms.h~ b/contrib/cm3/DGterms.h~
deleted file mode 100644
index de33311651..0000000000
--- a/contrib/cm3/DGterms.h~
+++ /dev/null
@@ -1,331 +0,0 @@
-//
-// 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_
diff --git a/contrib/cm3/DgC0PlateElementaryTerms1ddl.h~ b/contrib/cm3/DgC0PlateElementaryTerms1ddl.h~
deleted file mode 100644
index 1fa6c24644..0000000000
--- a/contrib/cm3/DgC0PlateElementaryTerms1ddl.h~
+++ /dev/null
@@ -1,69 +0,0 @@
-// 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;
-}
-
-
diff --git a/contrib/cm3/mainElasticity.cpp~ b/contrib/cm3/mainElasticity.cpp~
deleted file mode 100644
index 3ac287aedc..0000000000
--- a/contrib/cm3/mainElasticity.cpp~
+++ /dev/null
@@ -1,108 +0,0 @@
-#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();
-
-*/
-
-
-
-
-
--
GitLab