Collapsed element export bug fix, plus draft cavity volume calculations

src/boundary_condition_routines.f90

@@ -157,7 +157,9 @@ MODULE BOUNDARY_CONDITIONS_ROUTINES
 
   PUBLIC BoundaryConditions_ConstrainNodeDofsEqual
 
-CONTAINS  
+  PUBLIC BoundaryConditions_ConstrainNodeDofLinear
+
+CONTAINS
 
   !
   !================================================================================================================================
@@ -3049,6 +3051,68 @@ SUBROUTINE BoundaryConditions_ConstrainNodeDofsEqual( &
     RETURN 1
   END SUBROUTINE BoundaryConditions_ConstrainNodeDofsEqual
 
+  !
+  !================================================================================================================================
+  !
+
+  !>Constrain a nodal equations dependent field DOF to be another single solver DOF with a mapping coeff in the solver equations
+  SUBROUTINE BoundaryConditions_ConstrainNodeDofLinear(boundaryConditions,field,fieldVariableType, &
+      & versionNumbers,derivativeNumbers,componentNumbers,nodeNumbers,coefficient,err,error,*)
+
+    !Argument variables
+    TYPE(BOUNDARY_CONDITIONS_TYPE), POINTER, INTENT(IN) :: boundaryConditions !<The solver equations boundary conditions to constrain the DOFs for.
+    TYPE(FIELD_TYPE), POINTER, INTENT(IN) :: field !<The equations dependent field containing the field DOFs to be constrained.
+    INTEGER(INTG), INTENT(IN) :: fieldVariableType !<The field variable type of the DOFs to be constrained. \see OPENCMISS_FieldVariableTypes
+    INTEGER(INTG), INTENT(IN) :: versionNumbers(:) !<The derivative version number.
+    INTEGER(INTG), INTENT(IN) :: derivativeNumbers(:) !<The derivative number.
+    INTEGER(INTG), INTENT(IN) :: componentNumbers(:) !<The field component number of the DOFs to be constrained.
+    INTEGER(INTG), INTENT(IN) :: nodeNumbers(:) !<The user numbers of the nodes to be constrained to be equal.
+    REAL(DP), INTENT(IN) :: coefficient !<The mapping coefficient
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code.
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error message.
+    !Local variables
+    TYPE(FIELD_VARIABLE_TYPE), POINTER :: fieldVariable
+    INTEGER(INTG) :: numberOfDofs,localDofNumber,dofIdx
+    INTEGER(INTG), ALLOCATABLE :: globalDofs(:)
+
+    CALL Enters("BoundaryConditions_ConstrainNodeDofLinear",err,error,*998)
+
+    NULLIFY(fieldVariable)
+
+    IF(.NOT.ASSOCIATED(boundaryConditions)) THEN
+      CALL FlagError("Boundary conditions are not associated.",err,error,*998)
+    END IF
+
+    ! check if all input arguments have the same length
+    numberOfDofs=2
+    IF(SIZE(versionNumbers,1)==numberOfDofs .AND. SIZE(derivativeNumbers,1)==numberOfDofs .AND. &
+        & SIZE(componentNumbers,1)==numberOfDofs .AND. SIZE(nodeNumbers,1)==numberOfDofs) THEN
+      ALLOCATE(globalDofs(numberOfDofs),stat=err)
+      IF(err/=0) CALL FlagError("Could not allocate equal global DOFs array.",err,error,*998)
+      !Get field DOFs for nodes
+      DO dofIdx=1,numberOfDofs
+        CALL FIELD_COMPONENT_DOF_GET_USER_NODE(field,fieldVariableType,versionNumbers(dofIdx),derivativeNumbers(dofIdx), &
+          & nodeNumbers(dofIdx),componentNumbers(dofIdx),localDofNumber,globalDofs(dofIdx),err,error,*999)
+      ENDDO !dofIdx
+      !Get the field variable and boundary conditions variable for the field
+      CALL FIELD_VARIABLE_GET(field,fieldVariableType,fieldVariable,err,error,*999)
+      !Now set DOF constraint
+      CALL BoundaryConditions_DofConstraintSet( &
+        & boundaryConditions,fieldVariable,globalDofs(2),[globalDofs(1)],[coefficient],err,error,*999)
+      DEALLOCATE(globalDofs)
+    ELSE
+      CALL FlagError("Input version, derivative node and component numbers must be arrays of length 2.",err,error,*998)
+    ENDIF
+
+    CALL Exits("BoundaryConditions_ConstrainNodeDofLinear")
+    RETURN
+999 IF(ALLOCATED(globalDofs)) DEALLOCATE(globalDofs)
+998 CALL Errors("BoundaryConditions_ConstrainNodeDofLinear",err,error)
+    CALL Exits("BoundaryConditions_ConstrainNodeDofLinear")
+    RETURN 1
+  END SUBROUTINE BoundaryConditions_ConstrainNodeDofLinear
+
+
   !
   !================================================================================================================================
   !

src/field_IO_routines.f90

@@ -3207,7 +3207,7 @@ SUBROUTINE FIELD_IO_EXPORT_ELEMENTAL_GROUP_HEADER_FORTRAN( global_number, MAX_NO
                   DO nx=1,NodesX
                     nn=nn+1
                     IF (nx==1) THEN
-                      NODE_NUMBER_COLLAPSED=NodesX*(ny-1)+ny
+                      NODE_NUMBER_COLLAPSED=NodesX*(ny-1)+nx
                       NODE_NUMBER=NODE_NUMBER_COLLAPSED
                     ELSE
                       IF (NODE_NUMBER_COUNTER<NODE_NUMBER_COLLAPSED) THEN
@@ -3241,12 +3241,13 @@ SUBROUTINE FIELD_IO_EXPORT_ELEMENTAL_GROUP_HEADER_FORTRAN( global_number, MAX_NO
                 ENDDO
               ENDDO
             ELSE IF(BASIS%COLLAPSED_XI(1)==BASIS_COLLAPSED_AT_XI1) THEN
+              NODE_NUMBER_COLLAPSED = 0
               DO nz=1,NodesZ
                 DO ny=1,NodesY
                   DO nx=1,NodesX
                     nn=nn+1
                     IF (nx==NodesX) THEN
-                      NODE_NUMBER_COLLAPSED=NodesX*(ny-1)+ny+NodesX-1
+                      NODE_NUMBER_COLLAPSED=NodesX*(ny-1)+nx
                       NODE_NUMBER=NODE_NUMBER_COLLAPSED
                     ELSE
                       IF (NODE_NUMBER_COUNTER<NODE_NUMBER_COLLAPSED) THEN

src/field_routines.f90

@@ -749,6 +749,8 @@ MODULE FIELD_ROUTINES
   PUBLIC FIELD_VARIABLE_TYPES_SET,FIELD_VARIABLE_TYPES_SET_AND_LOCK
 
   PUBLIC MESH_EMBEDDING_PUSH_DATA, MESH_EMBEDDING_PULL_GAUSS_POINT_DATA, FIELD_PARAMETER_SET_GET_GAUSS_POINT_COORD
+
+  PUBLIC Field_CalculateEnclosedVolume
 
 CONTAINS
 
@@ -11553,6 +11555,88 @@ END SUBROUTINE FIELD_GEOMETRIC_PARAMETERS_SCALE_FACTORS_UPDATE
   !================================================================================================================================
   !
 
+  !> Evaluate volume enclosed by specified external faces of a geometric field by integrating the product of surface
+  ! area and distance from a fixed point within cavity space - approximating cone shapes.
+  SUBROUTINE Field_CalculateEnclosedVolume(field, point, elems, faces, volume, err, error,*)
+
+    ! Argument variables
+    TYPE(FIELD_TYPE), POINTER :: field !<Evaluate enclosed volume for this field.
+    REAL(DP), INTENT(IN) :: point(:) !<Global coordinates of arbitrary point within the enclosed volume.
+    INTEGER(INTG), INTENT(IN) :: elems(:) !<Array of elements which face the enclosed volume to be evaluated.
+    INTEGER(INTG), INTENT(IN) :: faces(:) !<Face number of each element specified in elems array which faces the enclosed volume.
+    REAL(DP), INTENT(OUT) :: volume !<Total volume of enclosed space.
+    INTEGER(INTG), INTENT(OUT) :: err !<The error code.
+    TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string.
+
+    ! Local Variables
+    REAL(DP) :: xi_gauss(10), w_gauss(10), xi(1), w, r
+    INTEGER(INTG) :: numberOfGaussPoints=4, ng, elemIndex
+    INTEGER(INTG) :: gaussStart(4)=[0,1,3,6]
+    TYPE(COORDINATE_SYSTEM_TYPE), POINTER :: coordinateSystem
+    TYPE(FIELD_INTERPOLATED_POINT_PTR_TYPE), POINTER :: interpolatedPoint(:)
+    TYPE(FIELD_INTERPOLATED_POINT_METRICS_PTR_TYPE), POINTER :: interpolatedPointMetrics(:)
+    TYPE(FIELD_INTERPOLATION_PARAMETERS_PTR_TYPE), POINTER :: interpolationParameters(:)
+
+    xi_gauss = [ 0.500000000000000_DP, &
+      &     0.211324865405187_DP,0.788675134594813_DP, &
+      &     0.112701665379258_DP,0.500000000000000_DP,0.887298334620742_DP, &
+      &     0.06943184420297349_DP,0.330009478207572_DP,0.669990521792428_DP,0.930568155797026_DP ]
+    w_gauss = [ 1.000000000000000_DP, &
+      &     0.500000000000000_DP,0.500000000000000_DP, &
+      &     0.277777777777778_DP,0.444444444444444_DP,0.277777777777778_DP, &
+      &     0.173927422568727_DP,0.326072577431273_DP,0.326072577431273_DP,0.173927422568727_DP ]
+    NULLIFY(interpolatedPoint)
+    NULLIFY(interpolatedPointMetrics)
+    NULLIFY(interpolationParameters)
+
+    CALL Enters("FIeld_CalculateEnclosedVolume",err,error,*999)
+
+    IF(.NOT.ASSOCIATED(field)) THEN
+      CALL FLAG_ERROR("Field is not associated.",err,error,*999)
+    END IF
+    IF(.NOT.field%FIELD_FINISHED) THEN
+      CALL FLAG_ERROR("Field has not been finished.",err,error,*999)
+    END IF
+    IF(.NOT.field%type==FIELD_GEOMETRIC_TYPE) THEN
+      CALL FLAG_ERROR("Field does not have a geometric type, cannot evaluate enclosed volume.",err,error,*999)
+    END IF
+    IF(.NOT.ASSOCIATED(field%GEOMETRIC_FIELD_PARAMETERS)) THEN
+      CALL FLAG_ERROR("Field geometric field parameters are not associated.",err,error,*999)
+    END IF
+
+    NULLIFY(coordinateSystem)
+    CALL FIELD_COORDINATE_SYSTEM_GET(field,coordinateSystem,err,error,*999)
+    CALL FIELD_INTERPOLATION_PARAMETERS_INITIALISE(field,interpolationParameters,err,error,*999)
+    CALL FIELD_INTERPOLATED_POINTS_INITIALISE(interpolationParameters,interpolatedPoint,err,error,*999)
+    CALL FIELD_INTERPOLATED_POINTS_METRICS_INITIALISE(interpolatedPoint,interpolatedPointMetrics,err,error,*999)
+
+    !Loop over elements
+    volume=0.0_DP
+    DO elemIndex=1,SIZE(elems(:))
+      CALL FIELD_INTERPOLATION_PARAMETERS_FACE_GET(FIELD_VALUES_SET_TYPE,faces(elemIndex), &
+        & interpolationParameters(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
+      DO ng=1,numberOfGaussPoints
+        xi=xi_gauss(gaussStart(numberOfGaussPoints)+ng)
+        w=w_gauss(gaussStart(numberOfGaussPoints)+ng)
+        CALL FIELD_INTERPOLATE_XI(FIRST_PART_DERIV,xi,interpolatedPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999)
+        CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(COORDINATE_JACOBIAN_AREA_TYPE, &
+          & interpolatedPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) !Evaluate area at current gauss point.
+        !Evaluate distance from fixed point to global gauss point position.
+        r=2.0_DP
+      ENDDO
+    ENDDO
+
+    CALL Exits("Field_CalculateEnclosedVolume")
+    RETURN
+
+999 CALL Errors("Field_CalculateEnclosedVolume",err,error)
+    CALL Exits("Field_CalculateEnclosedVolume")
+    RETURN 1
+  END SUBROUTINE Field_CalculateEnclosedVolume
+
+  !
+  !================================================================================================================================
+  !
   !>Gets the field label for a field for character labels. \see OPENCMISS::CMISSFieldLabelGet
   SUBROUTINE FIELD_LABEL_GET_C(FIELD,LABEL,ERR,ERROR,*)
 

src/finite_elasticity_routines.f90

@@ -633,7 +633,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT
     REAL(DP), PARAMETER :: TWOTHIRDS_UNITY(6) = [TWOTHIRDS,TWOTHIRDS,TWOTHIRDS,0.0_DP,0.0_DP,0.0_DP] !Rank 2 unit tensor times 2/3 in Voigt form.
     REAL(DP), PARAMETER :: UNITY_DIAGONAL(6)=[1.0_DP,1.0_DP,1.0_DP,0.5_DP,0.5_DP,0.5_DP] !Diagonal of rank 4 unit tensor in Voigt form.
     REAL(DP), POINTER :: C(:) !Parameters for constitutive laws
-    REAL(DP) :: MOD_DZDNU(3,3),MOD_DZDNUT(3,3),AZL(3,3),AZU(3,3)
+    REAL(DP) :: MOD_DZDNU(3,3),MOD_DZDNUT(3,3),AZL(3,3),AZU(3,3), DZDNUT(3,3)
     REAL(DP) :: TRACE,TWOTHIRDS_TRACE
     REAL(DP) :: B(6),E(6),DQ_DE(6)
     TYPE(FIELD_VARIABLE_TYPE), POINTER :: FIELD_VARIABLE
@@ -652,6 +652,8 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT
     MOD_DZDNU=DZDNU*Jznu**(-ONETHIRD)
     CALL MATRIX_TRANSPOSE(MOD_DZDNU,MOD_DZDNUT,ERR,ERROR,*999)
     CALL MATRIX_PRODUCT(MOD_DZDNUT,MOD_DZDNU,AZL,ERR,ERROR,*999)
+    !CALL MATRIX_TRANSPOSE(DZDNU,DZDNUT,ERR,ERROR,*999)
+    !CALL MATRIX_PRODUCT(DZDNUT,DZDNU,AZL,ERR,ERROR,*999)
 
     C=>MATERIALS_INTERPOLATED_POINT%VALUES(:,NO_PART_DERIV)
 
@@ -706,6 +708,8 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT
       ! Do push-forward of 2nd Piola tensor and the material elasticity tensor. 
       CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,MOD_DZDNU,Jznu,ERR,ERROR,*999)
       CALL FINITE_ELASTICITY_PUSH_ELASTICITY_TENSOR(ELASTICITY_TENSOR,MOD_DZDNU,Jznu,ERR,ERROR,*999)
+      !CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,DZDNU,Jznu,ERR,ERROR,*999)
+      !CALL FINITE_ELASTICITY_PUSH_ELASTICITY_TENSOR(ELASTICITY_TENSOR,DZDNU,Jznu,ERR,ERROR,*999)
 
       TRACE=SUM(STRESS_TENSOR(1:3))
       !Calculate isochoric Cauchy tensor (the deviatoric part).
@@ -782,6 +786,8 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT
       !Do push-forward of 2nd Piola tensor and the material elasticity tensor. 
       CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,MOD_DZDNU,Jznu,ERR,ERROR,*999)
       CALL FINITE_ELASTICITY_PUSH_ELASTICITY_TENSOR(ELASTICITY_TENSOR,MOD_DZDNU,Jznu,ERR,ERROR,*999)
+      !CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,DZDNU,Jznu,ERR,ERROR,*999)
+      !CALL FINITE_ELASTICITY_PUSH_ELASTICITY_TENSOR(ELASTICITY_TENSOR,DZDNU,Jznu,ERR,ERROR,*999)
 
       TRACE=SUM(STRESS_TENSOR(1:3))
       !Calculate isochoric Cauchy tensor (the deviatoric part) and volumetric part (hydrostatic pressure).
@@ -4114,7 +4120,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
     REAL(DP) :: ONETHIRD_TRACE
     TYPE(VARYING_STRING) :: LOCAL_ERROR
     TYPE(FIELD_VARIABLE_TYPE), POINTER :: FIELD_VARIABLE
-    REAL(DP) :: MOD_DZDNU(3,3),MOD_DZDNUT(3,3),AZL(3,3),AZU(3,3)
+    REAL(DP) :: MOD_DZDNU(3,3),MOD_DZDNUT(3,3),AZL(3,3),AZU(3,3), DZDNUT(3,3)
     REAL(DP) :: B(6),E(6),DQ_DE(6)
     REAL(DP), POINTER :: C(:) !Parameters for constitutive laws
 
@@ -4128,6 +4134,8 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
     MOD_DZDNU=DZDNU*Jznu**(-1.0_DP/3.0_DP)
     CALL MATRIX_TRANSPOSE(MOD_DZDNU,MOD_DZDNUT,ERR,ERROR,*999)
     CALL MATRIX_PRODUCT(MOD_DZDNUT,MOD_DZDNU,AZL,ERR,ERROR,*999)
+    !CALL MATRIX_TRANSPOSE(DZDNU,DZDNUT,ERR,ERROR,*999)
+    !CALL MATRIX_PRODUCT(DZDNUT,DZDNU,AZL,ERR,ERROR,*999)
     C=>MATERIALS_INTERPOLATED_POINT%VALUES(:,NO_PART_DERIV)
 
     SELECT CASE(EQUATIONS_SET%SUBTYPE)
@@ -4165,6 +4173,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
 
       !Do push-forward of 2nd Piola tensor. 
       CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,MOD_DZDNU,Jznu,ERR,ERROR,*999)
