fix bug on _lastNormalContact

NonLinearSolver/contact/contactDomain.h

@@ -226,7 +226,7 @@ class rigidPlaneContactDomain : public contactDomain{
 class rigidSphereContactDomain : public contactDomain{
 
  protected:
-  MVertex *_vergc;  // vertex of gravity center (center of sphere)
+  MVertex *_vergc;  // vertex of gravity center of sphere
   double _radius; // radius of sphere;
 //  double _thickContact; // (?) use for shell (contact with neutral axis of external fiber is !=0)
   double _density; // density of sphere (? Not a material law for now ?)
@@ -238,7 +238,7 @@ class rigidSphereContactDomain : public contactDomain{
                            const double penalty, const double rho, elementFilter *filSlave=NULL);
 
   rigidSphereContactDomain(const rigidSphereContactDomain& src): contactDomain(src),_vergc(src._vergc),
-                           _radius(src._radius),_density(src._density){}; // (? 1.a necessary constructor; 2.src ?)
+                           _radius(src._radius),_density(src._density){};
   virtual ~rigidSphereContactDomain(){}
   virtual MVertex* getGC() const{return _vergc;};
   virtual void setDomainAndFunctionSpace(partDomain *dom)=0; // (?)

NonLinearSolver/contact/rigidSphereContactTerms.cpp

@@ -61,12 +61,12 @@ void forceRigidSphereContact::get(const MVertex *ver, const double disp[6],doubl
   A = ver->point();
   Adisp.setPosition(disp[0],disp[1],disp[2]);
   A+=Adisp;
-  SVector3 BA(B,A);
-  double d = BA.norm();
+  SVector3 _lastNormalContact(B,A); // direction : from Master to Slave
+  double d = _lastNormalContact.norm();
 
   // test to know if contact
   if(d < _radius){
-    _lastNormalContact = BA.normalize(); // direction : from Master to Slave
+    _lastNormalContact.normalize(); 
     double penetration = _radius-d;
     double penpen = penetration*_penalty;
     for(int j=0;j<3;j++){

dG3D/benchmarks/panelSphereContact/simplifypanelSphereContact.geo → dG3D/benchmarks/cubeSphereContact/cubeSphereContact.geo

 // GEOMETRY AND MESH PARAMETERS
-lays = 10; // number of layers of the composite (must be an integer)
+lays = 5; // number of layers of the composite (must be an integer, in z direction)
 mm=1e-3; // unit for geometry
-Lx = 150*mm; // panel length in x direction
-Ly = 150*mm; // panel length in y direction
-Lz = 100*mm; // length in z direction
-R = 20*mm; // sphere radius
+Lx = 100*mm; // dimension in x direction
+Ly = 100*mm; // dimension in y direction
+Lz = 50*mm; // dimension in z direction
+R = 15*mm; // sphere radius
 dis = 0.1*mm; // initial contact distance 
-slp = 10*mm; // panel characteristic length
+slc = 10*mm; // cube characteristic length
 sls = 10*mm; // sphere characteristic length
-pT = 15; // panel transfinite line value (control nodes number by line)
+cT = 10; // panel transfinite line value (control nodes number by line)
 sT = 3; // sphere transfinite line value (control nodes number by line)
 //
 //////////////////////////
-//    CompositePanel	//
+//    CompositeCube	//
 //////////////////////////
 // GEOMETRY
-// panel
-Point(1)={-Lx/2,-Ly/2,0,slp};
-Point(2)={-Lx/2,Ly/2,0,slp};
-Point(3)={Lx/2,Ly/2,0,slp};
-Point(4)={Lx/2,-Ly/2,0,slp};
+Point(1)={-Lx/2,-Ly/2,0,slc};
+Point(2)={-Lx/2,Ly/2,0,slc};
+Point(3)={Lx/2,Ly/2,0,slc};
+Point(4)={Lx/2,-Ly/2,0,slc};
 
 Line(1)={1,2};
 Line(2)={2,3};
@@ -28,11 +27,10 @@ Line(4)={4,1};
 
 // surfaces
 Line Loop(1)={1,2,3,4};
-
 Plane Surface(1)={1};
 
 // replace unstructured grid by structured grid (with controlled nodes number by line)
-Transfinite Line {1,2,3,4} = pT Using Progression 1;
+Transfinite Line {1,2,3,4} = cT Using Progression 1;
 // transfinite surface guarantee a structured grid on a surface
 Transfinite Surface {1};
 
@@ -52,14 +50,14 @@ my_plate[]=Extrude {0, 0, -Lz} {
 
 // PHYSICAL GROUPS
 Physical Surface(81) = {13};
-Physical Surface(82) = {21};
-Physical Surface(83) = {17};
+Physical Surface(82) = {17};
+Physical Surface(83) = {21};
 Physical Surface(84) = {25};
-Physical Volume(85) = { my_plate[] }; //
+Physical Volume(85) = { my_plate[] };
 
 //
 //////////////////////////
-//        Sphere    	//
+//      RigidSphere  	//
 //////////////////////////
 // GEOMETRY
 Point(41) = {0,0,R+dis,sls};
@@ -134,6 +132,5 @@ Transfinite Surface {73,74,75,76,77,78,79,80};
 // Recombine Volume {3};
 
 // PHYSICAL GROUPS
-Physical Point(87) ={41}; // sphere center
-Physical Surface(80) ={78, 77, 80, 79, 74, 73, 76, 75};
-
+Physical Point(86) ={41}; // sphere center
+Physical Surface(87) ={78, 77, 80, 79, 74, 73, 76, 75}; // sphere surface

dG3D/benchmarks/panelSphereContact/simplifypanelSphereContact.py → dG3D/benchmarks/cubeSphereContact/cubeSphereContact.py

@@ -4,30 +4,31 @@ from gmshpy import *
 from dG3DpyDebug import*
 
 #script to launch PBC problem with a python script 
+#(! radius need to be the same value as in .geo !)
 
 # parameters for material law
-lawnum = 11 # unique number of law (?)
+lawnum = 11 # unique number of law
+rho = 100 # density (? unit ?)
 E = 3.e9 # Young's modulus Pa
 nu= 0.4 # Poisson ratio
 sy0 = 120e6 # initial yield stress Pa
 h = E/10. # (?)
-rho = 100 # density (? unit ?)
 
 # creation of material law
-law1 = J2LinearDG3DMaterialLaw(lawnum,rho,E,nu,sy0,h) # (?)
+law1 = J2LinearDG3DMaterialLaw(lawnum,rho,E,nu,sy0,h) # (dG3D/src/dG3DMaterialLaw.h/.cpp)
 
 # geometry & mesh
-geofile="simplifypanelSphereContact.geo"  # name of the geometry with physical groups
-meshfile="simplifypanelSphereContact.msh" # name of mesh file (to save directly in gmsh shift+ctrl+s)
+geofile="cubeSphereContact.geo"  # name of the geometry with physical groups
+meshfile="cubeSphereContact.msh" # name of mesh file (to save directly in gmsh shift+ctrl+s)
 
 # parameters for part domain
-nfield = 85 # number of the field (physical number of entity)
+nfield = 85 # number of the field (physical number of entity defined in gmsh)
 #myfield1 = FSDomain(1000,nfield,lawnum,3) # (?)
-fullDg = 0 # not a Discontinuous Galerkin part Domain (true=1/false=0 ? if plies =1 ?)
+fullDg = 0 # not a Discontinuous Galerkin part Domain
 beta1 = 1e2 # stability parameters (?)
 
 # creation of part domain
-myfield1 = dG3DDomain(1000,nfield,0,lawnum,fullDg,3) # (tag=1000 : tag number for slave domains with no linked ddls ; ws=0 : which state == initial)
+myfield1 = dG3DDomain(1000,nfield,0,lawnum,fullDg,3) # (dG3D/src/dG3DDomain.h/.cpp : tag=1000 : tag number for slave domains with no linked ddls ; ws=0 : functionSpaceType::Lagrange ; dim=3)
 myfield1.stabilityParameters(beta1)
 
 # parameters for solver
@@ -65,26 +66,26 @@ mysolver.displacementBC("Face",84,1,0.)
 mysolver.displacementBC("Face",84,2,0.)
 
 # parameters for contact definition
-radius=8.e-3
-penalty=1.e4 # (? unit ?)
+radius=15.e-3 # ! same value as in .geo !
+penalty=1.e8 # (? unit ?)
 rhocontact=7800.
 sphereDispl = -20.e-3
 contactDom1 = 1000 # (tag number for master domain)
 
 # creation of contact definition
 flaw = CoulombFrictionLaw(2,0.,0.,penalty,penalty/1000.) # normal to the contact, static friction coeff, dynamic friction coeff
-contact1 = dG3DRigidSphereContactDomain(contactDom1, 2, 80, 3, nfield, 87,radius,penalty,rhocontact) # (? 1.dimMaster=0; 2.physpt=center of sphere ?)
+contact1 = dG3DRigidSphereContactDomain(contactDom1, 2, 87, 3, nfield, 86,radius,penalty,rhocontact)
 contact1.setFriction(bool(0))
 contact1.addFrictionLaw(flaw)
 mysolver.contactInteraction(contact1)
 
-# BC sphere (? Physical Surface(86) or Physical Volume)
-mysolver.displacementRigidContactBC(80,0,0.)
-mysolver.displacementRigidContactBC(80,1,0.)
-mysolver.displacementRigidContactBC(80,2,sphereDispl)
-#mysolver.displacementRigidContactBC(87,0,0.)
-#mysolver.displacementRigidContactBC(87,1,0.)
-#mysolver.displacementRigidContactBC(87,2,sphereDispl)
+# BC sphere
+mysolver.displacementRigidContactBC(87,0,0.)
+mysolver.displacementRigidContactBC(87,1,0.)
+mysolver.displacementRigidContactBC(87,2,sphereDispl)
+#mysolver.displacementRigidContactBC(86,0,0.)
+#mysolver.displacementRigidContactBC(86,1,0.)
+#mysolver.displacementRigidContactBC(86,2,sphereDispl)
 
 #
 ## Outputs
@@ -108,7 +109,7 @@ mysolver.internalPointBuildView("sig_VM",IPField.SVM, 1, 1);
 mysolver.internalPointBuildView("Green-Lagrange equivalent strain",IPField.GL_EQUIVALENT_STRAIN, 1, 1);
 mysolver.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN, 1, 1);
 
-mysolver.archivingRigidContact(87,2,0, nstepArch);
+mysolver.archivingRigidContact(86,2,0, nstepArch);
 
 # mysolver.archivingIPOnPhysicalGroup
 # mysolver.archivingNodeDisplacement

dG3D/benchmarks/panelSphereContact/panelSphereContact.geo

-// GEOMETRY AND MESH PARAMETERS
-lays = 8; // number of layers of the composite (must be an integer)
-mm=1e-3; // unit for geometry
-Lx = 150*mm; // panel length in x direction (coarse mesh)
-Ly = 100*mm; // panel length in y direction (coarse mesh)
-Lz = (3.08/lays)*mm; // length in z direction
-lx = 45*mm; // sub-panel length in x direction (fine mesh)
-ly = 41*mm; // sub-panel length in y direction (fine mesh)
-R = 8*mm; // sphere radius
-dis = 1*mm; // initial contact distance 
-sl1 = 20*mm; // panel characteristic length
-sl2 = 5*mm;  // sub-panel characteristic length
-sl3 = 3*mm; // sphere characteristic length
-// sl3 = R;
-pT = 20; // panel transfinite line value (control nodes number by line)
-sT = 3; // sphere transfinite line value (control nodes number by line)
-//
-//////////////////////////
-//    CompositePanel	//
-//////////////////////////
-// GEOMETRY
-// panel
-Point(1)={-Lx/2,-Ly/2,0,sl1};
-Point(2)={-Lx/2,Ly/2,0,sl1};
-Point(3)={Lx/2,Ly/2,0,sl1};
-Point(4)={Lx/2,-Ly/2,0,sl1};
-
-Line(1)={1,2};
-Line(2)={2,3};
-Line(3)={3,4};
-Line(4)={4,1};
-
-// sub-panel
-Point(5)={-lx/2,-ly/2,0,sl2};
-Point(6)={-lx/2,ly/2,0,sl2};
-Point(7)={lx/2,ly/2,0,sl2};
-Point(8)={lx/2,-ly/2,0,sl2};
-
-Line(5)={5,6};
-Line(6)={6,7};
-Line(7)={7,8};
-Line(8)={8,5};
-
-// surfaces
-Line Loop(1)={1,2,3,4};
-Line Loop(2)={5,6,7,8};
-
-Plane Surface(1)={1,2};
-Plane Surface(2)={2};
-
-// replace unstructured grid by structured grid (with controlled nodes number by line)
-Transfinite Line {5,6,7,8} = pT Using Progression 1;
-// ? transfinite surface is used to define the ordering of the nodes/elements in the mesh ?
-Transfinite Surface {2};
-
-// replace tri/tetra by quad/hexa elements
-Recombine Surface{1,2};
-
-// extrusion of surfaces along z direction
-my_plate[]=Extrude {0, 0, -Lz} {
-  Surface{1}; Surface{2};
-  Layers{lays};
-  Recombine; 
-};
-
-// ? transfinite volume is used to define the ordering of the nodes/elements in the mesh ?
-// ? transfinite volume not needed when using extrusion ? (no effet)
-// Transfinite Volume {1,2};
-
-// PHYSICAL GROUPS
-Physical Surface(81) = {21};
-Physical Surface(82) = {33};
-Physical Surface(83) = {29};
-Physical Surface(84) = {25};
-Physical Volume(85) = { my_plate[] }; //
-
-//
-//////////////////////////
-//        Sphere    	//
-//////////////////////////
-// GEOMETRY
-Point(41) = {0,0,R+dis,sl3};
-Point(42) = {R,0,R+dis,sl3};
-Point(43) = {-R,0,R+dis,sl3};
-Point(44) = {0,R,R+dis,sl3};
-Point(45) = {0,-R,R+dis,sl3};
-Point(46) = {0,0,2*R+dis,sl3};
-Point(47) = {0,0,dis,sl3};
-
-Circle(51) = {42,41,44};
-Circle(52) = {44,41,43};
-Circle(53) = {43,41,45};
-Circle(54) = {45,41,42};
-Circle(55) = {42,41,46};
-Circle(56) = {46,41,43};
-Circle(57) = {43,41,47};
-Circle(58) = {47,41,42};
-Circle(59) = {44,41,46};
-Circle(60) = {46,41,45};
-Circle(61) = {45,41,47};
-Circle(62) ={47,41,44};
-
-//+ surfaces
-Curve Loop(3) = {51, -62, 58};
-//+
-Surface(73) = {3};
-//+
-Curve Loop(4) = {51, 59, -55};
-//+
-Surface(74) = {4};
-//+
-Curve Loop(5) = {52, -56, -59};
-//+
-Surface(75) = {5};
-//+
-Curve Loop(6) = {52, 57, 62};
-//+
-Surface(76) = {6};
-//+
-Curve Loop(7) = {53, 61, -57};
-//+
-Surface(77) = {7};
-//+
-Curve Loop(8) = {53, -60, 56};
-//+
-Surface(78) = {8};
-//+
-Curve Loop(9) = {54, 55, 60};
-//+
-Surface(79) = {9};
-//+
-Curve Loop(10) = {54, -58, -61};
-//+
-Surface(80) = {10};
-//+
-Surface Loop(1) = {78, 77, 80, 79, 74, 73, 76, 75};
-//+ volume
-Volume(3) = {1};
-
-//replace unstructured grid by structured grid (with controlled nodes number by line)
-Transfinite Line {51,52,53,54,55,56,57,58,59,60,61,62} = sT Using Progression 1;
-
-// ? transfinite surface and volume is used to define the ordering of the nodes/elements in the mesh ?
-Transfinite Surface {73,74,75,76,77,78,79,80};
-// transfinite algorithm only available for 5- and 6-face volumes
-// Transfinite Volume {3};
-
-// replace tri/tetra by quad/hexa elements
-Recombine Surface {73,74,75,76,77,78,79,80};
-Recombine Volume {3};
-
-// PHYSICAL GROUPS
-Physical Surface(86) = {73,74,75,76,77,78,79,80}; // sphere
-Physical Point(87) ={41}; // sphere center
-
-

dG3D/benchmarks/panelSphereContact/panelSphereContact.py

-#coding-Utf-8-*-
-from gmshpy import *
-#from dG3Dpy import*
-from dG3DpyDebug import*
-
-#script to launch PBC problem with a python script
-
-# material law
-lawnum = 11 # unique number of law (?)
-E = 3.e9 # Young's modulus Pa
-nu= 0.4 # Poisson ratio
-sy0 = 120e6 # initial yield stress Pa (?)
-h = E/10. # (?)
-rho = 100 # (?)
-
-# creation of material law
-law1 = J2LinearDG3DMaterialLaw(lawnum,rho,E,nu,sy0,h) # (?)
-
-# geometry & mesh
-geofile="panelSphereContact.geo"  # name of the geometry with physical groups
-meshfile="panelSphereContact.msh" # name of mesh file (to save directly in gmsh shift+ctrl+s)
-
-# creation of part Domain
-nfield = 85 # number of the field (physical number of entity)
-#myfield1 = FSDomain(1000,nfield,lawnum,3) # (?)
-fullDg = 0 # not a Discontinuous Galerkin part Domain (?)
-beta1 = 1e2 # stability parameters
-myfield1 = dG3DDomain(1000,nfield,0,lawnum,fullDg,3) (? 1000 ?)
-myfield1.stabilityParameters(beta1)
-
-# solver parameters
-sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2 (?)
-soltype =1 # StaticLinear=0 (default) StaticNonLinear=1
-nstep = 50  # number of step (used only if soltype=1)
-ftime =1.   # Final time (used only if soltype=1)
-tol=1.e-4  # relative tolerance for NR scheme (used only if soltype=1)
-tolAbs = 1e-12
-nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
-
-# creation of Solver
-mysolver = nonLinearMechSolver(1000) # tag=1000 (?)
-mysolver.loadModel(meshfile)
-mysolver.addDomain(myfield1)
-mysolver.addMaterialLaw(law1)
-
-mysolver.Scheme(soltype)
-mysolver.Solver(sol)
-mysolver.snlData(nstep,ftime,tol,tolAbs)
-mysolver.stepBetweenArchiving(nstepArch)
-
-# BC plate
-mysolver.displacementBC("Face",81,0,0.)
-mysolver.displacementBC("Face",81,1,0.)
-mysolver.displacementBC("Face",81,2,0.) # (?) zero displacement in z direction 
-mysolver.displacementBC("Face",82,0,0.)
-mysolver.displacementBC("Face",82,1,0.)
-mysolver.displacementBC("Face",82,2,0.) # (?) zero displacement in z direction 
-mysolver.displacementBC("Face",83,0,0.)
-mysolver.displacementBC("Face",83,1,0.)
-mysolver.displacementBC("Face",83,2,0.) # (?) zero displacement in z direction 
-mysolver.displacementBC("Face",84,0,0.)
-mysolver.displacementBC("Face",84,1,0.)
-mysolver.displacementBC("Face",84,2,0.) # (?) zero displacement in z direction
-
-# Contact definition
-penalty=1.e4 # (?)
-hcontact=0.5e-3 # (?) semi height of a virtual box placed on the contact surface and where contact is taken into account
-rhocontact=7800. # (?)
-sphereDispl = -1.e-3 # (?)
-
-contactDom1= 6000 # (?)
-
-# flaw1 = CoulombFrictionLaw(1,0.,1.,1e11,1e10)
-
-contact1 = dG3DRigidSphereContact(contactDom1, 2, 86, 3, nfield, 87,penalty,hcontact,rhocontact) # (? 1.dimMaster=2/3; 2.physpt=center of sphere ?)
-# (?) contact1.setFriction(bool(0))
-# (?) contact1.addFrictionLaw(flaw)
-mysolver.contactInteraction(contact1) # (?)
-
-# BC sphere (? Physical Surface(86) or Physical Volume)
-mysolver.displacementRigidContactBC(86,0,0.)
-mysolver.displacementRigidContactBC(86,1,0.)
-mysolver.displacementRigidContactBC(86,2,sphereDispl)
-
-#
-## Outputs
-mysolver.internalPointBuildView("Green-Lagrange_xx",IPField.STRAIN_XX, 1, 1);
-mysolver.internalPointBuildView("Green-Lagrange_yy",IPField.STRAIN_YY, 1, 1);
-mysolver.internalPointBuildView("Green-Lagrange_zz",IPField.STRAIN_ZZ, 1, 1);
-mysolver.internalPointBuildView("Green-Lagrange_xy",IPField.STRAIN_XY, 1, 1);
-mysolver.internalPointBuildView("Green-Lagrange_yz",IPField.STRAIN_YZ, 1, 1);
-mysolver.internalPointBuildView("Green-Lagrange_xz",IPField.STRAIN_XZ, 1, 1);
-mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
-mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
-mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
-mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
-mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
-mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
-mysolver.internalPointBuildView("sig_VM",IPField.SVM, 1, 1);
-mysolver.internalPointBuildView("Green-Lagrange equivalent strain",IPField.GL_EQUIVALENT_STRAIN, 1, 1);
-mysolver.internalPointBuildView("Equivalent plastic strain",IPField.PLASTICSTRAIN, 1, 1);
-
-# mysolver.archivingIPOnPhysicalGroup
-# mysolver.archivingNodeDisplacement
-# mysolver.archivingDispOnPhysicalGroup
-# mysolver.archivingForceOnPhysicalGroup
-# mysolver.archivingRigidContact
-# mysolver.archivingRigidContactForce
-
-
-mysolver.solve()
-
-# check = TestCheck() # (?)
-

dG3D/src/dG3DContactDomain.cpp

@@ -43,8 +43,7 @@ void dG3DRigidPlaneContactDomain::setDomainAndFunctionSpace(partDomain *dom){
 
 dG3DRigidSphereContactDomain::dG3DRigidSphereContactDomain(const int tag, const int dimMaster, const int physMaster,
                              const int dimSlave, const int physSlave, const int physPointBase,
-                             const double radius, const double penalty, const double rho, 
-                             elementFilter *fslave):
+                             const double radius, const double penalty, const double rho, elementFilter *fslave):                            
                              rigidSphereContactDomain(tag, dimMaster, physMaster, dimSlave, physSlave, physPointBase,
                              radius, penalty, rho, fslave){};