From d45d05e64ea49d4591ee53a4f32767d25380f82f Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Wed, 8 Jul 2020 12:46:57 -0600 Subject: [PATCH] Fix bug that could cause incorrect augmented mesh in L2-L2 or L2-P load (#488) Line2-to-Line2 or Line2-to-Point load transfers in some instances could have created augmented meshes with nodes in incorrect places, and possibly cause incorrect loads transfer in the process. Also fixed the mesh debugging code so that if a mesh has scalar quantities, it stores the number of scalars that it is writing to the binary file. --- modules/nwtc-library/src/ModMesh.f90 | 3 +- modules/nwtc-library/src/ModMesh_Mapping.f90 | 80 +++++++++----------- modules/nwtc-library/src/ModMesh_Types.f90 | 4 + 3 files changed, 40 insertions(+), 47 deletions(-) diff --git a/modules/nwtc-library/src/ModMesh.f90 b/modules/nwtc-library/src/ModMesh.f90 index e5af024392..4048677f02 100644 --- a/modules/nwtc-library/src/ModMesh.f90 +++ b/modules/nwtc-library/src/ModMesh.f90 @@ -91,6 +91,7 @@ SUBROUTINE MeshWrBin ( UnIn, M, ErrStat, ErrMsg, FileName) WRITE (UnIn, IOSTAT=ErrStat2) M%fieldmask ! BJJ: do we need to verify that this is size B4Ki? WRITE (UnIn, IOSTAT=ErrStat2) INT(M%Nnodes,B4Ki) WRITE (UnIn, IOSTAT=ErrStat2) INT(M%nelemlist,B4Ki) + if (M%Fieldmask(MASKID_SCALAR)) WRITE (UnIn, IOSTAT=ErrStat2) INT(M%nScalars,B4Ki) !........... @@ -2411,7 +2412,7 @@ SUBROUTINE MeshCommit( Mesh, ErrStat, ErrMess ) Mesh%ElemTable(ELEMENT_LINE2)%Elements(J)%det_jac = 0.5_ReKi * TwoNorm( n1_n2_vector ) ! = L / 2 - IF ( EqualRealNos( 2.0_ReKi*Mesh%ElemTable(ELEMENT_LINE2)%Elements(J)%det_jac, 0.0_Reki ) ) THEN + IF ( 2.0_ReKi*Mesh%ElemTable(ELEMENT_LINE2)%Elements(J)%det_jac < MIN_LINE2_ELEMENT_LENGTH ) THEN ErrStat = ErrID_Fatal ErrMess = trim(ErrMess)//"MeshCommit: Line2 element "//TRIM(Num2Lstr(j))//" has 0 length."//NewLine// & " n2 = n("//TRIM(Num2Lstr(n2))//") = ("//TRIM(Num2Lstr(Mesh%Position(1,n2)))//','//TRIM(Num2Lstr(mesh%position(2,n2)))//','//TRIM(Num2Lstr(mesh%position(3,n2))) //')'//NewLine// & diff --git a/modules/nwtc-library/src/ModMesh_Mapping.f90 b/modules/nwtc-library/src/ModMesh_Mapping.f90 index c2a13e4ee0..de3fdb7d0f 100644 --- a/modules/nwtc-library/src/ModMesh_Mapping.f90 +++ b/modules/nwtc-library/src/ModMesh_Mapping.f90 @@ -79,7 +79,7 @@ MODULE ModMesh_Mapping REAL(R8Ki), ALLOCATABLE :: LoadLn2_M(:,:) !< The 3-components of the moments for each node of an element in the point-to-line load mapping (for each element) TYPE(MeshMapLinearizationType) :: dM !< type that contains information for linearization matrices, partial M partial u (or y) - END TYPE + END TYPE MeshMapType ! note that these parameters must be negative (positive indicates the node/element it is mapped to) INTEGER(IntKi), PARAMETER :: NODE_NOT_MAPPED = -1 !< constant that indicates a node is not mapped @@ -391,7 +391,7 @@ SUBROUTINE MeshMapWrBin( UnIn, Src, Dest, MeshMap, ErrStat, ErrMsg, FileName ) IF ( HasLoadFields(Src) ) THEN WRITE (UnIn, IOSTAT=ErrStat2) INT(1,B4Ki) ! contains a load field - MapMat = 0.0_ReKi + MapMat = 0.0_ReKi IF ( Src%ElemTable(ELEMENT_POINT)%Nelem > 0 ) THEN DO i=1,Src%Nnodes j=MeshMap%MapLoads(i)%OtherMesh_Element @@ -4217,6 +4217,7 @@ SUBROUTINE Create_Augmented_Ln2_Src_Mesh(Src, Dest, MeshMap, Dest_TYPE, ErrStat, INTEGER(IntKi) :: Aug_NElem, curr_Aug_NElem INTEGER(IntKi) :: n1, n2 REAL(ReKi) :: p_ED(3), p_ES(3), n1S_nD_vector(3), position(3) + REAL(ReKi) :: p_ED_orig(3), denom_orig REAL(R8Ki) :: RefOrientation(3,3) REAL(DbKi) :: TmpVec(3), RefOrientationD(3,3), FieldValue(3,2) ! values for interpolating direction cosine matrices REAL(ReKi) :: denom, elem_position @@ -4320,6 +4321,7 @@ SUBROUTINE Create_Augmented_Ln2_Src_Mesh(Src, Dest, MeshMap, Dest_TYPE, ErrStat, IF ( Dest_TYPE == ELEMENT_LINE2 ) THEN p_eD = dest%Position(:, dest%ElemTable(Dest_TYPE)%Elements(jElem)%ElemNodes(2)) & - dest%Position(:, dest%ElemTable(Dest_TYPE)%Elements(jElem)%ElemNodes(1)) + p_eD_orig = p_eD ! save for later calculations (to allow point elements, too) END IF @@ -4337,11 +4339,11 @@ SUBROUTINE Create_Augmented_Ln2_Src_Mesh(Src, Dest, MeshMap, Dest_TYPE, ErrStat, IF ( .NOT. EqualRealNos( denom, 0.0_ReKi) ) THEN ! we ignore source elements that are parallel to the destination element (i.e., denom == 0) - DO jNode = j, NumNodes( Dest_TYPE ) + DO jNode = j, NumNodes( Dest_TYPE ) ! check each node of the destination element n1S_nD_vector = dest%Position(:, dest%ElemTable(Dest_TYPE)%Elements(jElem)%ElemNodes(jNode)) & - Temp_Ln2_Src%Position(:, Temp_Ln2_Src%ElemTable(ELEMENT_LINE2)%Elements(iElem)%ElemNodes(1)) - elem_position = DOT_PRODUCT( p_eD, n1S_nD_vector ) / denom + elem_position = DOT_PRODUCT( p_eD, n1S_nD_vector ) / denom ! Eq 37 (AIAA 2015 paper) !bjj: todo: we need to set this TOL based on actual distances, not relative values (0,1) on an element.... ! for now, I've calculated the element length inside this tolerance and reserve the right to reject new nodes that create 0-length elements. @@ -4356,14 +4358,14 @@ SUBROUTINE Create_Augmented_Ln2_Src_Mesh(Src, Dest, MeshMap, Dest_TYPE, ErrStat, ! calculate the position and orientation relative to the *original* source element: n1=Src%ElemTable(ELEMENT_LINE2)%Elements(Original_Src_Element(iElem))%ElemNodes(1) n2=Src%ElemTable(ELEMENT_LINE2)%Elements(Original_Src_Element(iElem))%ElemNodes(2) - + p_eS = Src%Position(:, n2) - Src%Position(:, n1) - IF ( Dest_TYPE == ELEMENT_POINT ) p_eD = p_eS - denom = DOT_PRODUCT( p_eD , p_eS ) ! we don't need to check that this is zero because it's just a shorter version of the temp Temp_Ln2_Src element + IF ( Dest_TYPE == ELEMENT_POINT ) p_eD_orig = p_eS + denom_orig = DOT_PRODUCT( p_eD_orig , p_eS ) ! we don't need to check that this is zero because it's just a shorter version of the temp Temp_Ln2_Src element n1S_nD_vector = dest%Position(:, dest%ElemTable(Dest_TYPE)%Elements(jElem)%ElemNodes(jNode)) & - Src%Position(:, n1 ) - shape_fn2(Aug_Nnodes) = DOT_PRODUCT( p_eD, n1S_nD_vector ) / denom ! save this for later, when we need to map the mesh fields... + shape_fn2(Aug_Nnodes) = DOT_PRODUCT( p_eD_orig, n1S_nD_vector ) / denom_orig ! save this for later, when we need to map the mesh fields... ! interpolate position on the original souce element: @@ -4372,19 +4374,19 @@ SUBROUTINE Create_Augmented_Ln2_Src_Mesh(Src, Dest, MeshMap, Dest_TYPE, ErrStat, ! let's just verify that this new node (n1) doesn't give us zero-length elements: ! (note we use the NEW (not original) source element, which may have been split) - p_eS = position - Temp_Ln2_Src%Position(:, Temp_Ln2_Src%ElemTable(ELEMENT_LINE2)%Elements(iElem)%ElemNodes(1) ) + p_eS = position - Temp_Ln2_Src%Position(:, Temp_Ln2_Src%ElemTable(ELEMENT_LINE2)%Elements(iElem)%ElemNodes(1) ) L = SQRT(dot_product(p_eS,p_eS)) ! length of new element - - IF ( L < TOL ) THEN ! this element is basically zero length + + IF ( L < MIN_LINE2_ELEMENT_LENGTH ) THEN ! this element is basically zero length ! for numerical reasons, we really didn't want this node.... - Aug_Nnodes = Aug_Nnodes - 1 + Aug_Nnodes = Aug_Nnodes - 1 ELSE ! let's verify the other node (n2) of this element doesn't give zero-length: - p_eS = position - Temp_Ln2_Src%Position(:, Temp_Ln2_Src%ElemTable(ELEMENT_LINE2)%Elements(iElem)%ElemNodes(2)) + p_eS = position - Temp_Ln2_Src%Position(:, Temp_Ln2_Src%ElemTable(ELEMENT_LINE2)%Elements(iElem)%ElemNodes(2)) L = SQRT(dot_product(p_eS,p_eS)) ! length of new element - IF ( L < TOL ) THEN ! this element is basically zero length + IF ( L < MIN_LINE2_ELEMENT_LENGTH ) THEN ! this element is basically zero length ! for numerical reasons, we really didn't want this node.... Aug_Nnodes = Aug_Nnodes - 1 ELSE @@ -4398,45 +4400,31 @@ SUBROUTINE Create_Augmented_Ln2_Src_Mesh(Src, Dest, MeshMap, Dest_TYPE, ErrStat, RETURN END IF - Original_Src_Element( Aug_NElem ) = Original_Src_Element( iElem ) ! this node can now be mapped to original source mesh element + Original_Src_Element( Aug_NElem ) = Original_Src_Element( iElem ) ! this node can now be mapped to original source mesh element - ! get the Reference orientation for this new node -#ifdef __NN_ORIENTATIONS - ! set the RefOrientation based on proximity to original element's nodes: - IF ( NINT( shape_fn2(Aug_Nnodes) ) .EQ. 0 ) THEN - n = n1 - ELSE - n = n2 - END IF - RefOrientation = Src%RefOrientation(:, :, n) -#else - - ! convert DCMs to tensors: + ! get the Reference orientation for this new node + ! convert DCMs to tensors: - RefOrientationD = Src%RefOrientation(:, :, n1) - CALL DCM_logmap( RefOrientationD, FieldValue(:,1), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev) RETURN + RefOrientationD = Src%RefOrientation(:, :, n1) + CALL DCM_logmap( RefOrientationD, FieldValue(:,1), ErrStat, ErrMsg ) + IF (ErrStat >= AbortErrLev) RETURN - RefOrientationD = Src%RefOrientation(:, :, n2) - CALL DCM_logmap( RefOrientationD, FieldValue(:,2), ErrStat, ErrMsg ) - IF (ErrStat >= AbortErrLev) RETURN + RefOrientationD = Src%RefOrientation(:, :, n2) + CALL DCM_logmap( RefOrientationD, FieldValue(:,2), ErrStat, ErrMsg ) + IF (ErrStat >= AbortErrLev) RETURN - CALL DCM_SetLogMapForInterp( FieldValue ) ! make sure we don't cross a 2pi boundary + CALL DCM_SetLogMapForInterp( FieldValue ) ! make sure we don't cross a 2pi boundary - ! interpolate tensors: - TmpVec = (1.0_ReKi - shape_fn2(Aug_Nnodes)) * FieldValue(:, 1) & - + shape_fn2(Aug_Nnodes) * FieldValue(:, 2) + ! interpolate tensors: + TmpVec = (1.0_ReKi - shape_fn2(Aug_Nnodes)) * FieldValue(:, 1) & + + shape_fn2(Aug_Nnodes) * FieldValue(:, 2) - ! convert back to DCM: - RefOrientationD = DCM_exp( TmpVec ) - RefOrientation = REAL(RefOrientationD, R8Ki) + ! convert back to DCM: + RefOrientationD = DCM_exp( TmpVec ) + RefOrientation = REAL(RefOrientationD, R8Ki) - -#endif - - CALL MeshPositionNode ( Mesh = Temp_Ln2_Src & ,INode = Aug_Nnodes & ,Pos = position & @@ -4446,10 +4434,10 @@ SUBROUTINE Create_Augmented_Ln2_Src_Mesh(Src, Dest, MeshMap, Dest_TYPE, ErrStat, CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) IF (ErrStat >= AbortErrLev) THEN CALL CleanUp() - RETURN + RETURN END IF - ! if we have to check a second node, we need to first recalculate p_eS and denom on Temp_Ln2_Src: + ! if we have to check a second node, we need to first recalculate p_eS and denom on Temp_Ln2_Src: IF ( jNode < NumNodes( Dest_TYPE )) THEN j = jNode+1 ! start on the next node CYCLE Src_Elements diff --git a/modules/nwtc-library/src/ModMesh_Types.f90 b/modules/nwtc-library/src/ModMesh_Types.f90 index f4c82f3977..88b3eca620 100644 --- a/modules/nwtc-library/src/ModMesh_Types.f90 +++ b/modules/nwtc-library/src/ModMesh_Types.f90 @@ -71,6 +71,10 @@ MODULE ModMesh_Types LOGICAL, PARAMETER :: mesh_debug = .FALSE. + +! REAL(ReKi), PARAMETER :: MIN_LINE2_ELEMENT_LENGTH = 0.001 ! 1 millimeter + REAL(ReKi), PARAMETER :: MIN_LINE2_ELEMENT_LENGTH = sqrt(epsilon(1.0_ReKi)) ! old length + !> element record type: fields for a particular element TYPE, PUBLIC :: ElemRecType ! note: any fields added to this type must be copied in Mesh_MoveAlloc_ElemRecType (modmesh_types::mesh_movealloc_elemrectype)