+      !CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,DZDNU,Jznu,ERR,ERROR,*999)
       !Calculate isochoric Cauchy tensor (the deviatoric part) and add the volumetric part (the hydrostatic pressure).
       ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:3))/3.0_DP
       STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)-ONETHIRD_TRACE+P
@@ -4181,7 +4190,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
         !add active contraction stress values
         !the active stress is stored inside the independent field that has been set up in the user program.
         !for generality we could set up 3 components in independent field for 3 different active stress components
-        !!!!! Be aware for modified DZDNU, check if this the right way to do it?
+        !!!! Be aware for modified DZDNU, check if this the right way to do it?
         CALL FIELD_VARIABLE_GET(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE,FIELD_VARIABLE,ERR,ERROR,*999)
         DO component_idx=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS
           dof_idx=FIELD_VARIABLE%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP% &
@@ -4193,6 +4202,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO
       ENDIF
       ! Do push-forward of 2nd Piola tensor. 
       CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,MOD_DZDNU,Jznu,ERR,ERROR,*999)
+      !CALL FINITE_ELASTICITY_PUSH_STRESS_TENSOR(STRESS_TENSOR,DZDNU,Jznu,ERR,ERROR,*999)
       !Calculate isochoric Cauchy tensor (the deviatoric part) and add the volumetric part (the hydrostatic pressure).
       ONETHIRD_TRACE=SUM(STRESS_TENSOR(1:3))/3.0_DP
       STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)-ONETHIRD_TRACE+P
