Merge pull request ElmerCSC#608 from ElmerCSC/MixedContact

Mixed contact

fem/src/modules/ElasticSolve.F90

@@ -110,7 +110,6 @@ SUBROUTINE ElasticSolver_Init( Model,Solver,dt,Transient )
       DOFs = dim
       CALL ListAddString( SolverParams, 'Variable', 'Displacement' )
     END IF
-
     CALL ListAddInteger( SolverParams, 'Variable DOFs', DOFs )
   END IF
 
@@ -2591,8 +2590,8 @@ SUBROUTINE NeoHookeanLocalMatrix( MassMatrix,DampMatrix,StiffMatrix,ForceVector,
     !------------------------------------------------------------------------------
     IF( MixedFormulation ) THEN
 
-      IF (PlaneStress) CALL Warn( Caller,  &
-          'Mixed formulation does not support plane stress: plane strain assumed instead' )
+      IF (PlaneStress) CALL Fatal( Caller,  &
+          'Mixed formulation does not support plane stress')
 
       DOFs = cdim + 1
       ! To reuse the code, set the lambda parameter to zero and instead

fem/tests/ContactPatch2DIncompressible/CMakeLists.txt

@@ -0,0 +1,8 @@
+INCLUDE(test_macros)
+INCLUDE_DIRECTORIES(${CMAKE_BINARY_DIR}/fem/src)
+
+CONFIGURE_FILE(case.sif case.sif COPYONLY)
+
+file(COPY ELMERSOLVER_STARTINFO squares.grd DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/")
+
+ADD_ELMER_TEST(ContactPatch2DIncompressible LABELS quick contact mortar elasticsolve)

fem/tests/ContactPatch2DIncompressible/ELMERSOLVER_STARTINFO

@@ -0,0 +1 @@
+case.sif

fem/tests/ContactPatch2DIncompressible/case.sif

@@ -0,0 +1,157 @@
+!------------------------------------------------------------------
+! 2d patch test.
+! Variation for incompressible case with dofs=dim+1.
+!
+! Peter Råback / 7.11.2024
+!------------------------------------------------------------------
+
+Header
+  CHECK KEYWORDS Warn
+  Mesh DB "." "squares"
+  Include Path ""
+  Results Directory ""
+End
+
+$fileid="b"
+
+Simulation
+  Max Output Level = 7
+  Coordinate System = Cartesian
+  Coordinate Mapping(3) = 1 2 3
+  Simulation Type = Steady State
+
+  Steady State Min Iterations = 1
+  Steady State Max Iterations = 1
+
+  ascii output = true
+  Post File = case_$fileid$.vtu
+  Save Geometry Ids = Logical True
+
+! The ElasticSolver does not really like the Dirichlet conditions at the start 
+! of the nonlinear iteration. 
+  Initialize Dirichlet Conditions = False
+End
+
+Equation 1
+  Name = "Deformation"
+  Active Solvers(1) = 1
+End
+
+Body 1
+  Name = "Lower block"
+  Target Bodies(1) = 1
+  Equation = 1
+  Material = 1
+End
+
+Body 2
+  Name = "Upper block"
+  Target Bodies(1) = 2
+  Equation = 1
+  Material = 1
+End
+
+Material 1
+  Name = "Ideal"
+  Youngs modulus = 90.0
+  Density = 10.0
+  Poisson ratio = 0.25
+End
+
+Solver 1
+  Equation = "Nonlinear elasticity"
+  Procedure = "ElasticSolve" "ElasticSolver"
+!  Variable = -dofs 2 Displacement
+
+! Have pressure as variable also!
+  Mixed Formulation = Logical True
+  Neo-Hookean Material = Logical True
+  Bubbles in Global System = False
+
+  Element = p:1 b:3
+
+  Nonlinear System Convergence Tolerance = 1.0e-5
+  Nonlinear System Max Iterations = 10
+  Nonlinear System Relaxation Factor = 1.0
+
+  Linear System Solver = "Iterative"
+  Linear System Iterative Method = "BiCGStab"
+  Linear System Abort Not Converged = True
+  Linear System Preconditioning = "ILU2"
+  Linear System Residual Output = 100
+  Linear System Max Iterations = 5000
+  BiCGStabl Polynomial Degree = 4
+
+  Linear System Convergence Tolerance = 1.0e-10
+
+  Apply Contact BCs = Logical True
+! Save Contact = Logical True
+
+! Restore the linear solution
+! Elasticity Solver Linear = Logical True
+
+  Calculate Stresses = Logical True
+! Optimize Bandwidth = False
+
+  Displace Mesh = Logical True
+
+! Do not include constraints when analyzing the convergence and norm of a solution
+  Nonlinear System Convergence Without Constraints = Logical True
+End
+
+Solver 2
+  Exec Solver = never
+  Equation = "SaveLine"
+  Procedure = "SaveData" "SaveLine"
+  Filename = f_$fileid$.dat
+End
+
+Boundary Condition 1
+  Name = "Support"
+  Target Boundaries(1) = 1
+  Displacement 2 = Real 0.0
+  Disp 2 = Real 0.0
+End
+
+Boundary Condition 2
+  Name = "Lower surface of upper block"
+  Target Boundaries(1) = 5
+
+  Mortar BC = Integer 3
+  Mortar BC Nonlinear = Logical True
+  Contact Depth Offset Initial = Real 1.0e-3
+  !Contact Active Set Minimum = Integer 1
+  !Contact No-Slip = Logical True
+
+! Create a strong projector for the line setting y-coordinate to zero
+  Flat Projector = Logical True
+
+! a) Use weak projector or not
+  Galerkin Projector = Logical True
+
+! b) Use more tailored projector able to do accurate integration
+  Level Projector = Logical True
+  Level Projector Generic = True
+End
+
+Boundary Condition 3
+  Name = "Upper surface of lower block"
+  Target Boundaries(1) = 3
+End
+
+Boundary Condition 4
+  Name = "Pressure load the upper surface of upper block"
+  Target Boundaries(1) = 7
+  Normal Surface Traction = -1.0
+End
+
+Boundary Condition 5
+  Name = "Symmetry"
+  Target Boundaries(2) = 4 8
+  Displacement 1 = 0.0
+  Disp 1 = Real 0.0
+End
+
+Solver 1 :: Reference Norm = 1.44268191E-01
+Solver 1 :: Reference Norm Tolerance = 1.0e-5
+

