Fix bug in NLS: creation of interface ans interdomain, Now no Dof type

suposition in solver
Update all benchmark lua --> python

NonLinearSolver/BoundaryConditions/nonLinearBC.h

@@ -56,7 +56,7 @@ class nonLinearDirichletBC : public nonLinearBoundaryCondition
 class rigidContactBC : public nonLinearBoundaryCondition
 {
  public:
-  rigidContactSpace *space;
+  rigidContactSpaceBase *space;
   int _comp; // component
   simpleFunctionTime<double> _f;
   rigidContactBC(const int physMaster) : nonLinearBoundaryCondition(){
@@ -70,7 +70,7 @@ class rigidContactBC : public nonLinearBoundaryCondition
     _f = source._f;
   }
   ~rigidContactBC(){}
-  void setSpace(rigidContactSpace *sp){space = sp;}
+  void setSpace(rigidContactSpaceBase *sp){space = sp;}
 };
 #endif
 class nonLinearNeumannBC  : public nonLinearBoundaryCondition

NonLinearSolver/Domain/partDomain.h

@@ -116,10 +116,6 @@ class partDomain{
 
 //  static void registerBindings(binding *b);
   virtual void computeIPVariable(AllIPState *aips,const unknownField *ufield,const IPStateBase::whichState ws)=0;
-  virtual void computeIpv(AllIPState *aips,MInterfaceElement *ie, IntPt *GP,const IPStateBase::whichState ws,
-                           partDomain* efMinus, partDomain *efPlus,materialLaw *mlawminus,
-                           materialLaw *mlawplus,fullVector<double> &dispm,
-                           fullVector<double> &dispp,const bool virt,const bool checkfrac=true)=0;
   virtual void computeIpv(AllIPState *aips,MElement *e, IPStateBase::whichState ws,
                            materialLaw *mlaw,fullVector<double> &disp)=0;
   virtual void setGaussIntegrationRule()=0;

NonLinearSolver/Interface/MInterfaceElement.h

@@ -24,7 +24,7 @@ class MInterfaceElement{
   virtual int getEdgeOrFaceNumber(const int index) const=0;
   virtual bool isSameDirection(const int index) const=0;
   // should return the element number !!
-  virtual int getNum() const=0;
+  virtual int getNumber() const=0; // {return{this->getNum();} in your derived class it derived from an MElement* too !!
   // compute the characteritic size of one element (This function can be defined as a method of MElement) ??
   static double characSize(MElement *e); // Area/perimeter
 };

NonLinearSolver/Interface/MInterfaceLine.cpp

@@ -13,8 +13,7 @@
 
 #include "MInterfaceLine.h"
 MInterfaceLine::MInterfaceLine(std::vector<MVertex*> &v, int num, int part,
-                               MElement *e_minus, MElement *e_plus) : MLineN(v, num, part), MInterfaceElement(),
-                                                                      _num(MElement::getGlobalNumber()) //avoid this ??
+                               MElement *e_minus, MElement *e_plus) : MLineN(v, num, part), MInterfaceElement()
 {
   _numElem[0]=e_minus;
   _numElem[1]=e_plus;

NonLinearSolver/Interface/MInterfaceLine.h

@@ -18,8 +18,6 @@
 #include "MInterfaceElement.h"
 class MInterfaceLine : public MLineN, public MInterfaceElement{ // or don't derivate but in this case a vector with the vertices of interface element has to be save ??
   protected :
-    // NUMBER DUPLICATION HOW TO AVOID THIS (MInterfaceElement : public MELelement but 2 MElement at this time and MELement has to be a virtual class ??)
-    int _num;
     // table of pointer on the two elements linked to the interface element
     MElement *_numElem[2];
     // edge's number linked to interface element of minus and plus element
@@ -35,8 +33,7 @@ class MInterfaceLine : public MLineN, public MInterfaceElement{ // or don't deri
     ~MInterfaceLine(){}
 
     // Try to avoid this HOW
-//    int getNum() const{MLineN::getNum();}
-    int getNum() const{return _num;}
+    int getNumber() const{return this->getNum();}
     // Give the number of minus 0 or plus 1 element
     MElement* getElem(int index) const {return _numElem[index];}
 

NonLinearSolver/contact/contactDomain.cpp

@@ -47,11 +47,11 @@ MElement* contactDomain::getFirstElement() const {
   return (*it);
 }
 
-rigidCylinderContact::rigidCylinderContact(const int tag, const int physMaster, const int physSlave, const int physPt1,
+rigidCylinderContactDomain::rigidCylinderContactDomain(const int tag, const int physMaster, const int physSlave, const int physPt1,
                                            const int physPt2, const double penalty,
-                                           const double h) : contactDomain(tag,physMaster,physSlave,
+                                           const double h,const double rho) : contactDomain(tag,physMaster,physSlave,
                                                                            penalty,rigidCylinder,true),
-                                                               _thick(h){
+                                                               _thick(h), _density(rho){
 
   // void gauss integration
   _integ = new QuadratureVoid();
@@ -98,18 +98,10 @@ rigidCylinderContact::rigidCylinderContact(const int tag, const int physMaster,
   _axisDirector->normalize();
 }
 
-void rigidCylinderContact::initializeTerms(const unknownField *ufield){
-  rigidContactSpace *sp = dynamic_cast<rigidContactSpace*>(_space);
+void rigidCylinderContactDomain::initializeTerms(const unknownField *ufield){
+  rigidContactSpaceBase *sp = dynamic_cast<rigidContactSpaceBase*>(_space);
   _massterm = new massRigidCylinder(this,sp);
   _lterm = new forceRigidCylinderContact(this,sp,_thickContact,ufield);
   rigidContactLinearTermBase<double> *rlterm = dynamic_cast<rigidContactLinearTermBase<double>*>(_lterm);
   _bterm = new stiffnessRigidCylinderContact(sp,rlterm,ufield);
 }
-
-void rigidCylinderContact::setDomain(partDomain *dom){
-  _dom = dom;
-  FunctionSpace<double> *spdom = dynamic_cast<FunctionSpace<double>*>(dom->getFunctionSpace());
-  if(_space != NULL) delete _space;
-  _space = new rigidContactSpace(_phys,spdom,_vergc);
-}
-

NonLinearSolver/contact/contactDomain.h

@@ -52,7 +52,6 @@ class contactDomain{
   virtual void setTag(const int t){ _tag =t;}
   virtual void setPhys(const int p){_phys =p;}
   virtual void setPenalty(const double pe){_penalty=pe;}
-  virtual void setDomain(partDomain *dom)=0;
   virtual void setContactType(const int ct);
   virtual void setContactType(const contact ct){_contype =ct;}
   virtual int getTag() const{return _tag;}
@@ -69,11 +68,13 @@ class contactDomain{
   virtual FunctionSpaceBase* getSpace() {return _space;}
   virtual const FunctionSpaceBase* getSpace() const {return _space;}
   virtual QuadratureBase* getGaussIntegration() const{return _integ;}
+  virtual void setDomainAndFunctionSpace(partDomain *dom)=0;
+  virtual void initializeTerms(const unknownField *ufield)=0;
 //  static void registerBindings(binding *b);
 #endif
 };
 
-class rigidCylinderContact : public contactDomain{
+class rigidCylinderContactDomain : public contactDomain{
  protected:
   double _length; // length of cylinder;
   double _radius; // outer radius of cylinder;
@@ -82,27 +83,24 @@ class rigidCylinderContact : public contactDomain{
   double _density; // density of cylinder Not a material law for now
   MVertex *_vergc;  // vertex of gravity center
   SVector3 *_axisDirector; // normalized vector director of cylinder's axis
+#ifndef SWIG
  public:
-  rigidCylinderContact() : contactDomain(), _length(0.), _radius(0.), _thick(0.){
+  rigidCylinderContactDomain() : contactDomain(), _length(0.), _radius(0.), _thick(0.){
     _contype = rigidCylinder;
     _rigid = true;
     _integ = new QuadratureVoid();
   }
-  rigidCylinderContact(const int tag, const int physMaster, const int physSlave, const int physpt1,
-                       const int physpt2,const double penalty,const double h);
-  ~rigidCylinderContact(){delete _axisDirector;}
-  void density(const double rho){ _density = rho;}
-//  void setThicknessContact(const double h){_thickContact = h;}
-#ifndef SWIG
+  rigidCylinderContactDomain(const int tag, const int physMaster, const int physSlave, const int physpt1,
+                       const int physpt2,const double penalty,const double h,const double rho);
+  ~rigidCylinderContactDomain(){delete _axisDirector;}
+  virtual void initializeTerms(const unknownField *ufield);
+  SVector3* getAxisDirector() const{return _axisDirector;}
+  virtual void setDomainAndFunctionSpace(partDomain *dom)=0;
   virtual MVertex* getGC() const{return _vergc;}
   double getDensity() const{return _density;}
   double getLength() const{return _length;}
   double getRadius() const{return _radius;}
   double getThickness() const{return _thick;}
-  void initializeTerms(const unknownField *ufield);
-  SVector3* getAxisDirector() const{return _axisDirector;}
-  virtual void setDomain(partDomain *dom);
 #endif
-
 };
 #endif // CONTACTDOMAIN_H_

NonLinearSolver/contact/contactFunctionSpace.h

@@ -12,7 +12,6 @@
 #ifndef CONTACTSPACE_H_
 #define CONTACTSPACE_H_
 #include "functionSpace.h"
-#include "Dof3IntType.h"
 template <class T1> class ContactSpaceBase : public FunctionSpaceBase{
  protected:
   FunctionSpace<T1> *_space; // functional space of domain where contact is computed
@@ -42,44 +41,19 @@ template <class T1> class ContactSpaceBase : public FunctionSpaceBase{
   }
 };
 
-class rigidContactSpace : public ContactSpaceBase<double>{
+class rigidContactSpaceBase : public ContactSpaceBase<double>{
  protected:
   const MVertex *_gc; // Node of gravity center
  public:
-  rigidContactSpace(const int id, FunctionSpace<double> *sp, MVertex *ver) : ContactSpaceBase<double>(id,sp), _gc(ver){}
-  rigidContactSpace(const int id, FunctionSpace<double> *sp, MVertex *ver, const int comp1,
+  rigidContactSpaceBase(const int id, FunctionSpace<double> *sp, MVertex *ver) : ContactSpaceBase<double>(id,sp), _gc(ver){}
+  rigidContactSpaceBase(const int id, FunctionSpace<double> *sp, MVertex *ver, const int comp1,
                    const int comp2 = -1, const int comp3 =-1) : ContactSpaceBase<double>(id,sp,comp1,comp2,comp3){}
-  virtual void getKeys(MElement *ele,std::vector<Dof> &keys){
-    _space->getKeys(ele,keys);
-    for(int i=0;i<_comp.size();i++){
-      keys.push_back(Dof(_gc->getNum(),Dof3IntType::createTypeWithThreeInts(_comp[i],_id)));
-    }
-  }
-  virtual int getNumKeys(MElement *ele){
-    int nkeysele = _space->getNumKeys(ele);
-    return nkeysele + _comp.size();
-  }
-  virtual void getKeysOfGravityCenter(std::vector<Dof> &keys){
-    for(int i=0;i<_comp.size(); i++)
-      keys.push_back(Dof(_gc->getNum(),Dof3IntType::createTypeWithThreeInts(_comp[i],_id)));
-  }
-  virtual int getNumKeysOfGravityCenter(){
-    return _comp.size();
-  }
-  virtual void getKeys(MElement *ele, const int ind, std::vector<Dof> &keys){
-    // generate keys of element and select the good ones after LAGRANGE OK ?? CHANGE THIS HOW TODO
-    std::vector<Dof> tkeys;
-    this->getKeys(ele,tkeys);
-    int nkeys = this->getNumKeys(ele);
-    int ncomp = _comp.size();
-    int nbver = nkeys/ncomp-1; // because GC is in the nkeys
-    for(int i=0;i<ncomp; i++)
-      keys.push_back(tkeys[ind+i*nbver]);
-    this->getKeysOfGravityCenter(keys);
-  }
-  virtual int getNumKeys(MElement *ele, int ind){
-    return 2*_comp.size();
-  }
+  virtual void getKeys(MElement *ele,std::vector<Dof> &keys)=0;
+  virtual int getNumKeys(MElement *ele)=0;
+  virtual void getKeysOfGravityCenter(std::vector<Dof> &keys)=0;
+  virtual int getNumKeysOfGravityCenter()=0;
+  virtual void getKeys(MElement *ele, const int ind, std::vector<Dof> &keys)=0;
+  virtual int getNumKeys(MElement *ele, int ind)=0; // number key in one component + number of key of GC
 };
 #endif // CONTACTSPACE_H_
 

NonLinearSolver/contact/contactTerms.h

@@ -50,7 +50,7 @@ template <class T2=double>class  contactBilinearTermBase : public BilinearTermBa
 template<class T2=double> class rigidContactLinearTermBase : public contactLinearTermBase<T2>{
  protected:
   const unknownField *_ufield;
-  rigidContactSpace *_spc;
+  rigidContactSpaceBase *_spc;
   // data for get function (Allocated once)
   mutable SVector3 _lastNormalContact; // allow to use the normal in the 2 get functions not very elegant
   mutable int nbdofs,nbdofsGC, nbvertex,nbcomp,nbcomptimenbvertex;
@@ -60,7 +60,7 @@ template<class T2=double> class rigidContactLinearTermBase : public contactLinea
   mutable double mvalue[6];
   mutable MVertex *ver;
  public:
-  rigidContactLinearTermBase(rigidContactSpace *sp,
+  rigidContactLinearTermBase(rigidContactSpaceBase *sp,
                              const unknownField *ufield) : contactLinearTermBase<T2>(),
                                                                  _ufield(ufield), _spc(sp), _lastNormalContact(0.,0.,0.),
                                                                  nbdofs(0), nbdofsGC(0), nbvertex(0), nbcomp(0), nbcomptimenbvertex(0)

NonLinearSolver/contact/rigidCylinderContactTerms.cpp

@@ -11,8 +11,8 @@
 #include "rigidCylinderContactTerms.h"
 #include "contactDomain.h"
 #include "unknownField.h"
-massRigidCylinder::massRigidCylinder(rigidCylinderContact *cdom,
-                                     rigidContactSpace *spc) : _length(cdom->getLength()),_radius(cdom->getRadius()),
+massRigidCylinder::massRigidCylinder(rigidCylinderContactDomain *cdom,
+                                     rigidContactSpaceBase *spc) : _length(cdom->getLength()),_radius(cdom->getRadius()),
                                                                _thick(cdom->getThickness()),_rho(cdom->getDensity()),
                                                                _spc(spc),nbdofs(0.), masse(0.),Rint(0.){}
 
@@ -32,8 +32,8 @@ void massRigidCylinder::get(MElement *ele,int npts,IntPt *GP,fullMatrix<double>
       m(i,i) = masse;
 }
 
- forceRigidCylinderContact::forceRigidCylinderContact(const rigidCylinderContact *cdom,
-                            rigidContactSpace *sp,const double thick,
+ forceRigidCylinderContact::forceRigidCylinderContact(const rigidCylinderContactDomain *cdom,
+                            rigidContactSpaceBase *sp,const double thick,
                             const unknownField *ufield) : rigidContactLinearTermBase(sp,ufield), _cdom(cdom),
                                                                 _vergc(cdom->getGC()),
                                                                 _axisDirector(cdom->getAxisDirector()),

NonLinearSolver/contact/rigidCylinderContactTerms.h

@@ -13,18 +13,18 @@
 #define RIGIDCONTACTCYLINDER_H_
 #include "contactTerms.h"
 #include "contactFunctionSpace.h"
-class rigidCylinderContact;
+class rigidCylinderContactDomain;
 class massRigidCylinder : public BilinearTermBase{
  protected:
   double _length, _radius, _thick, _rho;
-  rigidContactSpace *_spc;
+  rigidContactSpaceBase *_spc;
   // Data for get function (Allocated once)
   mutable int nbdofs;
   mutable double masse;
   mutable double Rint;
  public:
-  massRigidCylinder(rigidCylinderContact *cdom,rigidContactSpace *spc);
-  massRigidCylinder(rigidContactSpace *spc,double leng, double radius,
+  massRigidCylinder(rigidCylinderContactDomain *cdom,rigidContactSpaceBase *spc);
+  massRigidCylinder(rigidContactSpaceBase *spc,double leng, double radius,
                     double thick,double rho) : _spc(spc), _length(leng), _radius(radius),_thick(thick), _rho(rho){}
   // Arguments are useless but it's allow to use classical assemble function
   void get(MElement *ele, int npts, IntPt *GP, fullMatrix<double> &m) const;
@@ -40,7 +40,7 @@ class massRigidCylinder : public BilinearTermBase{
 
 class forceRigidCylinderContact : public rigidContactLinearTermBase<double>{
  protected:
-  const rigidCylinderContact *_cdom;
+  const rigidCylinderContactDomain *_cdom;
   double _fdgfactor, _facGC;
   const MVertex *_vergc;
   const SVector3* _axisDirector;
@@ -54,8 +54,8 @@ class forceRigidCylinderContact : public rigidContactLinearTermBase<double>{
   mutable double d,penetration,penpen;
   mutable SVector3 dirAC;
  public:
-  forceRigidCylinderContact(const rigidCylinderContact *cdom,
-                            rigidContactSpace *sp,const double thick,
+  forceRigidCylinderContact(const rigidCylinderContactDomain *cdom,
+                            rigidContactSpaceBase *sp,const double thick,
                             const unknownField *ufield);
   ~forceRigidCylinderContact(){}
 //  void get(MElement *ele,int npts, IntPt *GP,fullVector<double> &m);
@@ -79,7 +79,7 @@ class stiffnessRigidCylinderContact : public contactBilinearTermBase<double>{
   void get(MElement *ele, const int verIndex, fullMatrix<double> &m) const;
  protected:
   const unknownField *_ufield;
-  rigidContactSpace *_spc;
+  rigidContactSpaceBase *_spc;
   const double _perturbation; // numerical perturbation
   const double _doublepertexpm1;
   // Data for get function (Allocated Once)
@@ -91,7 +91,7 @@ class stiffnessRigidCylinderContact : public contactBilinearTermBase<double>{
   mutable std::vector<double> disp;
   mutable double pertdisp[6];
  public:
-  stiffnessRigidCylinderContact(rigidContactSpace *space, contactLinearTermBase<double> *lterm, const unknownField *ufield) : contactBilinearTermBase(lterm),
+  stiffnessRigidCylinderContact(rigidContactSpaceBase *space, contactLinearTermBase<double> *lterm, const unknownField *ufield) : contactBilinearTermBase(lterm),
                                                                                                 _ufield(ufield),
                                                                                                 _perturbation(1.e-10),
                                                                                                 _doublepertexpm1(1./(2.e-10)),

NonLinearSolver/internalPoints/ipField.h

@@ -12,7 +12,6 @@
 # ifndef _IPFIELD_H_
 # define _IPFIELD_H_
 #include<vector>
-#include"Dof3IntType.h"
 #include"quadratureRules.h"
 #include "unknownField.h"
 #include "elementField.h"
@@ -121,7 +120,7 @@ class IPField : public elementField {
   template<class T1,class T2> void getIPv(const MInterfaceElement *iele, const T1** vipv_m, const T2** vipv_p,
                                           const IPStateBase::whichState ws=IPStateBase::current) const
   {
-    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNum());
+    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
     // SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
     int npts = vips->size()/2;
     for(int i=0;i<npts;i++){
@@ -134,7 +133,7 @@ class IPField : public elementField {
   template<class T1,class T2> void getIPv(const MInterfaceElement *iele, const T1** vipv_m, const T1** vipvprev_m,
                                              const T2** vipv_p, const T2** vipvprev_p) const
   {
-    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNum());
+    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
     // SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
     int npts = vips->size()/2;
     for(int i=0;i<npts;i++){
@@ -149,7 +148,7 @@ class IPField : public elementField {
   template<class T1,class T2> void getIPv(const MInterfaceElement *iele, T1** vipv_m, T2** vipv_p,
                                              const IPStateBase::whichState ws=IPStateBase::current) const
   {
-    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNum());
+    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
     // SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
     int npts = vips->size()/2;
     for(int i=0;i<npts;i++){
@@ -162,7 +161,7 @@ class IPField : public elementField {
   template<class T1,class T2> void getIPv(const MInterfaceElement *iele, T1** vipv_m, T2** vipvprev_m,
                                              T2** vipv_p, T2** vipvprev_p) const
   {
-    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNum());
+    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
     // SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
     int npts = vips->size()/2;
     for(int i=0;i<npts;i++){
@@ -179,7 +178,7 @@ class IPField : public elementField {
   template<class T1> void getIPv(const MInterfaceElement *iele, const T1** vipv_m,
                                    const IPStateBase::whichState ws=IPStateBase::current) const
   {
-    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNum());
+    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
     // SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
     int npts = vips->size();
     for(int i=0;i<npts;i++){
@@ -189,7 +188,7 @@ class IPField : public elementField {
   }
   template<class T1> void getIPv(const MInterfaceElement *iele, const T1** vipv_m, const T1** vipvprev_m) const
   {
-    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNum());
+    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
     // SAME NUMBER OF GAUSS POINTS ON BOTH SIDES
     int npts = vips->size();
     for(int i=0;i<npts;i++){
@@ -202,7 +201,7 @@ class IPField : public elementField {
   template<class T1> void getIPv(const MInterfaceElement *iele, T1** vipv_m,
                                    const IPStateBase::whichState ws=IPStateBase::current) const
   {
-    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNum());
+    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
     int npts = vips->size();
     for(int i=0;i<npts;i++){
       IPStateBase *ipsm = (*vips)[i];
@@ -211,7 +210,7 @@ class IPField : public elementField {
   }
   template<class T1> void getIPv(const MInterfaceElement *iele, T1** vipv_m, T1** vipvprev_m) const
   {
-    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNum());
+    AllIPState::ipstateElementContainer *vips = _AIPS->getIPstate(iele->getNumber());
     int npts = vips->size();
     for(int i=0;i<npts;i++){
       IPStateBase *ipsm = (*vips)[i];

NonLinearSolver/internalPoints/ipstate.cpp

@@ -77,7 +77,7 @@ AllIPState::AllIPState(GModel *pModel, std::vector<partDomain*> &vdom)
         for(int i=npts_inter/2;i<npts_inter;i++){
           mlawPlus->createIPState(tp[i],&state,ele,iele->getElem(1)->getNumVertices());
         }
-        _mapall.insert(ipstatePairType(iele->getNum(),tp));
+        _mapall.insert(ipstatePairType(iele->getNumber(),tp));
       }
       // Virtual interface element (no duplication)
       materialLaw *mlaw = dgdom->getMaterialLaw();
@@ -89,7 +89,7 @@ AllIPState::AllIPState(GModel *pModel, std::vector<partDomain*> &vdom)
         for(int i=0;i<npts_inter;i++){
           mlaw->createIPState(tp[i],&state,ele,iele->getElem(0)->getNumVertices());
         }
-        _mapall.insert(ipstatePairType(iele->getNum(),tp));
+        _mapall.insert(ipstatePairType(iele->getNumber(),tp));
       }
     }
     // bulk element

NonLinearSolver/nlsolver/CMakeLists.txt

@@ -9,7 +9,6 @@ set(SRC
   unknownField.cpp
   unknownStaticField.cpp
   # inline
-  Dof3IntType.h
   explicitDofManager.h
   explicitHulbertChung.h
   nlsolAlgorithms.h

NonLinearSolver/nlsolver/Dof3IntType.h

-//
-// Description: Derivate class of Dof used for FullDg formulation where the type includes 3 informations (field,local vertex number, dof)
-//
-//
-// Author:  <Gauthier BECKER>, (C) 2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-#ifndef DGC0DOF_H_
-#define DGC0DOF_H_
-#include "dofManager.h"
-// Two functions are added to dof to define type with three int
-// Use by high level so impossible to let this in the projects ?
-class Dof3IntType : public Dof{
-  public :
-    Dof3IntType(long int ent, int type) : Dof(ent,type){}
-    ~Dof3IntType(){}
-    inline static int createTypeWithThreeInts(int comp,int field, int num=0){return 100000*field+100*num+comp;}
-    inline static void getThreeIntsFromType(int t, int &comp, int &field, int &num){
-      field = t/100000;
-      int i1=t%100000;
-      num = i1/100;
-      comp= i1%100;
-    }
-    // This function is used for Cg/Dg case (Allow to call the same function for the Two formulations)
-    inline static void getThreeIntsFromType(int t, int &comp, int &field){
-      field = t/100000;
-      comp =t%100000;
-    }
-};
-#endif // Dof3IntType

NonLinearSolver/nlsolver/nlsolAlgorithms.h

@@ -82,7 +82,7 @@ template<class Iterator,class Assembler> void AssembleMass(BilinearTermBase *ter
 
 // Assemble mass for a rigid contact space. Only ddl of GC
 template<class Assembler>
-void AssembleMass(BilinearTermBase *term, rigidContactSpace *space,Assembler *assembler){
+void AssembleMass(BilinearTermBase *term, rigidContactSpaceBase *space,Assembler *assembler){
   fullMatrix<typename Assembler::dataMat> localMatrix;
   explicitHCDofManager<typename Assembler::dataVec>* expAss = dynamic_cast<explicitHCDofManager<typename Assembler::dataVec>*>(assembler);
   std::vector<Dof> R;
@@ -146,7 +146,7 @@ template<class Iterator,class Assembler> void FixNodalDofs(FunctionSpaceBase *sp
 }
 
 template<class Assembler>
-void FixNodalDofs(rigidContactSpace *space,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter, Assembler &assembler){
+void FixNodalDofs(rigidContactSpaceBase *space,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter, Assembler &assembler){
   std::vector<Dof> R;
   space->getKeysOfGravityCenter(R);
   for(int i=0;i<R.size();i++){
@@ -214,7 +214,7 @@ template<class Iterator,class Assembler> void SetInitialDofs(FunctionSpaceBase *
 
 // Function Numbering Dof for rigid contact (Create Three Dofs for GC of rigid bodies)
 template<class Assembler>
-void NumberDofs(rigidContactSpace &space, Assembler &assembler){
+void NumberDofs(rigidContactSpaceBase &space, Assembler &assembler){
   // get Dofs of GC
   std::vector<Dof> R;
   space.getKeysOfGravityCenter(R);

NonLinearSolver/nlsolver/nonLinearMechSolver.cpp

@@ -91,7 +91,7 @@ void nonLinearMechSolver::initContactInteraction(){
         MElement *ele = *ite;
         if(ele == elementcont){
           flagfind =true;
-          cdom->setDomain(dom);
+          cdom->setDomainAndFunctionSpace(dom);
 //          spaceManager->addSpace(cdom);
           break;
         }
@@ -107,7 +107,7 @@ void nonLinearMechSolver::initContactInteraction(){
       for(contactContainer::iterator itcdom = _allContact.begin(); itcdom!=_allContact.end(); ++itcdom){
         contactDomain *cdom = *itcdom;
         if(cdom->getPhys() == diri._tag){
-          rigidContactSpace *sp = dynamic_cast<rigidContactSpace*>(cdom->getSpace());
+          rigidContactSpaceBase *sp = dynamic_cast<rigidContactSpaceBase*>(cdom->getSpace());
           diri.setSpace(sp);
           break;
         }
@@ -342,7 +342,7 @@ void nonLinearMechSolver::numberDofs(linearSystem<double> *lsys){
   for(contactContainer::iterator it = _allContact.begin(); it!=_allContact.end(); ++it){
     contactDomain* cdom = *it;
     if(cdom->isRigid()){
-      rigidContactSpace *rcspace = dynamic_cast<rigidContactSpace*>(cdom->getSpace());
+      rigidContactSpaceBase *rcspace = dynamic_cast<rigidContactSpaceBase*>(cdom->getSpace());
       NumberDofs(*rcspace,*pAssembler);
     }
   }
@@ -566,13 +566,8 @@ void nonLinearMechSolver::initTerms(unknownField *ufield, IPField *ipf){
   // contact domain
   for(contactContainer::iterator it = _allContact.begin(); it!=_allContact.end(); ++it){
     contactDomain *cdom = *it;
-    if(cdom->getContactType() == contactDomain::rigidCylinder){
-      rigidCylinderContact *rcc = dynamic_cast<rigidCylinderContact*>(cdom);
-      rcc->initializeTerms(ufield);
+      cdom->initializeTerms(ufield);
   }
-    else Msg::Error("Terms are not defined for contact type %d",cdom->getContactType());
-  }
-
 }
 
 void nonLinearMechSolver::solveStaticLinar()
@@ -1017,9 +1012,9 @@ void nonLinearMechSolver::archivingNode(const int numphys, const int comp, initi
   std::vector<MVertex*> vv;
   pModel->getMeshVerticesForPhysicalGroup(0,numphys,vv);
 
-  if(vv.size() !=0){
-    std::pair< Dof,initialCondition::whichCondition > pai(Dof(vv[0]->getNum(),Dof3IntType::createTypeWithThreeInts(comp,_tag)),wc);
-    anoded.push_back(pai);
+  if(vv.size() == 1){
+//    std::pair< Dof,initialCondition::whichCondition > pai(Dof(vv[0]->getNum(),Dof3IntType::createTypeWithThreeInts(comp,_tag)),wc);
+    anoded.push_back(archiveDispNode(vv[0]->getNum(),comp,_tag,wc));
     // remove old file (linux only ??)
     std::ostringstream oss;
     oss << vv[0]->getNum(); // Change this
@@ -1042,7 +1037,10 @@ void nonLinearMechSolver::archivingNode(const int numphys, const int comp, initi
     system(rfname.c_str());
   }
   else{
+    if(vv.size()==0)
       Msg::Warning("No physical group %d So it is impossible to archive nodal displacement",numphys);
+    else
+      Msg::Warning("More than one node in physical group impssible to archive for now (no loop)",numphys);
   }
 }
 
@@ -1504,7 +1502,7 @@ void nonLinearMechSolver::solveExplicit(){
   for(nonLinearMechSolver::contactContainer::iterator itc = _allContact.begin(); itc!= _allContact.end(); ++itc){
     contactDomain *cdom = *itc;
     if(cdom->isRigid()){
-      rigidContactSpace *rcsp = dynamic_cast<rigidContactSpace*>(cdom->getSpace());
+      rigidContactSpaceBase *rcsp = dynamic_cast<rigidContactSpaceBase*>(cdom->getSpace());
       AssembleMass(cdom->getMassTerm(),rcsp,pAssembler);
     }
   }

NonLinearSolver/nlsolver/nonLinearMechSolver.h

@@ -42,6 +42,15 @@ struct archiveForce{
   ~archiveForce(){}
 };
 
+struct archiveDispNode{
+  int _vernum;
+  int _comp;
+  initialCondition::whichCondition _wc;
+  int _tag;
+  archiveDispNode(const int vernum, const int comp, const int tag,initialCondition::whichCondition wc) : _vernum(vernum), _comp(comp), _wc(wc), _tag(tag){}
+  ~archiveDispNode(){}
+};
+
 struct archiveRigidContactDisp{
   int _physmaster;
   int _comp;
@@ -124,7 +133,7 @@ class nonLinearMechSolver
   std::map<int,double> aefvalue;
   std::vector<fextPrescribedDisp> _allaef;
   // std vector to archive a node displacement
-  std::vector<std::pair<Dof,initialCondition::whichCondition> > anoded;
+  std::vector<archiveDispNode> anoded;
   // std::vector to archive a force (info musy be store before creation because functionSpace are not initialize)
   std::vector<archiveForce> vaf;
   // std vector to archive ipvariable

NonLinearSolver/nlsolver/unknownDynamicField.cpp

@@ -13,8 +13,9 @@
 #include "nlsolAlgorithms.h"
 unknownDynamicField::unknownDynamicField(dofManager<double> *pas, std::vector<partDomain*> &vdom,
                                                    nonLinearMechSolver::contactContainer *acontact,
-                                     const int nc, std::vector<std::pair<Dof,initialCondition::whichCondition> > &archiving,
-                                     std::vector<archiveRigidContactDisp> &contactarch, const bool view_, const std::string filen) : unknownField(pas,vdom,acontact,nc,archiving,contactarch,view_,filen)
+                                     const int nc, std::vector<archiveDispNode> &archiving,
+                                     std::vector<archiveRigidContactDisp> &contactarch, const bool view_, const std::string filen) : unknownField(pas,vdom,acontact,nc,archiving,
+                                                                                                                                                  contactarch,view_,filen)
 {
   // to avoid constant dynamic cast
   dynassembler = dynamic_cast<explicitHCDofManager<double>*>(pAssembler);
@@ -149,8 +150,7 @@ void unknownDynamicField::archiving(const double time){
       oss << it->nodenum;
       std::string s = oss.str();
       // component of displacement
-      int field,comp,num;
-      Dof3IntType::getThreeIntsFromType(it->D.getType(),comp,field,num);
+      int comp = it->_comp;
       oss.str("");
       oss << comp;
       std::string s2 = oss.str();

NonLinearSolver/nlsolver/unknownDynamicField.h

@@ -21,7 +21,7 @@ class unknownDynamicField : public unknownField{
  public:
   unknownDynamicField(dofManager<double> *pas, std::vector<partDomain*> &vdom,
                            nonLinearMechSolver::contactContainer *acontact,
-                                     const int nc, std::vector<std::pair<Dof,initialCondition::whichCondition> > &archiving,
+                                     const int nc, std::vector<archiveDispNode> &archiving,
                                      std::vector<archiveRigidContactDisp> &contactarch, const bool view_=true, const std::string filen="disp.msh");
 
   virtual void update();