@@ -7799,7 +7809,7 @@ SUBROUTINE FINITE_ELASTICITY_PRE_SOLVE(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
             CASE DEFAULT
               LOCAL_ERROR="Problem subtype "//TRIM(NUMBER_TO_VSTRING(CONTROL_LOOP%PROBLEM%SUBTYPE,"*",ERR,ERROR))// &
                 & " is not valid for a finite elasticity problem class."
-              CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999)
+             CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999)
           END SELECT
         ELSE
           CALL FLAG_ERROR("Problem is not associated.",ERR,ERROR,*999)

src/mesh_routines.f90

@@ -1704,7 +1704,7 @@ SUBROUTINE DECOMP_TOPOLOGY_ELEM_ADJACENT_ELEM_CALCULATE(TOPOLOGY,ERR,ERROR,*)
     INTEGER(INTG) :: j,ne,ne1,nep1,ni,nic,nn,nn1,nn2,nn3,np,np1,DUMMY_ERR,FACE_XI(2),FACE_XIC(3),NODE_POSITION_INDEX(4)
     INTEGER(INTG) :: xi_direction,direction_index,xi_dir_check,xi_dir_search,NUMBER_NODE_MATCHES
     INTEGER(INTG) :: candidate_idx,face_node_idx,node_idx,surrounding_el_idx,candidate_el,idx
-    INTEGER(INTG) :: SURROUNDING_ELEMENTS(100) !Fixed size array... should be large enough
+    INTEGER(INTG) :: SURROUNDING_ELEMENTS(300) !Fixed size array... should be large enough
     INTEGER(INTG) :: NUMBER_SURROUNDING,NUMBER_OF_NODES_XIC(4)
     INTEGER(INTG), ALLOCATABLE :: NODE_MATCHES(:),ADJACENT_ELEMENTS(:)
     LOGICAL :: XI_COLLAPSED,FACE_COLLAPSED(-3:3),SUBSET