fem/tests/ContactPatch2DIncompressible/runtest.cmake

@@ -0,0 +1,3 @@
+include(test_macros)
+execute_process(COMMAND ${ELMERGRID_BIN} 1 2 squares.grd)
+RUN_ELMER_TEST()

fem/tests/ContactPatch2DIncompressible/squares.grd

@@ -0,0 +1,43 @@
+#####  ElmerGrid input file for structured grid generation  ######
+Version = 210903
+Coordinate System = Cartesian 2D
+Subcell Divisions in 2D = 1 2 
+Subcell Sizes 1 = 1 
+Subcell Sizes 2 = 1 1
+Material Structure in 2D
+  2  
+  1  
+End
+Materials Interval = 1 1
+Boundary Definitions
+# type     out      int     
+  1       -1        1        0       
+  2       -2        1        0       
+  3       -3        1        0       
+  4       -4        1        0       
+End
+Element Degree = 1
+Element Innernodes = False
+Triangles = False
+Element Ratios 1 = 1 
+Element Ratios 2 = 1 1
+Element Divisions 1 = 6
+Element Divisions 2 = 6 0
+
+Start New Mesh
+
+Materials Interval = 2 2
+Boundary Definitions
+# type     out      int     
+  5       -1        2        0       
+  6       -2        2        0       
+  7       -3        2        0       
+  8       -4        2        0       
+End
+Element Degree = 1
+Element Innernodes = False
+Triangles = False
+Element Divisions 1 = 7
+Element Divisions 2 = 0 7
+
+Unite = True