diff --git a/src/FieldExport.c b/src/FieldExport.c index d3589531..cf34cdf8 100644 --- a/src/FieldExport.c +++ b/src/FieldExport.c @@ -83,6 +83,8 @@ #define FIELD_IO_INTERPOLATION_HEADER_SCALE 4 #define FIELD_IO_INTERPOLATION_HEADER_NODAL 5 #define FIELD_IO_INTERPOLATION_HEADER_GRID 6 +#define FIELD_IO_INTERPOLATION_HEADER_GAUSS 7 +#define FIELD_IO_INTERPOLATION_HEADER_CONSTANT 8 #define FIELD_GEOMETRIC_TYPE 1 //Geometric field \see FIELD_ROUTINES_FieldTypes,FIELD_ROUTINES @@ -221,6 +223,7 @@ static SessionListEntry *FieldExport_GetSession( const int handle ) return NULL; } +const int FieldExport_InterpolationType( const int interpType ); /* CMISS-formatted file export routines. @@ -291,8 +294,17 @@ static int FieldExport_File_InterpolationHeader( FileSession *const session, con { if( labelType == FIELD_IO_INTERPOLATION_HEADER_GRID ) { - scaleFactorCount *= 2; - label = "l.Lagrange"; + scaleFactorCount *= 2; + label = "l.Lagrange"; + } + else if( labelType == FIELD_IO_INTERPOLATION_HEADER_CONSTANT ) + { + label = "constant"; + } + else if( labelType == FIELD_IO_INTERPOLATION_HEADER_GAUSS ) + { + scaleFactorCount *= 2; + label = "l.Lagrange"; } else { @@ -380,6 +392,12 @@ static int FieldExport_File_InterpolationHeader( FileSession *const session, con case FIELD_IO_INTERPOLATION_HEADER_GRID: FieldExport_FPrintf( session, ", no modify, grid based.\n" ); break; + case FIELD_IO_INTERPOLATION_HEADER_CONSTANT: + FieldExport_FPrintf( session, ", no modify, grid based.\n" ); + break; + case FIELD_IO_INTERPOLATION_HEADER_GAUSS: + FieldExport_FPrintf( session, ", no modify, grid based.\n" ); + break; default: return FIELD_EXPORT_ERROR_UNKNOWN_LABEL_TYPE; } @@ -572,11 +590,12 @@ static int FieldExport_File_Variable( FileSession *const session, const char *va static int FieldExport_File_CoordinateComponent( FileSession *const session, int coordinateSystemType, - const int componentNumber, const int isNodal, const int numberOfXi, const int *const interpolationXi ) + const int componentNumber, const int interpType, const int numberOfXi, const int *const interpolationXi ) { const char * const componentLabel = FieldExport_GetCoordinateComponentLabel( coordinateSystemType, componentNumber ); - const int headerType = isNodal ? FIELD_IO_INTERPOLATION_HEADER_NODAL : FIELD_IO_INTERPOLATION_HEADER_GRID; + int headerType; + headerType = FieldExport_InterpolationType( interpType ); if( componentLabel == NULL ) { FieldExport_FPrintf( session, " %d. ", componentNumber ); @@ -596,9 +615,9 @@ static int FieldExport_File_CoordinateComponent( FileSession *const session, int static int FieldExport_File_Component( FileSession *const session, - const int componentNumber, const int isNodal, const int numberOfXi, const int *const interpolationXi ) + const int componentNumber, const int interpType, const int numberOfXi, const int *const interpolationXi ) { - const int headerType = isNodal ? FIELD_IO_INTERPOLATION_HEADER_NODAL : FIELD_IO_INTERPOLATION_HEADER_GRID; + const int headerType = FieldExport_InterpolationType( interpType ); if( FieldExport_FPrintf( session, " %d. ", componentNumber ) != FIELD_EXPORT_NO_ERROR ) { @@ -609,9 +628,36 @@ static int FieldExport_File_Component( FileSession *const session, } -static int FieldExport_File_ElementGridSize( FileSession *const session, const int numberOfXi ) +static int FieldExport_File_ElementGridSize( FileSession *const session, const int interpType, const int numberOfXi, const int *const numberGauss ) { - int i; + int i,numGrid[3]; + const int headerType = FieldExport_InterpolationType( interpType ); + + if( headerType == FIELD_IO_INTERPOLATION_HEADER_CONSTANT) + { + numGrid[0]=0; + numGrid[1]=0; + numGrid[2]=0; + } + else if( headerType == FIELD_IO_INTERPOLATION_HEADER_GRID ) + { + numGrid[0]=1; + numGrid[1]=1; + numGrid[2]=1; + } + else if( headerType == FIELD_IO_INTERPOLATION_HEADER_GAUSS ) + { + for( i = 0; i < numberOfXi; i++ ) + { + numGrid[i]=numberGauss[i]-1; + } + } + else + { + numGrid[0]=1; + numGrid[1]=1; + numGrid[2]=1; + } if( FieldExport_FPrintf( session, " " ) != FIELD_EXPORT_NO_ERROR ) { @@ -620,7 +666,7 @@ static int FieldExport_File_ElementGridSize( FileSession *const session, const i for( i = 0; i < numberOfXi; i++ ) { - if( FieldExport_FPrintf( session, "#xi%d=1", i+1 ) != FIELD_EXPORT_NO_ERROR ) + if( FieldExport_FPrintf( session, "#xi%d=%d", i+1, numGrid[i] ) != FIELD_EXPORT_NO_ERROR ) { return session->error; } @@ -738,10 +784,10 @@ static int FieldExport_File_ElementNodeScales( FileSession *session, const int i } -static int FieldExport_File_ElementGridValues( FileSession *session, const int isFirstSet, const int dimensionCount, const double value ) +static int FieldExport_File_ElementGridValues( FileSession *session, const int isFirstSet, const int valueCount, const double value ) { int i; - int valueCount = 1 << dimensionCount; + /* int valueCount = 1 << dimensionCount; */ if( isFirstSet ) { @@ -1045,6 +1091,41 @@ static int FieldExport_File_CoordinateDerivativeIndices( FileSession *session, c /* Public API implementation */ + +const int FieldExport_InterpolationType( const int interpType ) +{ + if(interpType == 1) + { + return FIELD_IO_INTERPOLATION_HEADER_CONSTANT; + } + else if(interpType == 2) + { + return FIELD_IO_INTERPOLATION_HEADER_CONSTANT; + } + else if(interpType == 3) + { + return FIELD_IO_INTERPOLATION_HEADER_NODAL; + } + else if(interpType == 4) + { + return FIELD_IO_INTERPOLATION_HEADER_GRID; + } + else if(interpType == 5) + { + return FIELD_IO_INTERPOLATION_HEADER_GAUSS; + } + else if(interpType == 6) + { + return FIELD_IO_INTERPOLATION_HEADER_NODAL; + } + else + { + return FIELD_IO_INTERPOLATION_HEADER_GRID; + } + +} + + int FieldExport_OpenSession( const int type, const char *const name, int * const handle ) { if( type == EXPORT_TYPE_FILE ) @@ -1213,7 +1294,7 @@ int FieldExport_Variable( const int handle, const char *variableName, const int int FieldExport_CoordinateComponent( const int handle, int coordinateSystemType, - const int componentNumber, const int isNodal, const int numberOfXi, const int * const interpolationXi ) + const int componentNumber, const int interpType, const int numberOfXi, const int * const interpolationXi ) { SessionListEntry *session = FieldExport_GetSession( handle ); @@ -1223,7 +1304,7 @@ int FieldExport_CoordinateComponent( const int handle, int coordinateSystemType, } else if( session->type == EXPORT_TYPE_FILE ) { - return FieldExport_File_CoordinateComponent( &session->fileSession, coordinateSystemType, componentNumber, isNodal, numberOfXi, interpolationXi ); + return FieldExport_File_CoordinateComponent( &session->fileSession, coordinateSystemType, componentNumber, interpType, numberOfXi, interpolationXi ); } else { @@ -1232,7 +1313,7 @@ int FieldExport_CoordinateComponent( const int handle, int coordinateSystemType, } -int FieldExport_Component( const int handle, const int componentNumber, const int isNodal, const int numberOfXi, const int * const interpolationXi ) +int FieldExport_Component( const int handle, const int componentNumber, const int interpType, const int numberOfXi, const int * const interpolationXi ) { SessionListEntry *session = FieldExport_GetSession( handle ); @@ -1242,7 +1323,7 @@ int FieldExport_Component( const int handle, const int componentNumber, const in } else if( session->type == EXPORT_TYPE_FILE ) { - return FieldExport_File_Component( &session->fileSession, componentNumber, isNodal, numberOfXi, interpolationXi ); + return FieldExport_File_Component( &session->fileSession, componentNumber, interpType, numberOfXi, interpolationXi ); } else { @@ -1251,7 +1332,7 @@ int FieldExport_Component( const int handle, const int componentNumber, const in } -int FieldExport_ElementGridSize( const int handle, const int numberOfXi ) +int FieldExport_ElementGridSize( const int handle, const int interpType, const int numberOfXi, const int *const numberGauss ) { SessionListEntry *session = FieldExport_GetSession( handle ); @@ -1261,7 +1342,7 @@ int FieldExport_ElementGridSize( const int handle, const int numberOfXi ) } else if( session->type == EXPORT_TYPE_FILE ) { - return FieldExport_File_ElementGridSize( &session->fileSession, numberOfXi ); + return FieldExport_File_ElementGridSize( &session->fileSession, interpType, numberOfXi, numberGauss ); } else { diff --git a/src/basis_routines.f90 b/src/basis_routines.f90 index b0668e9c..63fe535c 100644 --- a/src/basis_routines.f90 +++ b/src/basis_routines.f90 @@ -1462,790 +1462,750 @@ END SUBROUTINE BASIS_INTERPOLATION_XI_SET_PTR !>Creates and initialises a Lagrange-Hermite tensor product basis that has already been allocated BASIS_CREATE_START !> \see BASIS_ROUTINES::BASIS_CREATE_START - SUBROUTINE BASIS_LHTP_BASIS_CREATE(BASIS,ERR,ERROR,*) + SUBROUTINE Basis_LHTPBasisCreate(basis,err,error,*) !Argument variables - TYPE(BASIS_TYPE), POINTER :: BASIS !MAX_NUM_NODES) MAX_NUM_NODES=BASIS%NUMBER_OF_NODES_XIC(ni) - ENDDO !ni + numberOfNodes=numberOfNodes*basis%NUMBER_OF_NODES_XIC(xiIdx) + IF(basis%NUMBER_OF_NODES_XIC(xiIdx)>maximumNumberOfNodes) maximumNumberOfNodes=basis%NUMBER_OF_NODES_XIC(xiIdx) + ENDDO !xiIdx !If a degenerate (collapsed) basis recalculate the number of nodes from the maximum posible number of nodes - IF(BASIS%DEGENERATE) THEN + IF(basis%degenerate) THEN !Calculate the NODE_AT_COLLAPSE array. - ALLOCATE(NODE_AT_COLLAPSE(NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate at collapse",ERR,ERROR,*999) - POSITION=1 - BASIS%NUMBER_OF_NODES=0 + ALLOCATE(nodeAtCollapse(numberOfNodes),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate at collapse",err,error,*999) + position=1 + basis%NUMBER_OF_NODES=0 !Loop over the maximum number of nodes which is currently set for the basis - DO nn=1,NUMBER_OF_NODES - AT_COLLAPSE=.FALSE. - DO ni=1,BASIS%NUMBER_OF_XI - IF(BASIS%COLLAPSED_XI(ni)==BASIS_COLLAPSED_AT_XI0.AND.POSITION(ni)==1.OR. & - & BASIS%COLLAPSED_XI(ni)==BASIS_COLLAPSED_AT_XI1.AND.POSITION(ni)==BASIS%NUMBER_OF_NODES_XIC(ni)) & - & AT_COLLAPSE=.TRUE. - ENDDO !ni - IF(AT_COLLAPSE) THEN - IF(BASIS%NUMBER_OF_COLLAPSED_XI==1) THEN - FIRST_COLLAPSED_POSITION=POSITION(COLLAPSED_XI(1))==1 + DO localNodeIdx=1,numberOfNodes + atCollapse=.FALSE. + DO xiIdx=1,basis%NUMBER_OF_XI + IF(basis%COLLAPSED_XI(xiIdx)==BASIS_COLLAPSED_AT_XI0.AND.position(xiIdx)==1.OR. & + & basis%COLLAPSED_XI(xiIdx)==BASIS_COLLAPSED_AT_XI1.AND.position(xiIdx)==basis%NUMBER_OF_NODES_XIC(xiIdx)) & + & atCollapse=.TRUE. + ENDDO !xiIdx + IF(atCollapse) THEN + IF(basis%NUMBER_OF_COLLAPSED_XI==1) THEN + firstCollapsedPosition=position(collapsedXi(1))==1 ELSE - FIRST_COLLAPSED_POSITION=(POSITION(COLLAPSED_XI(1))==1).AND.(POSITION(COLLAPSED_XI(2))==1) + firstCollapsedPosition=(position(collapsedXi(1))==1).AND.(position(collapsedXi(2))==1) ENDIF - IF(FIRST_COLLAPSED_POSITION) THEN - BASIS%NUMBER_OF_NODES=BASIS%NUMBER_OF_NODES+1 - NODE_AT_COLLAPSE(BASIS%NUMBER_OF_NODES)=.TRUE. + IF(firstCollapsedPosition) THEN + basis%NUMBER_OF_NODES=basis%NUMBER_OF_NODES+1 + nodeAtCollapse(basis%NUMBER_OF_NODES)=.TRUE. ENDIF ELSE - BASIS%NUMBER_OF_NODES=BASIS%NUMBER_OF_NODES+1 - NODE_AT_COLLAPSE(BASIS%NUMBER_OF_NODES)=.FALSE. + basis%NUMBER_OF_NODES=basis%NUMBER_OF_NODES+1 + nodeAtCollapse(basis%NUMBER_OF_NODES)=.FALSE. ENDIF - POSITION(1)=POSITION(1)+1 - DO ni=1,BASIS%NUMBER_OF_XI - IF(POSITION(ni)>BASIS%NUMBER_OF_NODES_XIC(ni)) THEN - POSITION(ni)=1 - POSITION(ni+1)=POSITION(ni+1)+1 + position(1)=position(1)+1 + DO xiIdx=1,basis%NUMBER_OF_XI + IF(position(xiIdx)>basis%NUMBER_OF_NODES_XIC(xiIdx)) THEN + position(xiIdx)=1 + position(xiIdx+1)=position(xiIdx+1)+1 ENDIF - ENDDO !ni - ENDDO !nn - ALLOCATE(BASIS%NODE_AT_COLLAPSE(BASIS%NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate basis node at collapse",ERR,ERROR,*999) - BASIS%NODE_AT_COLLAPSE(1:BASIS%NUMBER_OF_NODES)=NODE_AT_COLLAPSE(1:BASIS%NUMBER_OF_NODES) - DEALLOCATE(NODE_AT_COLLAPSE) + ENDDO !xiIdx + ENDDO !localNodeIdx + CALL MOVE_ALLOC(nodeAtCollapse,basis%NODE_AT_COLLAPSE) ELSE - BASIS%NUMBER_OF_NODES=NUMBER_OF_NODES - ALLOCATE(BASIS%NODE_AT_COLLAPSE(BASIS%NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate basis node at collapse",ERR,ERROR,*999) - BASIS%NODE_AT_COLLAPSE=.FALSE. - COLLAPSED_XI(1)=1 + basis%NUMBER_OF_NODES=numberOfNodes + ALLOCATE(basis%NODE_AT_COLLAPSE(basis%NUMBER_OF_NODES),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate basis node at collapse.",err,error,*999) + basis%NODE_AT_COLLAPSE=.FALSE. + collapsedXi(1)=1 ENDIF - - ALLOCATE(BASIS%NODE_POSITION_INDEX(BASIS%NUMBER_OF_NODES,BASIS%NUMBER_OF_XI_COORDINATES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate NODE_POSITION_INDEX",ERR,ERROR,*999) - SELECT CASE(BASIS%NUMBER_OF_XI_COORDINATES) + + ALLOCATE(basis%NODE_POSITION_INDEX(basis%NUMBER_OF_NODES,basis%NUMBER_OF_XI_COORDINATES),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate basis node position index.",err,error,*999) + SELECT CASE(basis%NUMBER_OF_XI_COORDINATES) CASE(1) - ALLOCATE(BASIS%NODE_POSITION_INDEX_INV(MAX_NUM_NODES,1,1,1),STAT=ERR) + ALLOCATE(basis%NODE_POSITION_INDEX_INV(maximumNumberOfNodes,1,1,1),STAT=err) CASE(2) - ALLOCATE(BASIS%NODE_POSITION_INDEX_INV(MAX_NUM_NODES,MAX_NUM_NODES,1,1),STAT=ERR) + ALLOCATE(basis%NODE_POSITION_INDEX_INV(maximumNumberOfNodes,maximumNumberOfNodes,1,1),STAT=err) CASE(3) - ALLOCATE(BASIS%NODE_POSITION_INDEX_INV(MAX_NUM_NODES,MAX_NUM_NODES,MAX_NUM_NODES,1),STAT=ERR) + ALLOCATE(basis%NODE_POSITION_INDEX_INV(maximumNumberOfNodes,maximumNumberOfNodes,maximumNumberOfNodes,1),STAT=err) CASE DEFAULT - CALL FLAG_ERROR("Invalid number of coordinates",ERR,ERROR,*999) + CALL FlagError("Invalid number of xi coordinates.",err,error,*999) END SELECT - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate NODE_POSITION_INDEX_INV",ERR,ERROR,*999) - BASIS%NODE_POSITION_INDEX_INV=0 - + IF(ERR/=0) CALL FlagError("Could not allocate node position index inverse.",err,error,*999) + basis%NODE_POSITION_INDEX_INV=0 + !Determine the node position index and it's inverse - POSITION=1 - nn=0 - FIRST_COLLAPSED_POSITION=.TRUE. - DO nn1=1,NUMBER_OF_NODES - AT_COLLAPSE=.FALSE. - IF(BASIS%DEGENERATE) THEN - DO ni=1,BASIS%NUMBER_OF_XI - IF(BASIS%COLLAPSED_XI(ni)==BASIS_COLLAPSED_AT_XI0.AND.POSITION(ni)==1.OR. & - & BASIS%COLLAPSED_XI(ni)==BASIS_COLLAPSED_AT_XI1.AND.POSITION(ni)==BASIS%NUMBER_OF_NODES_XIC(ni)) & - & AT_COLLAPSE=.TRUE. - ENDDO !ni - IF(BASIS%NUMBER_OF_COLLAPSED_XI==1) THEN - FIRST_COLLAPSED_POSITION=POSITION(COLLAPSED_XI(1))==1 + position=1 + localNode=0 + firstCollapsedPosition=.TRUE. + DO localNodeIdx1=1,numberOfNodes + atCollapse=.FALSE. + IF(basis%degenerate) THEN + DO xiIdx=1,basis%NUMBER_OF_XI + IF(basis%COLLAPSED_XI(xiIdx)==BASIS_COLLAPSED_AT_XI0.AND.position(xiIdx)==1.OR. & + & basis%COLLAPSED_XI(xiIdx)==BASIS_COLLAPSED_AT_XI1.AND.position(xiIdx)==basis%NUMBER_OF_NODES_XIC(xiIdx)) & + & atCollapse=.TRUE. + ENDDO !xiIdx + IF(basis%NUMBER_OF_COLLAPSED_XI==1) THEN + firstCollapsedPosition=position(collapsedXi(1))==1 ELSE - FIRST_COLLAPSED_POSITION=(POSITION(COLLAPSED_XI(1))==1).AND.(POSITION(COLLAPSED_XI(2))==1) + firstCollapsedPosition=(position(collapsedXi(1))==1).AND.(position(collapsedXi(2))==1) ENDIF ENDIF - PROCESS_NODE=(AT_COLLAPSE.AND.FIRST_COLLAPSED_POSITION).OR.(.NOT.AT_COLLAPSE) - IF(PROCESS_NODE) THEN - nn=nn+1 - BASIS%NODE_POSITION_INDEX(nn,:)=POSITION(1:BASIS%NUMBER_OF_XI) -!!Should the inverse of the node position index be adjusted so that the collapsed positions point to the collapsed nn? + processNode=(atCollapse.AND.firstCollapsedPosition).OR.(.NOT.atCollapse) + IF(processNode) THEN + localNode=localNode+1 + basis%NODE_POSITION_INDEX(localNode,:)=position(1:basis%NUMBER_OF_XI) +!!Should the inverse of the node position index be adjusted so that the collapsed positions point to the collapsed localNode? !!At the moment this is not the case and they are just set to zero. - SELECT CASE(BASIS%NUMBER_OF_XI) + SELECT CASE(basis%NUMBER_OF_XI) CASE(1) - BASIS%NODE_POSITION_INDEX_INV(BASIS%NODE_POSITION_INDEX(nn,1),1,1,1)=nn + basis%NODE_POSITION_INDEX_INV(basis%NODE_POSITION_INDEX(localNode,1),1,1,1)=localNode CASE(2) - BASIS%NODE_POSITION_INDEX_INV(BASIS%NODE_POSITION_INDEX(nn,1),BASIS%NODE_POSITION_INDEX(nn,2),1,1)=nn + basis%NODE_POSITION_INDEX_INV(basis%NODE_POSITION_INDEX(localNode,1),basis%NODE_POSITION_INDEX(localNode,2),1,1)= & + & localNode CASE(3) - BASIS%NODE_POSITION_INDEX_INV(BASIS%NODE_POSITION_INDEX(nn,1),BASIS%NODE_POSITION_INDEX(nn,2), & - & BASIS%NODE_POSITION_INDEX(nn,3),1)=nn + basis%NODE_POSITION_INDEX_INV(basis%NODE_POSITION_INDEX(localNode,1),basis%NODE_POSITION_INDEX(localNode,2), & + & basis%NODE_POSITION_INDEX(localNode,3),1)=localNode CASE DEFAULT - CALL FLAG_ERROR("Invalid number of Xi directions",ERR,ERROR,*999) + CALL FlagError("Invalid number of xi directions.",err,error,*999) END SELECT ENDIF - POSITION(1)=POSITION(1)+1 - DO ni=1,BASIS%NUMBER_OF_XI - IF(POSITION(ni)>BASIS%NUMBER_OF_NODES_XIC(ni)) THEN - POSITION(ni)=1 - POSITION(ni+1)=POSITION(ni+1)+1 + position(1)=position(1)+1 + DO xiIdx=1,basis%NUMBER_OF_XI + IF(position(xiIdx)>basis%NUMBER_OF_NODES_XIC(xiIdx)) THEN + position(xiIdx)=1 + position(xiIdx+1)=position(xiIdx+1)+1 ENDIF - ENDDO !ni - ENDDO !nn1 - !Calculate the maximum number of derivatives and number of element - !parameters - BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES=-1 - BASIS%NUMBER_OF_ELEMENT_PARAMETERS=0 - DO nn=1,BASIS%NUMBER_OF_NODES - NUMBER_OF_DERIVATIVES=1 - DO ni=1,BASIS%NUMBER_OF_XI - IF((.NOT.BASIS%NODE_AT_COLLAPSE(nn).OR.BASIS%COLLAPSED_XI(ni)==BASIS_NOT_COLLAPSED).AND. & - & BASIS%INTERPOLATION_TYPE(ni)==BASIS_HERMITE_INTERPOLATION.AND. & - & (BASIS%INTERPOLATION_ORDER(ni)==BASIS_CUBIC_INTERPOLATION_ORDER.OR. & - & (BASIS%NODE_POSITION_INDEX(nn,ni)==1.AND. & - & BASIS%INTERPOLATION_ORDER(ni)==BASIS_QUADRATIC2_INTERPOLATION_ORDER).OR. & - & (BASIS%NODE_POSITION_INDEX(nn,ni)==2.AND. & - & BASIS%INTERPOLATION_ORDER(ni)==BASIS_QUADRATIC1_INTERPOLATION_ORDER))) THEN + ENDDO !xiIdx + ENDDO !localNodeIdx1 + !Calculate the maximum number of derivatives and the number of element parameters + basis%MAXIMUM_NUMBER_OF_DERIVATIVES=-1 + basis%NUMBER_OF_ELEMENT_PARAMETERS=0 + DO localNodeIdx=1,basis%NUMBER_OF_NODES + numberOfDerivatives=1 + DO xiIdx=1,basis%NUMBER_OF_XI + IF((.NOT.basis%NODE_AT_COLLAPSE(localNodeIdx).OR.basis%COLLAPSED_XI(xiIdx)==BASIS_NOT_COLLAPSED).AND. & + & basis%INTERPOLATION_TYPE(xiIdx)==BASIS_HERMITE_INTERPOLATION.AND. & + & (basis%INTERPOLATION_ORDER(xiIdx)==BASIS_CUBIC_INTERPOLATION_ORDER.OR. & + & (basis%NODE_POSITION_INDEX(localNodeIdx,xiIdx)==1.AND. & + & basis%INTERPOLATION_ORDER(xiIdx)==BASIS_QUADRATIC2_INTERPOLATION_ORDER).OR. & + & (basis%NODE_POSITION_INDEX(localNodeIdx,xiIdx)==2.AND. & + & basis%INTERPOLATION_ORDER(xiIdx)==BASIS_QUADRATIC1_INTERPOLATION_ORDER))) THEN !Derivative in this direction - NUMBER_OF_DERIVATIVES=NUMBER_OF_DERIVATIVES*2 + numberOfDerivatives=numberOfDerivatives*2 ENDIF - ENDDO !ni - BASIS%NUMBER_OF_ELEMENT_PARAMETERS=BASIS%NUMBER_OF_ELEMENT_PARAMETERS+NUMBER_OF_DERIVATIVES - IF(NUMBER_OF_DERIVATIVES>BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES) BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES= & - & NUMBER_OF_DERIVATIVES - ENDDO !nn + ENDDO !xiIdx + basis%NUMBER_OF_ELEMENT_PARAMETERS=basis%NUMBER_OF_ELEMENT_PARAMETERS+numberOfDerivatives + IF(numberOfDerivatives>basis%MAXIMUM_NUMBER_OF_DERIVATIVES) basis%MAXIMUM_NUMBER_OF_DERIVATIVES=numberOfDerivatives + ENDDO !localNodeIdx !Now set up the number of derivatives and derivative order index - ALLOCATE(BASIS%NUMBER_OF_DERIVATIVES(BASIS%NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate NUMBER_OF_DERIVATIVES",ERR,ERROR,*999) - ALLOCATE(BASIS%DERIVATIVE_ORDER_INDEX(BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES,BASIS%NUMBER_OF_NODES, & - & BASIS%NUMBER_OF_XI),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate DERIVATIVE_ORDER_INDEX",ERR,ERROR,*999) - ALLOCATE(BASIS%DERIVATIVE_ORDER_INDEX_INV(FIRST_PART_DERIV,FIRST_PART_DERIV,FIRST_PART_DERIV, & - & BASIS%NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate DERIVATIVE_ORDER_INDEX_INV",ERR,ERROR,*999) - ALLOCATE(BASIS%PARTIAL_DERIVATIVE_INDEX(BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES,BASIS%NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate PARTIAL_DERIVATIVE_INDEX",ERR,ERROR,*999) - ALLOCATE(BASIS%ELEMENT_PARAMETER_INDEX(BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES,BASIS%NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate ELEMENT_PARAMETER_INDEX",ERR,ERROR,*999) - ALLOCATE(BASIS%ELEMENT_PARAMETER_INDEX_INV(2,BASIS%NUMBER_OF_ELEMENT_PARAMETERS),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate ELEMENT_PARAMETER_INDEX_INV",ERR,ERROR,*999) + ALLOCATE(basis%NUMBER_OF_DERIVATIVES(basis%NUMBER_OF_NODES),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate number of derivatives.",err,error,*999) + ALLOCATE(basis%DERIVATIVE_ORDER_INDEX(basis%MAXIMUM_NUMBER_OF_DERIVATIVES,basis%NUMBER_OF_NODES, & + & basis%NUMBER_OF_XI),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate derivative order index.",err,error,*999) + ALLOCATE(basis%DERIVATIVE_ORDER_INDEX_INV(FIRST_PART_DERIV,FIRST_PART_DERIV,FIRST_PART_DERIV, & + & basis%NUMBER_OF_NODES),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate derivative order index inverse.",err,error,*999) + ALLOCATE(basis%PARTIAL_DERIVATIVE_INDEX(basis%MAXIMUM_NUMBER_OF_DERIVATIVES,basis%NUMBER_OF_NODES),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate partial derivative index.",err,error,*999) + ALLOCATE(basis%ELEMENT_PARAMETER_INDEX(basis%MAXIMUM_NUMBER_OF_DERIVATIVES,basis%NUMBER_OF_NODES),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate element parameter index.",err,error,*999) + ALLOCATE(basis%ELEMENT_PARAMETER_INDEX_INV(2,basis%NUMBER_OF_ELEMENT_PARAMETERS),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate element parameter index inverse.",err,error,*999) !Set the derivative order index and its inverse, the element parameter index and the partial derivative index. - ns=0 - BASIS%DERIVATIVE_ORDER_INDEX=0 - BASIS%DERIVATIVE_ORDER_INDEX_INV=0 - DO nn=1,BASIS%NUMBER_OF_NODES - BASIS%NUMBER_OF_DERIVATIVES(nn)=1 - DO ni=1,BASIS%NUMBER_OF_XI - IF((.NOT.BASIS%NODE_AT_COLLAPSE(nn).OR.BASIS%COLLAPSED_XI(ni)==BASIS_NOT_COLLAPSED).AND. & - & BASIS%INTERPOLATION_TYPE(ni)==BASIS_HERMITE_INTERPOLATION.AND. & - & (BASIS%INTERPOLATION_ORDER(ni)==BASIS_CUBIC_INTERPOLATION_ORDER.OR. & - & (BASIS%NODE_POSITION_INDEX(nn,ni)==1.AND. & - & BASIS%INTERPOLATION_ORDER(ni)==BASIS_QUADRATIC2_INTERPOLATION_ORDER).OR. & - & (BASIS%NODE_POSITION_INDEX(nn,ni)==2.AND. & - & BASIS%INTERPOLATION_ORDER(ni)==BASIS_QUADRATIC1_INTERPOLATION_ORDER))) THEN - OLD_NUMBER_OF_DERIVATIVES=BASIS%NUMBER_OF_DERIVATIVES(nn) - BASIS%NUMBER_OF_DERIVATIVES(nn)=BASIS%NUMBER_OF_DERIVATIVES(nn)*2 - DO nk=1,OLD_NUMBER_OF_DERIVATIVES - BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,ni)=NO_PART_DERIV - BASIS%DERIVATIVE_ORDER_INDEX(OLD_NUMBER_OF_DERIVATIVES+nk,nn,ni)=FIRST_PART_DERIV - DO ni2=1,ni-1 - BASIS%DERIVATIVE_ORDER_INDEX(OLD_NUMBER_OF_DERIVATIVES+nk,nn,ni2)=BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,ni2) - ENDDO !ni - ENDDO !nk + elementParameter=0 + basis%DERIVATIVE_ORDER_INDEX=0 + basis%DERIVATIVE_ORDER_INDEX_INV=0 + DO localNodeIdx=1,basis%NUMBER_OF_NODES + basis%NUMBER_OF_DERIVATIVES(localNodeIdx)=1 + DO xiIdx1=1,basis%NUMBER_OF_XI + IF((.NOT.basis%NODE_AT_COLLAPSE(localNodeIdx).OR.basis%COLLAPSED_XI(xiIdx1)==BASIS_NOT_COLLAPSED).AND. & + & basis%INTERPOLATION_TYPE(xiIdx1)==BASIS_HERMITE_INTERPOLATION.AND. & + & (basis%INTERPOLATION_ORDER(xiIdx1)==BASIS_CUBIC_INTERPOLATION_ORDER.OR. & + & (basis%NODE_POSITION_INDEX(localNodeIdx,xiIdx1)==1.AND. & + & basis%INTERPOLATION_ORDER(xiIdx1)==BASIS_QUADRATIC2_INTERPOLATION_ORDER).OR. & + & (basis%NODE_POSITION_INDEX(localNodeIdx,xiIdx1)==2.AND. & + & basis%INTERPOLATION_ORDER(xiIdx1)==BASIS_QUADRATIC1_INTERPOLATION_ORDER))) THEN + oldNumberOfDerivatives=basis%NUMBER_OF_DERIVATIVES(localNodeIdx) + basis%NUMBER_OF_DERIVATIVES(localNodeIdx)=basis%NUMBER_OF_DERIVATIVES(localNodeIdx)*2 + DO derivativeIdx=1,oldNumberOfDerivatives + basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,xiIdx1)=NO_PART_DERIV + basis%DERIVATIVE_ORDER_INDEX(oldNumberOfDerivatives+derivativeIdx,localNodeIdx,xiIdx1)=FIRST_PART_DERIV + DO xiIdx2=1,xiIdx1-1 + basis%DERIVATIVE_ORDER_INDEX(oldNumberOfDerivatives+derivativeIdx,localNodeIdx,xiIdx2)= & + & basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,xiIdx2) + ENDDO !xiIdx2 + ENDDO !derivativeIdx ELSE - DO nk=1,BASIS%NUMBER_OF_DERIVATIVES(nn) - BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,ni)=1 - ENDDO !nk + DO derivativeIdx=1,basis%NUMBER_OF_DERIVATIVES(localNodeIdx) + basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,xiIdx1)=1 + ENDDO !derivativeIdx ENDIF - ENDDO !ni - - DO nk=1,BASIS%NUMBER_OF_DERIVATIVES(nn) - ns=ns+1 - BASIS%ELEMENT_PARAMETER_INDEX(nk,nn)=ns - BASIS%ELEMENT_PARAMETER_INDEX_INV(1,ns)=nn - BASIS%ELEMENT_PARAMETER_INDEX_INV(2,ns)=nk - SELECT CASE(BASIS%NUMBER_OF_XI) + ENDDO !xiIdx1 + + DO derivativeIdx=1,basis%NUMBER_OF_DERIVATIVES(localNodeIdx) + elementParameter=elementParameter+1 + basis%ELEMENT_PARAMETER_INDEX(derivativeIdx,localNodeIdx)=elementParameter + basis%ELEMENT_PARAMETER_INDEX_INV(1,elementParameter)=localNodeIdx + basis%ELEMENT_PARAMETER_INDEX_INV(2,elementParameter)=derivativeIdx + SELECT CASE(basis%NUMBER_OF_XI) CASE(1) - BASIS%DERIVATIVE_ORDER_INDEX_INV(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,1),1,1,nn)=nk - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,1)) + basis%DERIVATIVE_ORDER_INDEX_INV(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,1),1,1,localNodeIdx)= & + & derivativeIdx + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,1)) CASE(NO_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=NO_PART_DERIV + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=NO_PART_DERIV CASE(FIRST_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S1 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S1 CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE(2) - BASIS%DERIVATIVE_ORDER_INDEX_INV(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,1), & - & BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,2),1,nn)=nk - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,1)) + basis%DERIVATIVE_ORDER_INDEX_INV(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,1), & + & basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,2),1,localNodeIdx)=derivativeIdx + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,1)) CASE(NO_PART_DERIV) - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,2)) + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,2)) CASE(NO_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=NO_PART_DERIV + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=NO_PART_DERIV CASE(FIRST_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S2 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S2 CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE(FIRST_PART_DERIV) - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,2)) + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,2)) CASE(NO_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S1 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S1 CASE(FIRST_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S1_S2 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S1_S2 CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE(3) - BASIS%DERIVATIVE_ORDER_INDEX_INV(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,1), & - & BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,2),BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,3),nn)=nk - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,1)) + basis%DERIVATIVE_ORDER_INDEX_INV(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,1), & + & basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,2), & + & basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,3),localNodeIdx)=derivativeIdx + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,1)) CASE(NO_PART_DERIV) - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,2)) + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,2)) CASE(NO_PART_DERIV) - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,3)) + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,3)) CASE(NO_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=NO_PART_DERIV + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=NO_PART_DERIV CASE(FIRST_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S3 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S3 CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE(FIRST_PART_DERIV) - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,3)) + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,3)) CASE(NO_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S2 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S2 CASE(FIRST_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S2_S3 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S2_S3 CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE(FIRST_PART_DERIV) - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,2)) + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,2)) CASE(NO_PART_DERIV) - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,3)) + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,3)) CASE(NO_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S1 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S1 CASE(FIRST_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S1_S3 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S1_S3 CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE(FIRST_PART_DERIV) - SELECT CASE(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn,3)) + SELECT CASE(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx,3)) CASE(NO_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S1_S2 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S1_S2 CASE(FIRST_PART_DERIV) - BASIS%PARTIAL_DERIVATIVE_INDEX(nk,nn)=PART_DERIV_S1_S2_S3 + basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx)=PART_DERIV_S1_S2_S3 CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE DEFAULT - CALL FLAG_ERROR("Invalid derivative order index.",ERR,ERROR,*999) + CALL FlagError("Invalid derivative order index.",err,error,*999) END SELECT CASE DEFAULT - CALL FLAG_ERROR("Invalid number of Xi direcions.",ERR,ERROR,*999) + CALL FlagError("Invalid number of xi direcions.",err,error,*999) END SELECT - ENDDO !nk - ENDDO !nn + ENDDO !derivativeIdx + ENDDO !localNodeIdx !Set up the line information - SELECT CASE(BASIS%NUMBER_OF_XI) + SELECT CASE(basis%NUMBER_OF_XI) CASE(1) !1 xi directions - NUMBER_OF_LOCAL_LINES=1 - BASIS%NUMBER_OF_LOCAL_LINES=1 - ALLOCATE(BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate number of nodes in local line.",ERR,ERROR,*999) - BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=BASIS%NUMBER_OF_NODES_XIC(1) - ALLOCATE(BASIS%LOCAL_LINE_XI_DIRECTION(NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate local line xi direction.",ERR,ERROR,*999) - BASIS%LOCAL_LINE_XI_DIRECTION(1)=1 - ALLOCATE(BASIS%NODE_NUMBERS_IN_LOCAL_LINE(BASIS%NUMBER_OF_NODES_XIC(1),NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node numbers in local line.",ERR,ERROR,*999) - ALLOCATE(BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(BASIS%NUMBER_OF_NODES_XIC(1),NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate derivative numbers in local line.",ERR,ERROR,*999) - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=NO_PART_DERIV - DO nn2=1,BASIS%NUMBER_OF_NODES_XIC(1) - DO nn1=1,BASIS%NUMBER_OF_NODES - IF(BASIS%NODE_POSITION_INDEX(nn1,1)==nn2) THEN - BASIS%NODE_NUMBERS_IN_LOCAL_LINE(nn2,1)=nn1 - DO nk=1,BASIS%NUMBER_OF_DERIVATIVES(nn2) - IF(BASIS%DERIVATIVE_ORDER_INDEX(nk,nn2,1)==FIRST_PART_DERIV) THEN - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(nn2,1)=nk + numberOfLocalLines=1 + basis%NUMBER_OF_LOCAL_LINES=1 + ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_LINE(numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate number of nodes in local line.",err,error,*999) + basis%NUMBER_OF_NODES_IN_LOCAL_LINE(1)=basis%NUMBER_OF_NODES_XIC(1) + ALLOCATE(basis%LOCAL_LINE_XI_DIRECTION(numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate local line xi direction.",err,error,*999) + basis%LOCAL_LINE_XI_DIRECTION(1)=1 + ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_LINE(basis%NUMBER_OF_NODES_XIC(1),numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate node numbers in local line.",err,error,*999) + ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(basis%NUMBER_OF_NODES_XIC(1),numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate derivative numbers in local line.",err,error,*999) + basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=NO_PART_DERIV + DO localNodeIdx2=1,basis%NUMBER_OF_NODES_XIC(1) + DO localNodeIdx1=1,basis%NUMBER_OF_NODES + IF(basis%NODE_POSITION_INDEX(localNodeIdx1,1)==localNodeIdx2) THEN + basis%NODE_NUMBERS_IN_LOCAL_LINE(localNodeIdx2,1)=localNodeIdx1 + DO derivativeIdx=1,basis%NUMBER_OF_DERIVATIVES(localNodeIdx2) + IF(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNodeIdx2,1)==FIRST_PART_DERIV) THEN + basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(localNodeIdx2,1)=derivativeIdx EXIT ENDIF - ENDDO !nk + ENDDO !derivativeIdx EXIT ENDIF - ENDDO !nn2 - ENDDO !nn1 + ENDDO !localNodeIdx2 + ENDDO !localNodeIdx1 CASE(2) !2 xi directions !Determine the maximum node extent of the basis - MAXIMUM_NODE_EXTENT(1)=MAXVAL(BASIS%NODE_POSITION_INDEX(:,1)) - MAXIMUM_NODE_EXTENT(2)=MAXVAL(BASIS%NODE_POSITION_INDEX(:,2)) + maximumNodeExtent(1)=MAXVAL(basis%NODE_POSITION_INDEX(:,1)) + maximumNodeExtent(2)=MAXVAL(basis%NODE_POSITION_INDEX(:,2)) !Allocate and calculate the lines - NUMBER_OF_LOCAL_LINES=4-BASIS%NUMBER_OF_COLLAPSED_XI - ALLOCATE(BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate number of nodes in local line.",ERR,ERROR,*999) - BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE=0 - ALLOCATE(BASIS%LOCAL_LINE_XI_DIRECTION(NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate local line xi direction.",ERR,ERROR,*999) - ALLOCATE(BASIS%NODE_NUMBERS_IN_LOCAL_LINE(MAXVAL(BASIS%NUMBER_OF_NODES_XIC),NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node numbers in local line",ERR,ERROR,*999) - BASIS%NODE_NUMBERS_IN_LOCAL_LINE=0 - ALLOCATE(BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(MAXVAL(BASIS%NUMBER_OF_NODES_XIC),NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate derivative numbers in local line.",ERR,ERROR,*999) - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=NO_PART_DERIV - ALLOCATE(BASIS%LOCAL_XI_NORMAL(NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate local xi normal.",ERR,ERROR,*999) + numberOfLocalLines=4-basis%NUMBER_OF_COLLAPSED_XI + ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_LINE(numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate number of nodes in local line.",err,error,*999) + basis%NUMBER_OF_NODES_IN_LOCAL_LINE=0 + ALLOCATE(basis%LOCAL_LINE_XI_DIRECTION(numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate local line xi direction.",err,error,*999) + ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_LINE(MAXVAL(basis%NUMBER_OF_NODES_XIC),numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate node numbers in local line",err,error,*999) + basis%NODE_NUMBERS_IN_LOCAL_LINE=0 + ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(MAXVAL(basis%NUMBER_OF_NODES_XIC),numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate derivative numbers in local line.",err,error,*999) + basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=NO_PART_DERIV + ALLOCATE(basis%LOCAL_XI_NORMAL(numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate local xi normal.",err,error,*999) !Find the lines - BASIS%NUMBER_OF_LOCAL_LINES=0 - DO ni1=1,2 - ni2=OTHER_XI_DIRECTIONS2(ni1) - !We are looking for lines in the ni1 direction from the direction of ni1=0 - !Loop over the element extremes in the ni2 direction - DO nn2=1,MAXIMUM_NODE_EXTENT(ni2),MAXIMUM_NODE_EXTENT(ni2)-1 - NODE_COUNT=0 - SPECIAL_NODE_COUNT=0 - NODES_IN_LINE=0 - DO nn1=1,BASIS%NUMBER_OF_NODES - IF(BASIS%COLLAPSED_XI(ni1)/=BASIS_NOT_COLLAPSED) THEN - !The current xi direction, ni1, is in a degenerate plane - IF(BASIS%COLLAPSED_XI(ni2)==BASIS_XI_COLLAPSED) THEN + basis%NUMBER_OF_LOCAL_LINES=0 + DO xiIdx1=1,2 + xiIdx2=OTHER_XI_DIRECTIONS2(xiIdx1) + !We are looking for lines in the xiIdx1 direction from the direction of xiIdx1=0 + !Loop over the element extremes in the xiIdx2 direction + DO localNodeIdx2=1,maximumNodeExtent(xiIdx2),maximumNodeExtent(xiIdx2)-1 + nodeCount=0 + specialNodeCount=0 + nodesInLine=0 + DO localNodeIdx1=1,basis%NUMBER_OF_NODES + IF(basis%COLLAPSED_XI(xiIdx1)/=BASIS_NOT_COLLAPSED) THEN + !The current xi direction, xiIdx1, is in a degenerate plane + IF(basis%COLLAPSED_XI(xiIdx2)==BASIS_XI_COLLAPSED) THEN !The other xi direction is collapsed (must be the case) - IF(BASIS%COLLAPSED_XI(ni1)==BASIS_COLLAPSED_AT_XI0) THEN !Collapsed at the xi=0 end - IF(BASIS%NODE_POSITION_INDEX(nn1,ni2)==nn2.OR.BASIS%NODE_POSITION_INDEX(nn1,ni1)==1) THEN - NODE_COUNT=NODE_COUNT+1 - NODES_IN_LINE(NODE_COUNT)=nn1 + IF(basis%COLLAPSED_XI(xiIdx1)==BASIS_COLLAPSED_AT_XI0) THEN !Collapsed at the xi=0 end + IF(basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx2)==localNodeIdx2.OR. & + & basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx1)==1) THEN + nodeCount=nodeCount+1 + nodesInLine(nodeCount)=localNodeIdx1 ENDIF ELSE !Collapsed at the xi=1 end - IF(BASIS%NODE_POSITION_INDEX(nn1,ni2)==nn2) THEN - NODE_COUNT=NODE_COUNT+1 - NODES_IN_LINE(NODE_COUNT)=nn1 - ELSE IF(BASIS%NODE_POSITION_INDEX(nn1,ni1)==MAXIMUM_NODE_EXTENT(ni1)) THEN - IF(ni1<2) THEN !Special case - put the collapsed node at the end of the line - SPECIAL_NODE_COUNT=SPECIAL_NODE_COUNT+1 - NODES_IN_LINE(MAXIMUM_NODE_EXTENT(ni1))=nn1 + IF(basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx2)==localNodeIdx2) THEN + nodeCount=nodeCount+1 + nodesInLine(nodeCount)=localNodeIdx1 + ELSE IF(basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx1)==maximumNodeExtent(xiIdx1)) THEN + IF(xiIdx1<2) THEN !Special case - put the collapsed node at the end of the line + specialNodeCount=specialNodeCount+1 + nodesInLine(maximumNodeExtent(xiIdx1))=localNodeIdx1 ELSE - NODE_COUNT=NODE_COUNT+1 - NODES_IN_LINE(NODE_COUNT)=nn1 + nodeCount=nodeCount+1 + nodesInLine(nodeCount)=localNodeIdx1 ENDIF ENDIF ENDIF ELSE !The current xi direction must be collapsed - IF(BASIS%NODE_POSITION_INDEX(nn1,ni2)==nn2) THEN - NODE_COUNT=NODE_COUNT+1 - NODES_IN_LINE(NODE_COUNT)=nn1 + IF(basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx2)==localNodeIdx2) THEN + nodeCount=nodeCount+1 + nodesInLine(nodeCount)=localNodeIdx1 ENDIF ENDIF ELSE - !The current xi direction, ni1, is not involved in any collapsed (degenerate) planes - IF(BASIS%NODE_POSITION_INDEX(nn1,ni2)==nn2) THEN - NODE_COUNT=NODE_COUNT+1 - NODES_IN_LINE(NODE_COUNT)=nn1 + !The current xi direction, xiIdx1, is not involved in any collapsed (degenerate) planes + IF(basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx2)==localNodeIdx2) THEN + nodeCount=nodeCount+1 + nodesInLine(nodeCount)=localNodeIdx1 ENDIF ENDIF ENDDO !nn1 - IF((NODE_COUNT+SPECIAL_NODE_COUNT)>1) THEN !More than one node so it is a proper line - BASIS%NUMBER_OF_LOCAL_LINES=BASIS%NUMBER_OF_LOCAL_LINES+1 - BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(BASIS%NUMBER_OF_LOCAL_LINES)=NODE_COUNT+SPECIAL_NODE_COUNT - BASIS%NODE_NUMBERS_IN_LOCAL_LINE(1:BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(BASIS%NUMBER_OF_LOCAL_LINES), & - & BASIS%NUMBER_OF_LOCAL_LINES)=NODES_IN_LINE(1:BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(BASIS%NUMBER_OF_LOCAL_LINES)) - DO nnl=1,BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(BASIS%NUMBER_OF_LOCAL_LINES) - DO nk=1,BASIS%NUMBER_OF_DERIVATIVES(BASIS%NODE_NUMBERS_IN_LOCAL_LINE(nnl,BASIS%NUMBER_OF_LOCAL_LINES)) - IF(BASIS%DERIVATIVE_ORDER_INDEX(nk,BASIS%NODE_NUMBERS_IN_LOCAL_LINE(nnl,BASIS%NUMBER_OF_LOCAL_LINES),ni1)== & - & FIRST_PART_DERIV) THEN - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(nnl,BASIS%NUMBER_OF_LOCAL_LINES)=nk + IF((nodeCount+specialNodeCount)>1) THEN !More than one node so it is a proper line + basis%NUMBER_OF_LOCAL_LINES=basis%NUMBER_OF_LOCAL_LINES+1 + basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES)=nodeCount+specialNodeCount + basis%NODE_NUMBERS_IN_LOCAL_LINE(1:basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES), & + & basis%NUMBER_OF_LOCAL_LINES)=nodesInLine(1:basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES)) + DO localLineNodeIdx=1,basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES) + DO derivativeIdx=1,basis%NUMBER_OF_DERIVATIVES(basis%NODE_NUMBERS_IN_LOCAL_LINE(localLineNodeIdx, & + & basis%NUMBER_OF_LOCAL_LINES)) + IF(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,basis%NODE_NUMBERS_IN_LOCAL_LINE(localLineNodeIdx, & + & basis%NUMBER_OF_LOCAL_LINES),xiIdx1)==FIRST_PART_DERIV) THEN + basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(localLineNodeIdx,basis%NUMBER_OF_LOCAL_LINES)=derivativeIdx EXIT ENDIF - ENDDO !nk - ENDDO !nnl - BASIS%LOCAL_LINE_XI_DIRECTION(BASIS%NUMBER_OF_LOCAL_LINES)=ni1 - IF(nn2==1) THEN - BASIS%LOCAL_XI_NORMAL(BASIS%NUMBER_OF_LOCAL_LINES)=-ni2 + ENDDO !derivativeIdx + ENDDO !localLineNodeIdx + basis%LOCAL_LINE_XI_DIRECTION(basis%NUMBER_OF_LOCAL_LINES)=xiIdx1 + IF(localNodeIdx2==1) THEN + basis%LOCAL_XI_NORMAL(basis%NUMBER_OF_LOCAL_LINES)=-xiIdx2 ELSE - BASIS%LOCAL_XI_NORMAL(BASIS%NUMBER_OF_LOCAL_LINES)=ni2 + basis%LOCAL_XI_NORMAL(basis%NUMBER_OF_LOCAL_LINES)=xiIdx2 ENDIF ENDIF - ENDDO !nn2 - ENDDO !ni1 + ENDDO !localNodeIdx2 + ENDDO !localNodeIdx1 CASE(3) !3 xi directions !Determine the maximum node extent of the basis - MAXIMUM_NODE_EXTENT(1)=MAXVAL(BASIS%NODE_POSITION_INDEX(:,1)) - MAXIMUM_NODE_EXTENT(2)=MAXVAL(BASIS%NODE_POSITION_INDEX(:,2)) - MAXIMUM_NODE_EXTENT(3)=MAXVAL(BASIS%NODE_POSITION_INDEX(:,3)) + maximumNodeExtent(1)=MAXVAL(basis%NODE_POSITION_INDEX(:,1)) + maximumNodeExtent(2)=MAXVAL(basis%NODE_POSITION_INDEX(:,2)) + maximumNodeExtent(3)=MAXVAL(basis%NODE_POSITION_INDEX(:,3)) !Allocate and calculate the lines - IF(BASIS%NUMBER_OF_COLLAPSED_XI==1) THEN - NUMBER_OF_LOCAL_LINES=9 - NUMBER_OF_LOCAL_FACES=5 - BASIS%NUMBER_OF_LOCAL_FACES=5 - ELSE IF(BASIS%NUMBER_OF_COLLAPSED_XI==2) THEN - NUMBER_OF_LOCAL_LINES=8 - NUMBER_OF_LOCAL_FACES=5 - BASIS%NUMBER_OF_LOCAL_FACES=5 + IF(basis%NUMBER_OF_COLLAPSED_XI==1) THEN + numberOfLocalLines=9 + numberOfLocalFaces=5 + ELSE IF(basis%NUMBER_OF_COLLAPSED_XI==2) THEN + numberOfLocalLines=8 + numberOfLocalFaces=5 ELSE - NUMBER_OF_LOCAL_LINES=12 - NUMBER_OF_LOCAL_FACES=6 - BASIS%NUMBER_OF_LOCAL_FACES=6 + numberOfLocalLines=12 + numberOfLocalFaces=6 ENDIF + basis%NUMBER_OF_LOCAL_FACES=numberOfLocalFaces + + ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_LINE(numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate number of nodes in local line.",err,error,*999) + basis%NUMBER_OF_NODES_IN_LOCAL_LINE=0 + + ALLOCATE(basis%NUMBER_OF_NODES_IN_LOCAL_FACE(numberOfLocalFaces),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate number of nodes in local face.",err,error,*999) + basis%NUMBER_OF_NODES_IN_LOCAL_FACE=0 + + ALLOCATE(basis%LOCAL_LINE_XI_DIRECTION(numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate local line xi direction.",err,error,*999) + + ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_LINE(MAXVAL(basis%NUMBER_OF_NODES_XIC),numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate node numbers in local line.",err,error,*999) + basis%NODE_NUMBERS_IN_LOCAL_LINE=0 + + ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(MAXVAL(basis%NUMBER_OF_NODES_XIC),numberOfLocalLines),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate derivative numbers in local line.",err,error,*999) + basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=NO_PART_DERIV + + ALLOCATE(basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0:basis%MAXIMUM_NUMBER_OF_DERIVATIVES, & + & MAXVAL(basis%NUMBER_OF_NODES_XIC)**2,numberOfLocalFaces),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate derivative numbers in local face.",err,error,*999) + basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE=NO_PART_DERIV + basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0,:,:)=1 + + ALLOCATE(basis%NODE_NUMBERS_IN_LOCAL_FACE(MAX(maximumNodeExtent(2)*maximumNodeExtent(3), & + & maximumNodeExtent(3)*maximumNodeExtent(1),maximumNodeExtent(2)*maximumNodeExtent(1)), & + & numberOfLocalFaces),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate node numbers in local face.",err,error,*999) + basis%NODE_NUMBERS_IN_LOCAL_FACE=0 + + ALLOCATE(basis%LOCAL_XI_NORMAL(numberOfLocalFaces),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate local xi normal.",err,error,*999) - ALLOCATE(BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate number of nodes in local line",ERR,ERROR,*999) - BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE=0 - - ALLOCATE(BASIS%NUMBER_OF_NODES_IN_LOCAL_FACE(NUMBER_OF_LOCAL_FACES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate number of nodes in local face",ERR,ERROR,*999) - BASIS%NUMBER_OF_NODES_IN_LOCAL_FACE=0 - - ALLOCATE(BASIS%LOCAL_LINE_XI_DIRECTION(NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate local line xi direction",ERR,ERROR,*999) - - ALLOCATE(BASIS%NODE_NUMBERS_IN_LOCAL_LINE(MAXVAL(BASIS%NUMBER_OF_NODES_XIC),NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node numbers in local line",ERR,ERROR,*999) - BASIS%NODE_NUMBERS_IN_LOCAL_LINE=0 - - ALLOCATE(BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(MAXVAL(BASIS%NUMBER_OF_NODES_XIC),NUMBER_OF_LOCAL_LINES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate derivative numbers in local line",ERR,ERROR,*999) - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_LINE=NO_PART_DERIV + ALLOCATE(basis%LOCAL_FACE_XI_DIRECTION(numberOfLocalFaces),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate local face xi direction.",err,error,*999) - ALLOCATE(BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES, & - & MAXVAL(BASIS%NUMBER_OF_NODES_XIC)**2,NUMBER_OF_LOCAL_FACES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate derivative numbers in local face",ERR,ERROR,*999) - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE=NO_PART_DERIV - - ALLOCATE(BASIS%NODE_NUMBERS_IN_LOCAL_FACE(MAX(MAXIMUM_NODE_EXTENT(2)*MAXIMUM_NODE_EXTENT(3), & - & MAXIMUM_NODE_EXTENT(3)*MAXIMUM_NODE_EXTENT(1),MAXIMUM_NODE_EXTENT(2)*MAXIMUM_NODE_EXTENT(1)), & - & NUMBER_OF_LOCAL_FACES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node numbers in local face",ERR,ERROR,*999) - BASIS%NODE_NUMBERS_IN_LOCAL_FACE=0 - ALLOCATE(BASIS%LOCAL_XI_NORMAL(NUMBER_OF_LOCAL_FACES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate local xi normal.",ERR,ERROR,*999) - - ALLOCATE(BASIS%LOCAL_FACE_XI_DIRECTION(NUMBER_OF_LOCAL_FACES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate local face xi direction",ERR,ERROR,*999) - !Find the lines and faces - BASIS%NUMBER_OF_LOCAL_LINES=0 - DO ni1=1,3 - ni2=OTHER_XI_DIRECTIONS3(ni1,2,1) - ni3=OTHER_XI_DIRECTIONS3(ni1,3,1) - !We are looking for lines going in the ni1 direction, starting from ni1=0. - DO nn3=1,MAXIMUM_NODE_EXTENT(ni3),MAXIMUM_NODE_EXTENT(ni3)-1 - DO nn2=1,MAXIMUM_NODE_EXTENT(ni2),MAXIMUM_NODE_EXTENT(ni2)-1 - NODE_COUNT=0 - SPECIAL_NODE_COUNT=0 - NODES_IN_LINE=0 + basis%NUMBER_OF_LOCAL_LINES=0 + DO xiIdx1=1,3 + xiIdx2=OTHER_XI_DIRECTIONS3(xiIdx1,2,1) + xiIdx3=OTHER_XI_DIRECTIONS3(xiIdx1,3,1) + !We are looking for lines going in the xiIdx1 direction, starting from xiIdx1=0. + DO localNodeIdx3=1,maximumNodeExtent(xiIdx3),maximumNodeExtent(xiIdx3)-1 + DO localNodeIdx2=1,maximumNodeExtent(xiIdx2),maximumNodeExtent(xiIdx2)-1 + nodeCount=0 + specialNodeCount=0 + nodesInLine=0 !Iterate over nodes in the line of interest - DO nn1=1,BASIS%NUMBER_OF_NODES - IF(BASIS%COLLAPSED_XI(ni1)/=BASIS_NOT_COLLAPSED) THEN - !The current xi direction, ni1, is involved in a collapsed (degenerate) plane - IF(BASIS%COLLAPSED_XI(ni2)==BASIS_XI_COLLAPSED.AND.BASIS%COLLAPSED_XI(ni3)==BASIS_XI_COLLAPSED) THEN + DO localNodeIdx1=1,basis%NUMBER_OF_NODES + IF(basis%COLLAPSED_XI(xiIdx1)/=BASIS_NOT_COLLAPSED) THEN + !The current xi direction, xiIdx1, is involved in a collapsed (degenerate) plane + IF(basis%COLLAPSED_XI(xiIdx2)==BASIS_XI_COLLAPSED.AND.basis%COLLAPSED_XI(xiIdx3)==BASIS_XI_COLLAPSED) THEN !Both of the other two xi directions are collapsed - IF(BASIS%COLLAPSED_XI(ni1)==BASIS_COLLAPSED_AT_XI0) THEN !Collapsed at the xi=0 end - IF((BASIS%NODE_POSITION_INDEX(nn1,ni2)==nn2.OR.BASIS%NODE_POSITION_INDEX(nn1,ni1)==1).AND. & - & (BASIS%NODE_POSITION_INDEX(nn1,ni3)==nn3.OR.BASIS%NODE_POSITION_INDEX(nn1,ni1)==1)) THEN - NODE_COUNT=NODE_COUNT+1 - NODES_IN_LINE(NODE_COUNT)=nn1 + IF(basis%COLLAPSED_XI(xiIdx1)==BASIS_COLLAPSED_AT_XI0) THEN !Collapsed at the xi=0 end + IF((basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx2)==localNodeIdx2.OR. & + & basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx1)==1).AND. & + & (basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx3)==localNodeIdx3.OR. & + & basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx1)==1)) THEN + nodeCount=nodeCount+1 + nodesInLine(nodeCount)=localNodeIdx1 ENDIF ELSE !Collapsed at the xi=1 end - IF(BASIS%NODE_POSITION_INDEX(nn1,ni2)==nn2.AND.BASIS%NODE_POSITION_INDEX(nn1,ni3)==nn3) THEN - NODE_COUNT=NODE_COUNT+1 - NODES_IN_LINE(NODE_COUNT)=nn1 - ELSE IF(BASIS%NODE_POSITION_INDEX(nn1,ni1)==MAXIMUM_NODE_EXTENT(ni1)) THEN - IF(ni1<3) THEN !Special case - put the collapsed node at the end of the line - SPECIAL_NODE_COUNT=SPECIAL_NODE_COUNT+1 - NODES_IN_LINE(MAXIMUM_NODE_EXTENT(ni1))=nn1 + IF(basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx2)==localNodeIdx2.AND. & + & basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx3)==localNodeIdx3) THEN + nodeCount=nodeCount+1 + nodesInLine(nodeCount)=localNodeIdx1 + ELSE IF(basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx1)==maximumNodeExtent(xiIdx1)) THEN + IF(xiIdx1<3) THEN !Special case - put the collapsed node at the end of the line + specialNodeCount=specialNodeCount+1 + nodesInLine(maximumNodeExtent(xiIdx1))=localNodeIdx1 ELSE - NODE_COUNT=NODE_COUNT+1 - NODES_IN_LINE(NODE_COUNT)=nn1 + nodeCount=nodeCount+1 + nodesInLine(nodeCount)=localNodeIdx1 ENDIF ENDIF ENDIF ELSE - IF(BASIS%COLLAPSED_XI(ni2)==BASIS_XI_COLLAPSED) THEN - !The other ni2 xi direction is collapsed - IF(BASIS%COLLAPSED_XI(ni1)==BASIS_COLLAPSED_AT_XI0) THEN !Collapsed at the xi=0 end - IF((BASIS%NODE_POSITION_INDEX(nn1,ni2)==nn2.OR.BASIS%NODE_POSITION_INDEX(nn1,ni1)==1).AND. & - & BASIS%NODE_POSITION_INDEX(nn1,ni3)==nn3) THEN - NODE_COUNT=NODE_COUNT+1 - NODES_IN_LINE(NODE_COUNT)=nn1 + IF(basis%COLLAPSED_XI(xiIdx2)==BASIS_XI_COLLAPSED) THEN + !The other xiIdx2 xi direction is collapsed + IF(basis%COLLAPSED_XI(xiIdx1)==BASIS_COLLAPSED_AT_XI0) THEN !Collapsed at the xi=0 end + IF((basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx2)==localNodeIdx2.OR. & + & basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx1)==1).AND. & + & basis%NODE_POSITION_INDEX(localNodeIdx1,xiIdx3)==localNodeIdx3) THEN + nodeCount=nodeCount+1 + nodesInLine(nodeCount)=localNodeIdx1 ENDIF - ELSE IF(BASIS%COLLAPSED_XI(ni1)==BASIS_COLLAPSED_AT_XI1) THEN !Collapsed at the xi=1 end - IF(BASIS%NODE_POSITION_INDEX(nn1,ni2)==nn2.AND.BASIS%NODE_POSITION_INDEX(nn1,ni3)==nn3) THEN - NODE_COUNT=NODE_COUNT+1 - NODES_IN_LINE(NODE_COUNT)=nn1 - ELSE IF(BASIS%NODE_POSITION_INDEX(nn1,ni1)==MAXIMUM_NODE_EXTENT(ni1).AND. & - & BASIS%NODE_POSITION_INDEX(nn1,ni3)==nn3) THEN - IF(ni11) THEN !More than one node so it is a proper line - BASIS%NUMBER_OF_LOCAL_LINES=BASIS%NUMBER_OF_LOCAL_LINES+1 - BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(BASIS%NUMBER_OF_LOCAL_LINES)=NODE_COUNT+SPECIAL_NODE_COUNT - BASIS%NODE_NUMBERS_IN_LOCAL_LINE(1:BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(BASIS%NUMBER_OF_LOCAL_LINES), & - & BASIS%NUMBER_OF_LOCAL_LINES)=NODES_IN_LINE(1:BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE( & - & BASIS%NUMBER_OF_LOCAL_LINES)) - DO nnl=1,BASIS%NUMBER_OF_NODES_IN_LOCAL_LINE(BASIS%NUMBER_OF_LOCAL_LINES) - DO nk=1,BASIS%NUMBER_OF_DERIVATIVES(BASIS%NODE_NUMBERS_IN_LOCAL_LINE(nnl,BASIS%NUMBER_OF_LOCAL_LINES)) - IF(BASIS%DERIVATIVE_ORDER_INDEX(nk,BASIS%NODE_NUMBERS_IN_LOCAL_LINE(nnl,BASIS%NUMBER_OF_LOCAL_LINES),ni1)== & - & FIRST_PART_DERIV) THEN - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(nnl,BASIS%NUMBER_OF_LOCAL_LINES)=nk + ENDDO !localNodeIdx1 + IF((nodeCount+specialNodeCount)>1) THEN !More than one node so it is a proper line + basis%NUMBER_OF_LOCAL_LINES=basis%NUMBER_OF_LOCAL_LINES+1 + basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES)=nodeCount+specialNodeCount + basis%NODE_NUMBERS_IN_LOCAL_LINE(1:basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES), & + & basis%NUMBER_OF_LOCAL_LINES)=nodesInLine(1:basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES)) + DO localLineNodeIdx=1,basis%NUMBER_OF_NODES_IN_LOCAL_LINE(basis%NUMBER_OF_LOCAL_LINES) + DO derivativeIdx=1,basis%NUMBER_OF_DERIVATIVES(basis%NODE_NUMBERS_IN_LOCAL_LINE(localLineNodeIdx, & + & basis%NUMBER_OF_LOCAL_LINES)) + IF(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,basis%NODE_NUMBERS_IN_LOCAL_LINE(localLineNodeIdx, & + & basis%NUMBER_OF_LOCAL_LINES),xiIdx1)==FIRST_PART_DERIV) THEN + basis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(localLineNodeIdx,basis%NUMBER_OF_LOCAL_LINES)=derivativeIdx EXIT ENDIF - ENDDO !nk - ENDDO !nnl - BASIS%LOCAL_LINE_XI_DIRECTION(BASIS%NUMBER_OF_LOCAL_LINES)=ni1 + ENDDO !derivativeIdx + ENDDO !localLineNodeIdx + basis%LOCAL_LINE_XI_DIRECTION(basis%NUMBER_OF_LOCAL_LINES)=xiIdx1 ENDIF - ENDDO !nn2 - ENDDO !nn3 - ENDDO !ni1 - - !Find the local nodes in each face and the local face xi direction - !\todo code below needs to be tested - !Find the faces - nn4=1 - ef=0 !element face counter - DO ni1=1,3 - !Go through all three directions - ni2=OTHER_XI_DIRECTIONS3(ni1,2,1) - ni3=OTHER_XI_DIRECTIONS3(ni1,3,1) - !Pointers for argument list - np1=>nn1 - np2=>nn2 - np3=>nn3 - np4=>nn4 - IF (ni1==1) THEN - ARGLIST=INTG_POINTER(np1,np2,np3,np4) - ELSE IF (ni1==2) THEN - ARGLIST=INTG_POINTER(np2,np1,np3,np4) - ELSE - ARGLIST=INTG_POINTER(np2,np3,np1,np4) - ENDIF - - nn1=MAXIMUM_NODE_EXTENT(ni1) - !start with eventually upper face of ni1 - LOCAL_NODE_COUNT=0 - - IF(BASIS%COLLAPSED_XI(ni1)/=BASIS_COLLAPSED_AT_XI1) THEN - ef=ef+1 - DO nn3=1,MAXIMUM_NODE_EXTENT(ni2) - DO nn2=1,MAXIMUM_NODE_EXTENT(ni3) - LOCAL_NODE_COUNT=LOCAL_NODE_COUNT+1 - BASIS%NODE_NUMBERS_IN_LOCAL_FACE(LOCAL_NODE_COUNT,ef)= & - & BASIS%NODE_POSITION_INDEX_INV(ARGLIST%a1,ARGLIST%a2,ARGLIST%a3,ARGLIST%a4) - IF (BASIS%NODE_NUMBERS_IN_LOCAL_FACE(LOCAL_NODE_COUNT,ef)==0) THEN - IF (BASIS%COLLAPSED_XI(ni1)==BASIS_XI_COLLAPSED) THEN - !set Arglist(ni1-direction)=1 - nn1=1 - BASIS%NODE_NUMBERS_IN_LOCAL_FACE(LOCAL_NODE_COUNT,ef)= & - & BASIS%NODE_POSITION_INDEX_INV(ARGLIST%a1,ARGLIST%a2,ARGLIST%a3,ARGLIST%a4) - nn1=MAXIMUM_NODE_EXTENT(ni1) - ELSE - LOCAL_NODE_COUNT=LOCAL_NODE_COUNT-1 - ENDIF - ENDIF - ENDDO - ENDDO - BASIS%NUMBER_OF_NODES_IN_LOCAL_FACE(ef)=LOCAL_NODE_COUNT - BASIS%LOCAL_FACE_XI_DIRECTION(ef)=ni1 - ENDIF - - nn1=1 - LOCAL_NODE_COUNT=0 - IF(BASIS%COLLAPSED_XI(ni1)/=BASIS_COLLAPSED_AT_XI0) THEN - ef=ef+1 - DO nn3=1,MAXIMUM_NODE_EXTENT(ni2) - DO nn2=1,MAXIMUM_NODE_EXTENT(ni3) - LOCAL_NODE_COUNT=LOCAL_NODE_COUNT+1 - BASIS%NODE_NUMBERS_IN_LOCAL_FACE(LOCAL_NODE_COUNT,ef)= & - & BASIS%NODE_POSITION_INDEX_INV(ARGLIST%a1,ARGLIST%a2,ARGLIST%a3,ARGLIST%a4) - IF (BASIS%NODE_NUMBERS_IN_LOCAL_FACE(LOCAL_NODE_COUNT,ef)==0) THEN - IF (BASIS%COLLAPSED_XI(ni1)==BASIS_XI_COLLAPSED) THEN - !set Arglist(ni1-direction)=Max extent - nn1=MAXIMUM_NODE_EXTENT(ni1) - BASIS%NODE_NUMBERS_IN_LOCAL_FACE(LOCAL_NODE_COUNT,ef)= & - & BASIS%NODE_POSITION_INDEX_INV(ARGLIST%a1,ARGLIST%a2,ARGLIST%a3,ARGLIST%a4) - nn1=1 - ELSE - LOCAL_NODE_COUNT=LOCAL_NODE_COUNT-1 - ENDIF - ENDIF - ENDDO - ENDDO - BASIS%NUMBER_OF_NODES_IN_LOCAL_FACE(ef)=LOCAL_NODE_COUNT - BASIS%LOCAL_FACE_XI_DIRECTION(ef)=-ni1 - ENDIF - ENDDO - - !For each face local derivative index set its corresponding global derivative index - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(1,:,:)=NO_GLOBAL_DERIV - IF(BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES>1) THEN - !Face 1 (xi1+) - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(2,:,1)=GLOBAL_DERIV_S2 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(3,:,1)=GLOBAL_DERIV_S3 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(4,:,1)=GLOBAL_DERIV_S2_S3 - !Face 2 (xi1-) - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(2,:,2)=GLOBAL_DERIV_S2 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(3,:,2)=GLOBAL_DERIV_S3 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(4,:,2)=GLOBAL_DERIV_S2_S3 - !Face 3 (xi2+) - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(2,:,3)=GLOBAL_DERIV_S1 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(3,:,3)=GLOBAL_DERIV_S3 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(4,:,3)=GLOBAL_DERIV_S1_S3 - !Face 4 (xi2-) - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(2,:,4)=GLOBAL_DERIV_S1 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(3,:,4)=GLOBAL_DERIV_S3 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(4,:,4)=GLOBAL_DERIV_S1_S3 - !Face 5 (xi3-) - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(2,:,5)=GLOBAL_DERIV_S1 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(3,:,5)=GLOBAL_DERIV_S2 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(4,:,5)=GLOBAL_DERIV_S1_S2 - !Face 6 (xi3-) - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(2,:,6)=GLOBAL_DERIV_S1 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(3,:,6)=GLOBAL_DERIV_S2 - BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(4,:,6)=GLOBAL_DERIV_S1_S2 - ENDIF - - CASE DEFAULT - CALL FLAG_ERROR("Invalid number of xi directions.",ERR,ERROR,*999) + ENDDO !localNodeIdx2 + ENDDO !localNodeIdx3 + ENDDO !xiIdx1 + + !Find the local nodes and derivatives in each face and the local face xi direction + localFaceIdx=0 + !Loop over the -'ve and +'ve xi direction + DO directionIdx=-1,1,2 + !Loop over the three xi directions + DO xiIdx1=1,3 + !xiIdx1 is the +/- face normal direction. xiIdx2 and xiIdx3 are the xi directions in the face. + xiIdx2=OTHER_XI_DIRECTIONS3(xiIdx1,2,1) + xiIdx3=OTHER_XI_DIRECTIONS3(xiIdx1,3,1) + IF(directionIdx==1) THEN + !The +'ve xi direction + localNodeIdx1=maximumNodeExtent(xiIdx1) + !Compute if the face in the +xiIdx1 direction is collapsed. + collapsedFace=basis%COLLAPSED_XI(xiIdx1)==BASIS_COLLAPSED_AT_XI1 + ELSE + !The -'ve xi direction + localNodeIdx1=1 + !Compute if the face in the +xiIdx1 direction is collapsed. + collapsedFace=basis%COLLAPSED_XI(xiIdx1)==BASIS_COLLAPSED_AT_XI0 + ENDIF + localNodeCount=0 + IF(.NOT.collapsedFace) THEN + !If the face has not been collapsed + localFaceIdx=localFaceIdx+1 + !Loop over the local nodes in the face + DO localNodeIdx3=1,maximumNodeExtent(xiIdx2) + DO localNodeIdx2=1,maximumNodeExtent(xiIdx3) + IF(xiIdx1==1) THEN + localNodeIdx=basis%NODE_POSITION_INDEX_INV(localNodeIdx1,localNodeIdx2,localNodeIdx3,1) + ELSE IF(xiIdx1==2) THEN + localNodeIdx=basis%NODE_POSITION_INDEX_INV(localNodeIdx2,localNodeIdx1,localNodeIdx3,1) + ELSE + localNodeIdx=basis%NODE_POSITION_INDEX_INV(localNodeIdx2,localNodeIdx3,localNodeIdx1,1) + ENDIF + IF(localNodeIdx/=0) THEN + !The node hasn't been collapsed + localNodeCount=localNodeCount+1 + basis%NODE_NUMBERS_IN_LOCAL_FACE(localNodeCount,localFaceIdx)=localNodeIdx + ENDIF + ENDDO !localNodeIdx3 + ENDDO !localNodexIdx2 + basis%NUMBER_OF_NODES_IN_LOCAL_FACE(localFaceIdx)=localNodeCount + basis%LOCAL_FACE_XI_DIRECTION(localFaceIdx)=directionIdx*xiIdx1 + !Compute derivatives in the face + DO localNodeIdx=1,basis%NUMBER_OF_NODES_IN_LOCAL_FACE(localFaceIdx) + localNode=basis%NODE_NUMBERS_IN_LOCAL_FACE(localNodeIdx,localFaceIdx) + localFaceDerivative=0 + DO derivativeIdx=1,basis%NUMBER_OF_DERIVATIVES(localNode) + IF(basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNode,xiIdx2)==FIRST_PART_DERIV.OR. & + & basis%DERIVATIVE_ORDER_INDEX(derivativeIdx,localNode,xiIdx3)==FIRST_PART_DERIV) THEN + localFaceDerivative=localFaceDerivative+1 + basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(localFaceDerivative,localNodeIdx,localFaceIdx)=derivativeIdx + ENDIF + ENDDO !derivativeIdx + basis%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0,localNodeIdx,localFaceIdx)=localFaceDerivative + ENDDO !localNodeIdx + ENDIF + ENDDO !xiIdx1 + ENDDO !directionIdx + + CASE DEFAULT + CALL FlagError("Invalid number of xi directions.",err,error,*999) END SELECT - CALL BASIS_QUADRATURE_CREATE(BASIS,ERR,ERROR,*999) - + CALL BASIS_QUADRATURE_CREATE(basis,err,error,*999) + ELSE - CALL FLAG_ERROR("Basis is not a Lagrange Hermite tensor product basis.",ERR,ERROR,*999) + CALL FlagError("Basis is not a Lagrange Hermite tensor product basis.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Basis is not associated",ERR,ERROR,*999) + CALL FlagError("Basis is not associated.",err,error,*999) ENDIF - - CALL EXITS("BASIS_LHTP_BASIS_CREATE") + CALL Exits("Basis_LHTPBasisCreate") RETURN -999 CALL ERRORS("BASIS_LHTP_BASIS_CREATE",ERR,ERROR) - CALL EXITS("BASIS_LHTP_BASIS_CREATE") +999 IF(ALLOCATED(nodeAtCollapse)) DEALLOCATE(nodeAtCollapse) + CALL Errors("Basis_LHTPBasisCreate",err,error) + CALL Exits("Basis_LHTPBasisCreate") RETURN 1 - END SUBROUTINE BASIS_LHTP_BASIS_CREATE + END SUBROUTINE Basis_LHTPBasisCreate ! ! @@ -2403,7 +2363,7 @@ SUBROUTINE BASIS_LHTP_FAMILY_CREATE(BASIS,ERR,ERROR,*) IF(ASSOCIATED(BASIS)) THEN !Create the main (parent) basis - CALL BASIS_LHTP_BASIS_CREATE(BASIS,ERR,ERROR,*999) + CALL Basis_LHTPBasisCreate(basis,err,error,*999) IF(BASIS%NUMBER_OF_XI>1) THEN !Create the line bases as sub-basis types ALLOCATE(BASIS%LINE_BASES(BASIS%NUMBER_OF_XI),STAT=ERR) @@ -2424,7 +2384,7 @@ SUBROUTINE BASIS_LHTP_FAMILY_CREATE(BASIS,ERR,ERROR,*) !Create the new sub-basis CALL BASIS_SUB_BASIS_CREATE(BASIS,1,[ni],NEW_SUB_BASIS,ERR,ERROR,*999) !Fill in the basis information - CALL BASIS_LHTP_BASIS_CREATE(NEW_SUB_BASIS,ERR,ERROR,*999) + CALL Basis_LHTPBasisCreate(NEW_SUB_BASIS,err,error,*999) BASIS%LINE_BASES(ni)%PTR=>NEW_SUB_BASIS ENDIF ENDDO !ni @@ -2459,7 +2419,7 @@ SUBROUTINE BASIS_LHTP_FAMILY_CREATE(BASIS,ERR,ERROR,*) !Create the new sub-basis CALL BASIS_SUB_BASIS_CREATE(BASIS,2,[FACE_XI(1),FACE_XI(2)],NEW_SUB_BASIS,ERR,ERROR,*999) !Fill in the basis information - CALL BASIS_LHTP_BASIS_CREATE(NEW_SUB_BASIS,ERR,ERROR,*999) + CALL Basis_LHTPBasisCreate(NEW_SUB_BASIS,err,error,*999) NEW_SUB_BASIS%LINE_BASES(1)%PTR=>BASIS%LINE_BASES(FACE_XI(1))%PTR NEW_SUB_BASIS%LINE_BASES(2)%PTR=>BASIS%LINE_BASES(FACE_XI(2))%PTR BASIS%FACE_BASES(ni)%PTR=>NEW_SUB_BASIS @@ -3108,7 +3068,7 @@ SUBROUTINE BASIS_QUADRATURE_CREATE(BASIS,ERR,ERROR,*) DO ng=1,SCHEME%NUMBER_OF_GAUSS CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Gauss point = ",ng,ERR,ERROR,*999) CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,BASIS%NUMBER_OF_XI_COORDINATES,3,3,SCHEME%GAUSS_POSITIONS(:,ng), & - & '(" POSITION(ni) :",3(X,F12.4))','(26X,3(X,F12.4))',ERR,ERROR,*999) + & '(" position(ni) :",3(X,F12.4))','(26X,3(X,F12.4))',ERR,ERROR,*999) CALL WRITE_STRING_FMT_VALUE(DIAGNOSTIC_OUTPUT_TYPE," WEIGHT : ",SCHEME%GAUSS_WEIGHTS(ng), & & "(F12.4)",ERR,ERROR,*999) ENDDO !ng @@ -4652,26 +4612,28 @@ SUBROUTINE BASIS_SIMPLEX_BASIS_CREATE(BASIS,ERR,ERROR,*) ALLOCATE(BASIS%NODE_NUMBERS_IN_LOCAL_FACE(3,BASIS%NUMBER_OF_LOCAL_FACES),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node numbers in local face.",ERR,ERROR,*999) !\todo Number of local face node derivatives currenlty set to 1 (calculation of BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE for simplex elements has not been implemented yet) - ALLOCATE(BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(1,3,BASIS%NUMBER_OF_LOCAL_FACES),STAT=ERR) + ALLOCATE(BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0:1,3,BASIS%NUMBER_OF_LOCAL_FACES),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate derivative numbers in local face.",ERR,ERROR,*999) CASE(BASIS_QUADRATIC_INTERPOLATION_ORDER) ALLOCATE(BASIS%NODE_NUMBERS_IN_LOCAL_FACE(6,BASIS%NUMBER_OF_LOCAL_FACES),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node numbers in local face.",ERR,ERROR,*999) !\todo Number of local face node derivatives currenlty set to 1 (calculation of BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE for simplex elements has not been implemented yet) - ALLOCATE(BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(1,6,BASIS%NUMBER_OF_LOCAL_FACES),STAT=ERR) + ALLOCATE(BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0:1,6,BASIS%NUMBER_OF_LOCAL_FACES),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate derivative numbers in local face.",ERR,ERROR,*999) CASE(BASIS_CUBIC_INTERPOLATION_ORDER) ALLOCATE(BASIS%NODE_NUMBERS_IN_LOCAL_FACE(10,BASIS%NUMBER_OF_LOCAL_FACES),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node numbers in local face.",ERR,ERROR,*999) !\todo Number of local face node derivatives currenlty set to 1 (calculation of BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE for simplex elements has not been implemented yet) - ALLOCATE(BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(1,10,BASIS%NUMBER_OF_LOCAL_FACES),STAT=ERR) + ALLOCATE(BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0:1,10,BASIS%NUMBER_OF_LOCAL_FACES),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate derivative numbers in local face.",ERR,ERROR,*999) CASE DEFAULT LOCAL_ERROR="Interpolation order "//TRIM(NUMBER_TO_VSTRING(BASIS%INTERPOLATION_ORDER(1),"*",ERR,ERROR))// & & " is invalid for a simplex basis type." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) END SELECT - + BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(1,:,:)=NO_PART_DERIV + BASIS%DERIVATIVE_NUMBERS_IN_LOCAL_FACE(0,:,:)=1 + ALLOCATE(BASIS%LOCAL_XI_NORMAL(BASIS%NUMBER_OF_LOCAL_LINES),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate local line normal.",ERR,ERROR,*999) diff --git a/src/data_projection_routines.f90 b/src/data_projection_routines.f90 index 68332378..c165ff98 100644 --- a/src/data_projection_routines.f90 +++ b/src/data_projection_routines.f90 @@ -2845,7 +2845,7 @@ SUBROUTINE DataProjection_ProjectionCandidatesSet(dataProjection,elementUserNumb INTEGER(INTG), INTENT(OUT) :: err !field variable component type pointer for IO TYPE MESH_ELEMENTS_TYPE_PTR_TYPE - TYPE(MeshComponentElementsType), POINTER :: PTR !< pointer field variable component + TYPE(MeshElementsType), POINTER :: PTR !< pointer field variable component END TYPE MESH_ELEMENTS_TYPE_PTR_TYPE !>field variable component type pointer for IO @@ -217,7 +217,7 @@ FUNCTION FieldExport_Variable( handle, variableName, variableNumber, fieldType, INTEGER(C_INT) :: FieldExport_Variable END FUNCTION FieldExport_Variable - FUNCTION FieldExport_CoordinateComponent( handle, coordinateSystemType, componentNumber, isNodal, & + FUNCTION FieldExport_CoordinateComponent( handle, coordinateSystemType, componentNumber, interpType, & & numberOfXi, interpolationXi ) & & BIND(C,NAME="FieldExport_CoordinateComponent") USE TYPES @@ -225,30 +225,32 @@ FUNCTION FieldExport_CoordinateComponent( handle, coordinateSystemType, componen INTEGER(C_INT), VALUE :: handle INTEGER(C_INT), VALUE :: coordinateSystemType INTEGER(C_INT), VALUE :: componentNumber - INTEGER(C_INT), VALUE :: isNodal + INTEGER(C_INT), VALUE :: interpType INTEGER(C_INT), VALUE :: numberOfXi TYPE(C_PTR), VALUE :: interpolationXi INTEGER(C_INT) :: FieldExport_CoordinateComponent END FUNCTION FieldExport_CoordinateComponent - FUNCTION FieldExport_Component( handle, componentNumber, isNodal, numberOfXi, interpolationXi ) & + FUNCTION FieldExport_Component( handle, componentNumber, interpType, numberOfXi, interpolationXi ) & & BIND(C,NAME="FieldExport_Component") USE TYPES USE ISO_C_BINDING INTEGER(C_INT), VALUE :: handle INTEGER(C_INT), VALUE :: componentNumber - INTEGER(C_INT), VALUE :: isNodal + INTEGER(C_INT), VALUE :: interpType INTEGER(C_INT), VALUE :: numberOfXi TYPE(C_PTR), VALUE :: interpolationXi INTEGER(C_INT) :: FieldExport_Component END FUNCTION FieldExport_Component - FUNCTION FieldExport_ElementGridSize( handle, numberOfXi ) & + FUNCTION FieldExport_ElementGridSize( handle, headerType, numberOfXi, numberGauss ) & & BIND(C,NAME="FieldExport_ElementGridSize") USE TYPES USE ISO_C_BINDING INTEGER(C_INT), VALUE :: handle + INTEGER(C_INT), VALUE :: headerType INTEGER(C_INT), VALUE :: numberOfXi + TYPE(C_PTR), VALUE :: numberGauss INTEGER(C_INT) :: FieldExport_ElementGridSize END FUNCTION FieldExport_ElementGridSize @@ -2562,7 +2564,7 @@ SUBROUTINE FIELD_IO_EXPORT_ELEMENTAL_GROUP_HEADER_FORTRAN( global_number, MAX_NO INTEGER(INTG), ALLOCATABLE :: GROUP_NODE(:), GROUP_VARIABLES(:) INTEGER(C_INT), TARGET :: INTERPOLATION_XI(3),ELEMENT_DERIVATIVES(64*64),NUMBER_OF_DERIVATIVES(64), NODE_INDEXES(128) INTEGER(INTG) :: nn, nx, ny, nz, NodesX, NodesY, NodesZ, mm, NUM_OF_VARIABLES, MAX_NUM_NODES !NUM_OF_NODES - INTEGER(INTG) :: local_number, isNodal, NODE_NUMBER, NODE_NUMBER_COUNTER, NODE_NUMBER_COLLAPSED, NUMBER_OF_ELEMENT_NODES + INTEGER(INTG) :: local_number, interpType, NODE_NUMBER, NODE_NUMBER_COUNTER, NODE_NUMBER_COLLAPSED, NUMBER_OF_ELEMENT_NODES INTEGER(INTG) :: num_scl, num_node, comp_idx, scaleIndex, scaleIndex1, var_idx, derivativeIndex !value_idx field_idx global_var_idx comp_idx1 ny2 LOGICAL :: SAME_SCALING_SET @@ -2757,40 +2759,54 @@ SUBROUTINE FIELD_IO_EXPORT_ELEMENTAL_GROUP_HEADER_FORTRAN( global_number, MAX_NO DOMAIN_ELEMENTS=>componentDomain%TOPOLOGY%ELEMENTS BASIS=>DOMAIN_ELEMENTS%ELEMENTS(GROUP_LOCAL_NUMBER(comp_idx))%BASIS - IF( component%INTERPOLATION_TYPE == FIELD_NODE_BASED_INTERPOLATION ) THEN - isNodal = 1 + SELECT CASE(component%INTERPOLATION_TYPE) + CASE(FIELD_CONSTANT_INTERPOLATION) + interpType = 1 + CASE(FIELD_ELEMENT_BASED_INTERPOLATION) + interpType = 2 + CASE(FIELD_NODE_BASED_INTERPOLATION) + interpType = 3 + CASE(FIELD_GRID_POINT_BASED_INTERPOLATION) + interpType = 4 + CASE(FIELD_GAUSS_POINT_BASED_INTERPOLATION) + interpType = 5 + CASE(FIELD_DATA_POINT_BASED_INTERPOLATION) + interpType = 6 + CASE DEFAULT + interpType = 0 + END SELECT + + IF(component%INTERPOLATION_TYPE==FIELD_GAUSS_POINT_BASED_INTERPOLATION) THEN + !TEMP HACK. Fake gauss point export as regular grid. Use interpolation xi to pass in number of Gauss. + INTERPOLATION_XI(1:BASIS%NUMBER_OF_XI)=BASIS%QUADRATURE%NUMBER_OF_GAUSS_XI(1:BASIS%NUMBER_OF_XI) ELSE - isNodal = 0 + !Copy interpolation xi to a temporary array that has the target attribute. gcc bug 38813 prevents using C_LOC with + !the array directly. nb using a fixed length array here which is dangerous but should suffice for now. + INTERPOLATION_XI(1:BASIS%NUMBER_OF_XI)=BASIS%INTERPOLATION_XI(1:BASIS%NUMBER_OF_XI) ENDIF - + IF( variable_ptr%FIELD%TYPE == FIELD_GEOMETRIC_TYPE .AND. & & variable_ptr%VARIABLE_TYPE == FIELD_U_VARIABLE_TYPE ) THEN !!TEMP !ERR = FieldExport_CoordinateComponent( sessionHandle, variable_ptr%FIELD%REGION%COORDINATE_SYSTEM, & ! & component%COMPONENT_NUMBER, basis%NUMBER_OF_XI, C_LOC( basis%INTERPOLATION_XI ) ) - !!Copy interpolation xi to a temporary array that has the target attribute. gcc bug 38813 prevents using C_LOC with - !!the array directly. nb using a fixed length array here which is dangerous but should suffice for now. - INTERPOLATION_XI(1:BASIS%NUMBER_OF_XI)=BASIS%INTERPOLATION_XI(1:BASIS%NUMBER_OF_XI) NULLIFY(COORDINATE_SYSTEM) CALL FIELD_COORDINATE_SYSTEM_GET(variable_ptr%FIELD,COORDINATE_SYSTEM,ERR,ERROR,*999) ERR = FieldExport_CoordinateComponent( sessionHandle, COORDINATE_SYSTEM%TYPE, & - & component%COMPONENT_NUMBER,isNodal,basis%NUMBER_OF_XI, C_LOC( INTERPOLATION_XI )) + & component%COMPONENT_NUMBER,interpType,basis%NUMBER_OF_XI, C_LOC( INTERPOLATION_XI )) ELSE !!TEMP !ERR = FieldExport_Component( sessionHandle, & ! & component%COMPONENT_NUMBER, basis%NUMBER_OF_XI, C_LOC( basis%INTERPOLATION_XI ) ) - !!Copy interpolation xi to a temporary array that has the target attribute. gcc bug 38813 prevents using C_LOC with - !!the array directly. nb using a fixed length array here which is dangerous but should suffice for now. - INTERPOLATION_XI(1:BASIS%NUMBER_OF_XI)=BASIS%INTERPOLATION_XI(1:BASIS%NUMBER_OF_XI) ERR = FieldExport_Component( sessionHandle, & - & component%COMPONENT_NUMBER,isNodal,basis%NUMBER_OF_XI, C_LOC( INTERPOLATION_XI ) ) + & component%COMPONENT_NUMBER,interpType,basis%NUMBER_OF_XI, C_LOC( INTERPOLATION_XI ) ) ENDIF IF(ERR/=0) THEN CALL FLAG_ERROR( "File write error during field export", ERR, ERROR,*999 ) ENDIF - IF( isNodal == 0 ) THEN - ERR = FieldExport_ElementGridSize( sessionHandle, basis%NUMBER_OF_XI ) + IF( interpType /= 3 .AND. interpType /= 6) THEN + ERR = FieldExport_ElementGridSize( sessionHandle, interpType, basis%NUMBER_OF_XI, C_LOC( INTERPOLATION_XI ) ) ELSE ! IF(.NOT.BASIS%DEGENERATE) THEN IF(LIST_COMP_SCALE(comp_idx)==1) THEN @@ -2817,9 +2833,21 @@ SUBROUTINE FIELD_IO_EXPORT_ELEMENTAL_GROUP_HEADER_FORTRAN( global_number, MAX_NO !!$ &ERR,ERROR,*999) NODE_NUMBER_COUNTER=0 - NodesX=BASIS%INTERPOLATION_XI(1)+1 - NodesY=BASIS%INTERPOLATION_XI(2)+1 - NodesZ=BASIS%INTERPOLATION_XI(3)+1 + IF(BASIS%INTERPOLATION_XI(1)>3) THEN + NodesX=2 + ELSE + NodesX=BASIS%INTERPOLATION_XI(1)+1 + ENDIF + IF(BASIS%INTERPOLATION_XI(2)>3) THEN + NodesY=2 + ELSE + NodesY=BASIS%INTERPOLATION_XI(2)+1 + ENDIF + IF(BASIS%INTERPOLATION_XI(3)>3) THEN + NodesZ=2 + ELSE + NodesZ=BASIS%INTERPOLATION_XI(3)+1 + ENDIF !The following if-sentences goes through all possible wedge formed elements and renumber the nodes in order to @@ -3453,7 +3481,7 @@ SUBROUTINE FIELD_IO_EXPORT_ELEMENTS_INTO_LOCAL_FILE(ELEMENTAL_INFO_SET, NAME, my NULLIFY(GEOMETRIC_PARAMETERS) CALL FIELD_PARAMETER_SET_DATA_GET(component%FIELD_VARIABLE%FIELD,& & component%FIELD_VARIABLE%VARIABLE_TYPE,FIELD_VALUES_SET_TYPE,GEOMETRIC_PARAMETERS,ERR,ERROR,*999) - ERR = FieldExport_ElementGridValues( sessionHandle, isFirstValueSet, BASIS%NUMBER_OF_XI, & + ERR = FieldExport_ElementGridValues( sessionHandle, isFirstValueSet, 2**BASIS%NUMBER_OF_XI, & & GEOMETRIC_PARAMETERS(component%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(local_number))) ELSE IF(component%FIELD_VARIABLE%DATA_TYPE==FIELD_INTG_TYPE) THEN NULLIFY(GEOMETRIC_PARAMETERS_INTG) @@ -3463,7 +3491,7 @@ SUBROUTINE FIELD_IO_EXPORT_ELEMENTS_INTO_LOCAL_FILE(ELEMENTAL_INFO_SET, NAME, my IF(ERR/=0) CALL FLAG_ERROR("Could not allocate geometric parameters dp", ERR, ERROR,*999 ) GEOMETRIC_PARAMETERS_DP(1:SIZE(GEOMETRIC_PARAMETERS_INTG))= & & REAL(GEOMETRIC_PARAMETERS_INTG(1:SIZE(GEOMETRIC_PARAMETERS_INTG))) - ERR = FieldExport_ElementGridValues( sessionHandle, isFirstValueSet, BASIS%NUMBER_OF_XI, & + ERR = FieldExport_ElementGridValues( sessionHandle, isFirstValueSet, 2**BASIS%NUMBER_OF_XI, & & GEOMETRIC_PARAMETERS_DP(component%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(local_number))) DEALLOCATE(GEOMETRIC_PARAMETERS_DP) ELSE @@ -3475,7 +3503,7 @@ SUBROUTINE FIELD_IO_EXPORT_ELEMENTS_INTO_LOCAL_FILE(ELEMENTAL_INFO_SET, NAME, my NULLIFY(GEOMETRIC_PARAMETERS) CALL FIELD_PARAMETER_SET_DATA_GET(COMPONENT%FIELD_VARIABLE%FIELD,COMPONENT%FIELD_VARIABLE%VARIABLE_TYPE, & & FIELD_VALUES_SET_TYPE,GEOMETRIC_PARAMETERS,ERR,ERROR,*999) - ERR = FieldExport_ElementGridValues( sessionHandle, isFirstValueSet, BASIS%NUMBER_OF_XI, & + ERR = FieldExport_ElementGridValues( sessionHandle, isFirstValueSet, 2**BASIS%NUMBER_OF_XI, & & GEOMETRIC_PARAMETERS(component%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP)) ELSE IF(component%FIELD_VARIABLE%DATA_TYPE==FIELD_INTG_TYPE) THEN NULLIFY(GEOMETRIC_PARAMETERS_INTG) @@ -3485,13 +3513,37 @@ SUBROUTINE FIELD_IO_EXPORT_ELEMENTS_INTO_LOCAL_FILE(ELEMENTAL_INFO_SET, NAME, my IF(ERR/=0) CALL FLAG_ERROR("Could not allocate geometric parameters dp", ERR, ERROR,*999 ) GEOMETRIC_PARAMETERS_DP(1:SIZE(GEOMETRIC_PARAMETERS_INTG))= & & REAL(GEOMETRIC_PARAMETERS_INTG(1:SIZE(GEOMETRIC_PARAMETERS_INTG))) - ERR = FieldExport_ElementGridValues( sessionHandle, isFirstValueSet, BASIS%NUMBER_OF_XI, & + ERR = FieldExport_ElementGridValues( sessionHandle, isFirstValueSet, 2**BASIS%NUMBER_OF_XI, & & GEOMETRIC_PARAMETERS_DP(component%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP)) DEALLOCATE(GEOMETRIC_PARAMETERS_DP) ELSE CALL FLAG_ERROR( "Only INTG and REAL data types implemented.", ERR, ERROR,*999 ) ENDIF isFirstValueSet = 0 + ELSE IF( component%INTERPOLATION_TYPE == FIELD_GAUSS_POINT_BASED_INTERPOLATION) THEN + IF(component%FIELD_VARIABLE%DATA_TYPE==FIELD_DP_TYPE) THEN + NULLIFY(GEOMETRIC_PARAMETERS) + CALL FIELD_PARAMETER_SET_DATA_GET(COMPONENT%FIELD_VARIABLE%FIELD,COMPONENT%FIELD_VARIABLE%VARIABLE_TYPE, & + & FIELD_VALUES_SET_TYPE,GEOMETRIC_PARAMETERS,ERR,ERROR,*999) + ERR = FieldExport_ElementGridValues( sessionHandle, isFirstValueSet, BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP( & + & BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR%NUMBER_OF_GAUSS, & + & GEOMETRIC_PARAMETERS(component%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(1,local_number))) + ELSE IF(component%FIELD_VARIABLE%DATA_TYPE==FIELD_INTG_TYPE) THEN + NULLIFY(GEOMETRIC_PARAMETERS_INTG) + CALL FIELD_PARAMETER_SET_DATA_GET(COMPONENT%FIELD_VARIABLE%FIELD,COMPONENT%FIELD_VARIABLE%VARIABLE_TYPE, & + & FIELD_VALUES_SET_TYPE,GEOMETRIC_PARAMETERS_INTG,ERR,ERROR,*999) + ALLOCATE(GEOMETRIC_PARAMETERS_DP(SIZE(GEOMETRIC_PARAMETERS_INTG))) + IF(ERR/=0) CALL FLAG_ERROR("Could not allocate geometric parameters dp", ERR, ERROR,*999 ) + GEOMETRIC_PARAMETERS_DP(1:SIZE(GEOMETRIC_PARAMETERS_INTG))= & + & REAL(GEOMETRIC_PARAMETERS_INTG(1:SIZE(GEOMETRIC_PARAMETERS_INTG))) + ERR = FieldExport_ElementGridValues( sessionHandle, isFirstValueSet, BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP( & + & BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR%NUMBER_OF_GAUSS, & + & GEOMETRIC_PARAMETERS_DP(component%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(1,local_number))) + DEALLOCATE(GEOMETRIC_PARAMETERS_DP) + ELSE + CALL FLAG_ERROR( "Only INTG and REAL data types implemented.", ERR, ERROR,*999 ) + ENDIF + isFirstValueSet = 0 ENDIF IF(ERR/=0) THEN @@ -4187,6 +4239,9 @@ SUBROUTINE FIELD_IO_COMPARE_INFO_SET_DERIVATIVES( SET1, SET2, my_computational_n doesMatch = .FALSE. EXIT !out of loop-component_idx=1,SET1%NUMBER_OF_COMPONENTS ENDIF + ELSE + doesMatch = .FALSE. + EXIT ENDIF ! Check that the nodes have the same number of versions, otherwise they must be grouped separately diff --git a/src/field_routines.f90 b/src/field_routines.f90 index 476d09f3..e0cc9c7b 100644 --- a/src/field_routines.f90 +++ b/src/field_routines.f90 @@ -1249,7 +1249,7 @@ SUBROUTINE Field_componentDofGetUserDataPoint(field,variableType,userDataPointNu decompositionTopology=>decomposition%TOPOLOGY userDataPointExists=.TRUE. IF(ASSOCIATED(decompositionTopology)) THEN - CALL DecompositionTopology_DataPointCheckExists(decompositionTopology,userDataPointNumber,userDataPointExists, & + CALL DecompositionTopologyDataPointCheckExists(decompositionTopology,userDataPointNumber,userDataPointExists, & & decompositionLocalDataPointNumber,ghostDataPoint,err,error,*999) IF(userDataPointExists) THEN localDof=fieldVariable%COMPONENTS(componentNumber)%PARAM_TO_DOF_MAP% & @@ -18549,7 +18549,7 @@ SUBROUTINE Field_ParameterSetGetDataPointIntg(field,variableType,fieldSetType,us IF(ASSOCIATED(decomposition)) THEN decompositionTopology=>decomposition%TOPOLOGY IF(ASSOCIATED(decompositionTopology)) THEN - CALL DecompositionTopology_DataPointCheckExists(decompositionTopology,userDataPointNumber, & + CALL DecompositionTopologyDataPointCheckExists(decompositionTopology,userDataPointNumber, & & userDataPointExists,decompositionLocalDataPointNumber,ghostDataPoint,err,error,*999) IF(userDataPointExists) THEN DofIdx=fieldVariable%COMPONENTS(componentNumber)%PARAM_TO_DOF_MAP% & @@ -18701,7 +18701,7 @@ SUBROUTINE Field_ParameterSetGetDataPointSP(field,variableType,fieldSetType,user IF(ASSOCIATED(decomposition)) THEN decompositionTopology=>decomposition%TOPOLOGY IF(ASSOCIATED(decompositionTopology)) THEN - CALL DecompositionTopology_DataPointCheckExists(decompositionTopology,userDataPointNumber, & + CALL DecompositionTopologyDataPointCheckExists(decompositionTopology,userDataPointNumber, & & userDataPointExists,decompositionLocalDataPointNumber,ghostDataPoint,err,error,*999) IF(userDataPointExists) THEN DofIdx=fieldVariable%COMPONENTS(componentNumber)%PARAM_TO_DOF_MAP% & @@ -18853,7 +18853,7 @@ SUBROUTINE Field_ParameterSetGetDataPointDP(field,variableType,fieldSetType,user IF(ASSOCIATED(decomposition)) THEN decompositionTopology=>decomposition%TOPOLOGY IF(ASSOCIATED(decompositionTopology)) THEN - CALL DecompositionTopology_DataPointCheckExists(decompositionTopology,userDataPointNumber, & + CALL DecompositionTopologyDataPointCheckExists(decompositionTopology,userDataPointNumber, & & userDataPointExists,decompositionLocalDataPointNumber,ghostDataPoint,err,error,*999) IF(userDataPointExists) THEN DofIdx=fieldVariable%COMPONENTS(componentNumber)%PARAM_TO_DOF_MAP% & @@ -19005,7 +19005,7 @@ SUBROUTINE Field_ParameterSetGetDataPointL(field,variableType,fieldSetType,userD IF(ASSOCIATED(decomposition)) THEN decompositionTopology=>decomposition%TOPOLOGY IF(ASSOCIATED(decompositionTopology)) THEN - CALL DecompositionTopology_DataPointCheckExists(decompositionTopology,userDataPointNumber, & + CALL DecompositionTopologyDataPointCheckExists(decompositionTopology,userDataPointNumber, & & userDataPointExists,decompositionLocalDataPointNumber,ghostDataPoint,err,error,*999) IF(userDataPointExists) THEN DofIdx=fieldVariable%COMPONENTS(componentNumber)%PARAM_TO_DOF_MAP% & @@ -23224,7 +23224,7 @@ SUBROUTINE Field_ParameterSetUpdateDataPointIntg(field,variableType,fieldSetType IF(ASSOCIATED(decomposition)) THEN decompositionTopology=>decomposition%TOPOLOGY IF(ASSOCIATED(decompositionTopology)) THEN - CALL DecompositionTopology_DataPointCheckExists(decompositionTopology,userDataPointNumber, & + CALL DecompositionTopologyDataPointCheckExists(decompositionTopology,userDataPointNumber, & & userDataPointExists,decompositionLocalDataPointNumber,ghostDataPoint,err,error,*999) IF(userDataPointExists) THEN DofIdx=fieldVariable%COMPONENTS(componentNumber)%PARAM_TO_DOF_MAP% & @@ -23377,7 +23377,7 @@ SUBROUTINE Field_ParameterSetUpdateDataPointSP(field,variableType,fieldSetType,u IF(ASSOCIATED(decomposition)) THEN decompositionTopology=>decomposition%TOPOLOGY IF(ASSOCIATED(decompositionTopology)) THEN - CALL DecompositionTopology_DataPointCheckExists(decompositionTopology,userDataPointNumber, & + CALL DecompositionTopologyDataPointCheckExists(decompositionTopology,userDataPointNumber, & & userDataPointExists,decompositionLocalDataPointNumber,ghostDataPoint,err,error,*999) IF(userDataPointExists) THEN DofIdx=fieldVariable%COMPONENTS(componentNumber)%PARAM_TO_DOF_MAP% & @@ -23530,7 +23530,7 @@ SUBROUTINE Field_ParameterSetUpdateDataPointDP(field,variableType,fieldSetType,u IF(ASSOCIATED(decomposition)) THEN decompositionTopology=>decomposition%TOPOLOGY IF(ASSOCIATED(decompositionTopology)) THEN - CALL DecompositionTopology_DataPointCheckExists(decompositionTopology,userDataPointNumber, & + CALL DecompositionTopologyDataPointCheckExists(decompositionTopology,userDataPointNumber, & & userDataPointExists,decompositionLocalDataPointNumber,ghostDataPoint,err,error,*999) IF(userDataPointExists) THEN DofIdx=fieldVariable%COMPONENTS(componentNumber)%PARAM_TO_DOF_MAP% & @@ -23683,7 +23683,7 @@ SUBROUTINE Field_ParameterSetUpdateDataPointL(field,variableType,fieldSetType,us IF(ASSOCIATED(decomposition)) THEN decompositionTopology=>decomposition%TOPOLOGY IF(ASSOCIATED(decompositionTopology)) THEN - CALL DecompositionTopology_DataPointCheckExists(decompositionTopology,userDataPointNumber, & + CALL DecompositionTopologyDataPointCheckExists(decompositionTopology,userDataPointNumber, & & userDataPointExists,decompositionLocalDataPointNumber,ghostDataPoint,err,error,*999) IF(userDataPointExists) THEN DofIdx=fieldVariable%COMPONENTS(componentNumber)%PARAM_TO_DOF_MAP% & @@ -30544,7 +30544,7 @@ SUBROUTINE MESH_EMBEDDING_PULL_GAUSS_POINT_DATA(MESH_EMBEDDING,PARENT_FIELD,PARE INTEGER(INTG), INTENT(IN) :: CHILD_COMPONENT ! 3 MATRICES diff --git a/src/generated_mesh_routines.f90 b/src/generated_mesh_routines.f90 index e5af0eb3..7bae8187 100644 --- a/src/generated_mesh_routines.f90 +++ b/src/generated_mesh_routines.f90 @@ -1415,7 +1415,7 @@ SUBROUTINE GENERATED_MESH_REGULAR_CREATE_FINISH(GENERATED_MESH,MESH_USER_NUMBER, TYPE(COORDINATE_SYSTEM_TYPE), POINTER :: COORDINATE_SYSTEM TYPE(GENERATED_MESH_REGULAR_TYPE), POINTER :: REGULAR_MESH TYPE(INTERFACE_TYPE), POINTER :: INTERFACE - TYPE(MeshComponentElementsType), POINTER :: MESH_ELEMENTS + TYPE(MeshElementsType), POINTER :: MESH_ELEMENTS TYPE(NODES_TYPE), POINTER :: NODES TYPE(REGION_TYPE), POINTER :: REGION TYPE(VARYING_STRING) :: LOCAL_ERROR @@ -2244,7 +2244,7 @@ SUBROUTINE GENERATED_MESH_ELLIPSOID_CREATE_FINISH(GENERATED_MESH,MESH_USER_NUMBE INTEGER(INTG), ALLOCATABLE :: APEX_ELEMENT_NODES_USER_NUMBERS(:), WALL_ELEMENT_NODES_USER_NUMBERS(:) INTEGER(INTG), ALLOCATABLE :: NIDX(:,:,:),CORNER_NODES(:,:,:),EIDX(:,:,:) REAL(DP) :: DELTA(3),DELTAi(3) - TYPE(MeshComponentElementsType), POINTER :: MESH_ELEMENTS + TYPE(MeshElementsType), POINTER :: MESH_ELEMENTS CALL ENTERS("GENERATED_MESH_ELLIPSOID_CREATE_FINISH",ERR,ERROR,*999) @@ -2498,7 +2498,7 @@ SUBROUTINE GENERATED_MESH_CYLINDER_CREATE_FINISH(GENERATED_MESH,MESH_USER_NUMBER INTEGER(INTG), ALLOCATABLE :: ELEMENT_NODES(:),ELEMENT_NODES_USER_NUMBERS(:) INTEGER(INTG), ALLOCATABLE :: NIDX(:,:,:),EIDX(:,:,:) REAL(DP) :: DELTA(3),DELTAi(3) - TYPE(MeshComponentElementsType), POINTER :: MESH_ELEMENTS + TYPE(MeshElementsType), POINTER :: MESH_ELEMENTS CALL ENTERS("GENERATED_MESH_CYLINDER_CREATE_FINISH",ERR,ERROR,*999) diff --git a/src/interface_routines.f90 b/src/interface_routines.f90 index e6986a9a..e5efecd5 100644 --- a/src/interface_routines.f90 +++ b/src/interface_routines.f90 @@ -1732,7 +1732,7 @@ SUBROUTINE InterfacePointsConnectivity_ElementNumberSet(pointsConnectivity,dataP INTEGER(INTG), INTENT(OUT) :: err !0)) THEN IF (ALLOCATED(pointsConnectivity%pointsConnectivity)) THEN - CALL MESH_TOPOLOGY_ELEMENT_CHECK_EXISTS(pointsConnectivity%INTERFACE%COUPLED_MESHES(coupledMeshIndexNumber)%PTR, & - & meshComponentNumber,coupledMeshUserElementNumber,elementExists,elementGlobalNumber,err,error,*999) !Make sure user element exists + CALL MeshTopologyElementCheckExists(pointsConnectivity%INTERFACE%COUPLED_MESHES(coupledMeshIndexNumber)%PTR, & + & meshComponentNumber,coupledMeshUserElementNumber,elementExists,elementMeshNumber,err,error,*999) !Make sure user element exists IF(elementExists) THEN pointsConnectivity%pointsConnectivity(dataPointGlobalNumber,coupledMeshIndexNumber)%coupledMeshElementNumber= & - & elementGlobalNumber + & elementMeshNumber ELSE CALL FLAG_ERROR("Element with user number ("//TRIM(NUMBER_TO_VSTRING & & (coupledMeshUserElementNumber,"*",err,error))//") does not exist.",err,error,*999) diff --git a/src/mesh_routines.f90 b/src/mesh_routines.f90 index 0071dddc..229d0013 100644 --- a/src/mesh_routines.f90 +++ b/src/mesh_routines.f90 @@ -99,6 +99,16 @@ MODULE MESH_ROUTINES MODULE PROCEDURE MESH_USER_NUMBER_FIND_REGION END INTERFACE !MESH_USER_NUMBER_FIND + INTERFACE MeshTopologyNodeCheckExists + MODULE PROCEDURE MeshTopologyNodeCheckExistsMesh + MODULE PROCEDURE MeshTopologyNodeCheckExistsMeshNodes + END INTERFACE MeshTopologyNodeCheckExists + + INTERFACE MeshTopologyElementCheckExists + MODULE PROCEDURE MeshTopologyElementCheckExistsMesh + MODULE PROCEDURE MeshTopologyElementCheckExistsMeshElements + END INTERFACE MeshTopologyElementCheckExists + PUBLIC DECOMPOSITION_ALL_TYPE,DECOMPOSITION_CALCULATED_TYPE,DECOMPOSITION_USER_DEFINED_TYPE PUBLIC DECOMPOSITIONS_INITIALISE,DECOMPOSITIONS_FINALISE @@ -115,15 +125,15 @@ MODULE MESH_ROUTINES PUBLIC DECOMPOSITION_NUMBER_OF_DOMAINS_GET,DECOMPOSITION_NUMBER_OF_DOMAINS_SET - PUBLIC DECOMPOSITION_TOPOLOGY_ELEMENT_CHECK_EXISTS,DecompositionTopology_DataPointCheckExists + PUBLIC DECOMPOSITION_TOPOLOGY_ELEMENT_CHECK_EXISTS,DecompositionTopologyDataPointCheckExists - PUBLIC DecompositionTopology_DataProjectionCalculate + PUBLIC DecompositionTopologyDataProjectionCalculate - PUBLIC DecompositionTopology_ElementDataPointLocalNumberGet + PUBLIC DecompositionTopologyElementDataPointLocalNumberGet - PUBLIC DecompositionTopology_ElementDataPointUserNumberGet + PUBLIC DecompositionTopologyElementDataPointUserNumberGet - PUBLIC DecompositionTopology_NumberOfElementDataPointsGet + PUBLIC DecompositionTopologyNumberOfElementDataPointsGet PUBLIC DECOMPOSITION_TYPE_GET,DECOMPOSITION_TYPE_SET @@ -135,7 +145,7 @@ MODULE MESH_ROUTINES PUBLIC DOMAIN_TOPOLOGY_NODE_CHECK_EXISTS - PUBLIC MESH_TOPOLOGY_ELEMENT_CHECK_EXISTS,MESH_TOPOLOGY_NODE_CHECK_EXISTS + PUBLIC MeshTopologyElementCheckExists,MeshTopologyNodeCheckExists PUBLIC MESH_CREATE_START,MESH_CREATE_FINISH @@ -161,9 +171,19 @@ MODULE MESH_ROUTINES PUBLIC MESH_TOPOLOGY_ELEMENTS_ELEMENT_USER_NUMBER_GET,MESH_TOPOLOGY_ELEMENTS_ELEMENT_USER_NUMBER_SET - PUBLIC Mesh_TopologyElementsUserNumbersAllSet + PUBLIC MeshTopologyElementsUserNumbersAllSet + + PUBLIC MeshTopologyDataPointsCalculateProjection + + PUBLIC MeshTopologyNodeDerivativesGet - PUBLIC Mesh_TopologyDataPointsCalculateProjection + PUBLIC MeshTopologyNodeNumberOfDerivativesGet + + PUBLIC MeshTopologyNodeNumberOfVersionsGet + + PUBLIC MeshTopologyNodesDestroy + + PUBLIC MeshTopologyNodesGet PUBLIC MESH_USER_NUMBER_FIND, MESH_USER_NUMBER_TO_MESH @@ -791,7 +811,7 @@ SUBROUTINE DECOMPOSITION_ELEMENT_DOMAIN_GET(DECOMPOSITION,USER_ELEMENT_NUMBER,DO TYPE(VARYING_STRING) :: LOCAL_ERROR INTEGER(INTG) :: GLOBAL_ELEMENT_NUMBER TYPE(TREE_NODE_TYPE), POINTER :: TREE_NODE - TYPE(MeshComponentElementsType), POINTER :: MESH_ELEMENTS + TYPE(MeshElementsType), POINTER :: MESH_ELEMENTS CALL ENTERS("DECOMPOSITION_ELEMENT_DOMAIN_GET",ERR,ERROR,*999) @@ -1262,14 +1282,14 @@ END SUBROUTINE DecompositionTopology_DataPointsCalculate ! !>Calculates the decomposition element topology for a data projection (for data projections on fields). - SUBROUTINE DecompositionTopology_DataProjectionCalculate(decompositionTopology,err,error,*) + SUBROUTINE DecompositionTopologyDataProjectionCalculate(decompositionTopology,err,error,*) !Argument variables TYPE(DECOMPOSITION_TOPOLOGY_TYPE), POINTER :: decompositionTopology !Gets the local data point number for data points projected on an element - SUBROUTINE DecompositionTopology_ElementDataPointLocalNumberGet(decompositionTopology,elementNumber,dataPointIndex, & + SUBROUTINE DecompositionTopologyElementDataPointLocalNumberGet(decompositionTopology,elementNumber,dataPointIndex, & & dataPointLocalNumber,err,error,*) !Argument variables @@ -1305,7 +1325,7 @@ SUBROUTINE DecompositionTopology_ElementDataPointLocalNumberGet(decompositionTop INTEGER(INTG) :: numberOfDataPoints TYPE(VARYING_STRING) :: localError - CALL ENTERS("DecompositionTopology_ElementDataPointLocalNumberGet",err,error,*999) + CALL ENTERS("DecompositionTopologyElementDataPointLocalNumberGet",err,error,*999) IF(ASSOCIATED(decompositionTopology)) THEN decompositionData=>decompositionTopology%dataPoints @@ -1325,19 +1345,19 @@ SUBROUTINE DecompositionTopology_ElementDataPointLocalNumberGet(decompositionTop CALL FLAG_ERROR("Decomposition topology is not associated.",err,error,*999) ENDIF - CALL EXITS("DecompositionTopology_ElementDataPointLocalNumberGet") + CALL EXITS("DecompositionTopologyElementDataPointLocalNumberGet") RETURN -999 CALL ERRORS("DecompositionTopology_ElementDataPointLocalNumberGet",err,error) - CALL EXITS("DecompositionTopology_ElementDataPointLocalNumberGet") +999 CALL ERRORS("DecompositionTopologyElementDataPointLocalNumberGet",err,error) + CALL EXITS("DecompositionTopologyElementDataPointLocalNumberGet") RETURN 1 - END SUBROUTINE DecompositionTopology_ElementDataPointLocalNumberGet + END SUBROUTINE DecompositionTopologyElementDataPointLocalNumberGet ! !================================================================================================================================ ! !>Gets the user (global) data point number for data points projected on an element - SUBROUTINE DecompositionTopology_ElementDataPointUserNumberGet(decompositionTopology,userElementNumber,dataPointIndex, & + SUBROUTINE DecompositionTopologyElementDataPointUserNumberGet(decompositionTopology,userElementNumber,dataPointIndex, & & dataPointUserNumber,err,error,*) !Argument variables @@ -1353,7 +1373,7 @@ SUBROUTINE DecompositionTopology_ElementDataPointUserNumberGet(decompositionTopo LOGICAL :: ghostElement,userElementExists TYPE(VARYING_STRING) :: localError - CALL ENTERS("DecompositionTopology_ElementDataPointUserNumberGet",err,error,*999) + CALL ENTERS("DecompositionTopologyElementDataPointUserNumberGet",err,error,*999) IF(ASSOCIATED(decompositionTopology)) THEN decompositionData=>decompositionTopology%dataPoints @@ -1389,19 +1409,19 @@ SUBROUTINE DecompositionTopology_ElementDataPointUserNumberGet(decompositionTopo CALL FLAG_ERROR("Decomposition topology is not associated.",err,error,*999) ENDIF - CALL EXITS("DecompositionTopology_ElementDataPointUserNumberGet") + CALL EXITS("DecompositionTopologyElementDataPointUserNumberGet") RETURN -999 CALL ERRORS("DecompositionTopology_ElementDataPointUserNumberGet",err,error) - CALL EXITS("DecompositionTopology_ElementDataPointUserNumberGet") +999 CALL ERRORS("DecompositionTopologyElementDataPointUserNumberGet",err,error) + CALL EXITS("DecompositionTopologyElementDataPointUserNumberGet") RETURN 1 - END SUBROUTINE DecompositionTopology_ElementDataPointUserNumberGet + END SUBROUTINE DecompositionTopologyElementDataPointUserNumberGet ! !================================================================================================================================ ! !>Gets the number of data points projected on an element - SUBROUTINE DecompositionTopology_NumberOfElementDataPointsGet(decompositionTopology,userElementNumber, & + SUBROUTINE DecompositionTopologyNumberOfElementDataPointsGet(decompositionTopology,userElementNumber, & & numberOfDataPoints,err,error,*) !Argument variables @@ -1416,7 +1436,7 @@ SUBROUTINE DecompositionTopology_NumberOfElementDataPointsGet(decompositionTopol LOGICAL :: ghostElement,userElementExists TYPE(VARYING_STRING) :: localError - CALL ENTERS("DecompositionTopology_NumberOfElementDataPointsGet",err,error,*999) + CALL ENTERS("DecompositionTopologyNumberOfElementDataPointsGet",err,error,*999) IF(ASSOCIATED(decompositionTopology)) THEN decompositionData=>decompositionTopology%dataPoints @@ -1444,19 +1464,19 @@ SUBROUTINE DecompositionTopology_NumberOfElementDataPointsGet(decompositionTopol CALL FLAG_ERROR("Decomposition topology is not associated.",err,error,*999) ENDIF - CALL EXITS("DecompositionTopology_NumberOfElementDataPointsGet") + CALL EXITS("DecompositionTopologyNumberOfElementDataPointsGet") RETURN -999 CALL ERRORS("DecompositionTopology_NumberOfElementDataPointsGet",err,error) - CALL EXITS("DecompositionTopology_NumberOfElementDataPointsGet") +999 CALL ERRORS("DecompositionTopologyNumberOfElementDataPointsGet",err,error) + CALL EXITS("DecompositionTopologyNumberOfElementDataPointsGet") RETURN 1 - END SUBROUTINE DecompositionTopology_NumberOfElementDataPointsGet + END SUBROUTINE DecompositionTopologyNumberOfElementDataPointsGet ! !================================================================================================================================ ! !>Checks that a user element number exists in a decomposition. - SUBROUTINE DecompositionTopology_DataPointCheckExists(decompositionTopology,userDataPointNumber,userDataPointExists, & + SUBROUTINE DecompositionTopologyDataPointCheckExists(decompositionTopology,userDataPointNumber,userDataPointExists, & & decompositionLocalDataPointNumber,ghostDataPoint,err,error,*) !Argument variables @@ -1471,7 +1491,7 @@ SUBROUTINE DecompositionTopology_DataPointCheckExists(decompositionTopology,user TYPE(DecompositionDataPointsType), POINTER :: decompositionData TYPE(TREE_NODE_TYPE), POINTER :: treeNode - CALL ENTERS("DecompositionTopology_DataPointCheckExists",ERR,error,*999) + CALL ENTERS("DecompositionTopologyDataPointCheckExists",ERR,error,*999) userDataPointExists=.FALSE. decompositionLocalDataPointNumber=0 @@ -1493,13 +1513,13 @@ SUBROUTINE DecompositionTopology_DataPointCheckExists(decompositionTopology,user CALL FLAG_ERROR("Decomposition topology is not associated.",err,error,*999) ENDIF - CALL EXITS("DecompositionTopology_DataPointCheckExists") + CALL EXITS("DecompositionTopologyDataPointCheckExists") RETURN -999 CALL ERRORS("DecompositionTopology_DataPointCheckExists",err,error) - CALL EXITS("DecompositionTopology_DataPointCheckExists") +999 CALL ERRORS("DecompositionTopologyDataPointCheckExists",err,error) + CALL EXITS("DecompositionTopologyDataPointCheckExists") RETURN 1 - END SUBROUTINE DecompositionTopology_DataPointCheckExists + END SUBROUTINE DecompositionTopologyDataPointCheckExists ! !================================================================================================================================ @@ -1978,7 +1998,7 @@ SUBROUTINE DECOMPOSITION_TOPOLOGY_ELEMENTS_CALCULATE(TOPOLOGY,ERR,ERROR,*) TYPE(DOMAIN_MAPPINGS_TYPE), POINTER :: DOMAIN_MAPPINGS TYPE(DOMAIN_TOPOLOGY_TYPE), POINTER :: DOMAIN_TOPOLOGY TYPE(MESH_TYPE), POINTER :: MESH - TYPE(MeshComponentElementsType), POINTER :: MESH_ELEMENTS + TYPE(MeshElementsType), POINTER :: MESH_ELEMENTS TYPE(MeshComponentTopologyType), POINTER :: MESH_TOPOLOGY CALL ENTERS("DECOMPOSITION_TOPOLOGY_ELEMENTS_CALCULATE",ERR,ERROR,*999) @@ -4015,8 +4035,8 @@ SUBROUTINE DOMAIN_MAPPINGS_ELEMENTS_CALCULATE(DOMAIN,ERR,ERROR,*) CALL LIST_ITEM_ADD(ADJACENT_DOMAINS_LIST,domain_no,ERR,ERROR,*999) DO nn=1,BASIS%NUMBER_OF_NODES np=MESH%TOPOLOGY(component_idx)%PTR%ELEMENTS%ELEMENTS(ne)%MESH_ELEMENT_NODES(nn) - DO no_adjacent_element=1,MESH%TOPOLOGY(component_idx)%PTR%NODES%NODES(np)%NUMBER_OF_SURROUNDING_ELEMENTS - adjacent_element=MESH%TOPOLOGY(component_idx)%PTR%NODES%NODES(np)%SURROUNDING_ELEMENTS(no_adjacent_element) + DO no_adjacent_element=1,MESH%TOPOLOGY(component_idx)%PTR%NODES%NODES(np)%numberOfSurroundingElements + adjacent_element=MESH%TOPOLOGY(component_idx)%PTR%NODES%NODES(np)%surroundingElements(no_adjacent_element) IF(DECOMPOSITION%ELEMENT_DOMAIN(adjacent_element)/=domain_no) THEN CALL LIST_ITEM_ADD(ADJACENT_ELEMENTS_LIST(domain_no)%PTR,adjacent_element,ERR,ERROR,*999) CALL LIST_ITEM_ADD(ADJACENT_DOMAINS_LIST,DECOMPOSITION%ELEMENT_DOMAIN(adjacent_element),ERR,ERROR,*999) @@ -4360,15 +4380,15 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) IF(ERR/=0) GOTO 999 !Calculate the local and global numbers and set up the mappings - ALLOCATE(NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(MESH_TOPOLOGY%NODES%NUMBER_OF_NODES),STAT=ERR) + ALLOCATE(NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(MESH_TOPOLOGY%NODES%numberOfNodes),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node mapping global to local map.",ERR,ERROR,*999) - NODES_MAPPING%NUMBER_OF_GLOBAL=MESH_TOPOLOGY%NODES%NUMBER_OF_NODES + NODES_MAPPING%NUMBER_OF_GLOBAL=MESH_TOPOLOGY%NODES%numberOfNodes ALLOCATE(LOCAL_NODE_NUMBERS(0:DECOMPOSITION%NUMBER_OF_DOMAINS-1),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate local node numbers.",ERR,ERROR,*999) LOCAL_NODE_NUMBERS=0 - ALLOCATE(DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(MESH_TOPOLOGY%DOFS%NUMBER_OF_DOFS),STAT=ERR) + ALLOCATE(DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(MESH_TOPOLOGY%dofs%numberOfDofs),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate dofs mapping global to local map.",ERR,ERROR,*999) - DOFS_MAPPING%NUMBER_OF_GLOBAL=MESH_TOPOLOGY%DOFS%NUMBER_OF_DOFS + DOFS_MAPPING%NUMBER_OF_GLOBAL=MESH_TOPOLOGY%DOFS%numberOfDofs ALLOCATE(LOCAL_DOF_NUMBERS(0:DECOMPOSITION%NUMBER_OF_DOMAINS-1),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate local dof numbers.",ERR,ERROR,*999) LOCAL_DOF_NUMBERS=0 @@ -4378,7 +4398,7 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) NULLIFY(GHOST_NODES_LIST(domain_idx)%PTR) CALL LIST_CREATE_START(GHOST_NODES_LIST(domain_idx)%PTR,ERR,ERROR,*999) CALL LIST_DATA_TYPE_SET(GHOST_NODES_LIST(domain_idx)%PTR,LIST_INTG_TYPE,ERR,ERROR,*999) - CALL LIST_INITIAL_SIZE_SET(GHOST_NODES_LIST(domain_idx)%PTR,INT(MESH_TOPOLOGY%NODES%NUMBER_OF_NODES/2), & + CALL LIST_INITIAL_SIZE_SET(GHOST_NODES_LIST(domain_idx)%PTR,INT(MESH_TOPOLOGY%NODES%numberOfNodes/2), & & ERR,ERROR,*999) CALL LIST_CREATE_FINISH(GHOST_NODES_LIST(domain_idx)%PTR,ERR,ERROR,*999) ENDDO !domain_idx @@ -4390,7 +4410,7 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) NUMBER_BOUNDARY_NODES=0 !For the first pass just determine the internal and boundary nodes - DO node_idx=1,MESH_TOPOLOGY%NODES%NUMBER_OF_NODES + DO node_idx=1,MESH_TOPOLOGY%NODES%numberOfNodes NULLIFY(ADJACENT_DOMAINS_LIST) CALL LIST_CREATE_START(ADJACENT_DOMAINS_LIST,ERR,ERROR,*999) CALL LIST_DATA_TYPE_SET(ADJACENT_DOMAINS_LIST,LIST_INTG_TYPE,ERR,ERROR,*999) @@ -4401,8 +4421,8 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) CALL LIST_DATA_TYPE_SET(ALL_ADJACENT_DOMAINS_LIST,LIST_INTG_TYPE,ERR,ERROR,*999) CALL LIST_INITIAL_SIZE_SET(ALL_ADJACENT_DOMAINS_LIST,DECOMPOSITION%NUMBER_OF_DOMAINS,ERR,ERROR,*999) CALL LIST_CREATE_FINISH(ALL_ADJACENT_DOMAINS_LIST,ERR,ERROR,*999) - DO no_adjacent_element=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%NUMBER_OF_SURROUNDING_ELEMENTS - adjacent_element=MESH_TOPOLOGY%NODES%NODES(node_idx)%SURROUNDING_ELEMENTS(no_adjacent_element) + DO no_adjacent_element=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%numberOfSurroundingElements + adjacent_element=MESH_TOPOLOGY%NODES%NODES(node_idx)%surroundingElements(no_adjacent_element) domain_no=DECOMPOSITION%ELEMENT_DOMAIN(adjacent_element) CALL LIST_ITEM_ADD(ADJACENT_DOMAINS_LIST,domain_no,ERR,ERROR,*999) DO domain_idx=1,ELEMENTS_MAPPING%GLOBAL_TO_LOCAL_MAP(adjacent_element)%NUMBER_OF_DOMAINS @@ -4433,9 +4453,9 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node global to local map domain number.",ERR,ERROR,*999) ALLOCATE(NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%LOCAL_TYPE(MAX_NUMBER_DOMAINS),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node global to local map local type.",ERR,ERROR,*999) - DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS - ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%DOF_INDEX(version_idx) + DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%numberOfDerivatives + DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%numberOfVersions + ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%dofIndex(version_idx) ALLOCATE(DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%LOCAL_NUMBER(MAX_NUMBER_DOMAINS),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate dof global to local map local number.",ERR,ERROR,*999) ALLOCATE(DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%DOMAIN_NUMBER(MAX_NUMBER_DOMAINS),STAT=ERR) @@ -4454,9 +4474,9 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%LOCAL_NUMBER(1)=-1 NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%DOMAIN_NUMBER(1)=DOMAINS(1) NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%LOCAL_TYPE(1)=DOMAIN_LOCAL_INTERNAL - DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS - ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%DOF_INDEX(version_idx) + DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%numberOfDerivatives + DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%numberOfVersions + ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%dofIndex(version_idx) DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%NUMBER_OF_DOMAINS=1 DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%LOCAL_NUMBER(1)=-1 DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%DOMAIN_NUMBER(1)=domain_no @@ -4466,9 +4486,9 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) ELSE !Node is on the boundary of computational domains NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%NUMBER_OF_DOMAINS=NUMBER_OF_DOMAINS - DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS - ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%DOF_INDEX(version_idx) + DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%numberOfDerivatives + DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%numberOfVersions + ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%dofIndex(version_idx) DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%NUMBER_OF_DOMAINS=NUMBER_OF_DOMAINS ENDDO !version_idx ENDDO !derivative_idx @@ -4479,9 +4499,9 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%LOCAL_NUMBER(domain_idx)=-1 NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%DOMAIN_NUMBER(domain_idx)=domain_no NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%LOCAL_TYPE(domain_idx)=DOMAIN_LOCAL_BOUNDARY - DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS - ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%DOF_INDEX(version_idx) + DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%numberOfDerivatives + DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%numberOfVersions + ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%dofIndex(version_idx) DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%LOCAL_NUMBER(domain_idx)=-1 DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%DOMAIN_NUMBER(domain_idx)=domain_no DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%LOCAL_TYPE(domain_idx)=DOMAIN_LOCAL_BOUNDARY @@ -4494,20 +4514,20 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) ENDDO !node_idx !For the second pass assign boundary nodes to one domain on the boundary and set local node numbers. - NUMBER_OF_NODES_PER_DOMAIN=FLOOR(REAL(MESH_TOPOLOGY%NODES%NUMBER_OF_NODES,DP)/ & + NUMBER_OF_NODES_PER_DOMAIN=FLOOR(REAL(MESH_TOPOLOGY%NODES%numberOfNodes,DP)/ & & REAL(DECOMPOSITION%NUMBER_OF_DOMAINS,DP)) - ALLOCATE(DOMAIN%NODE_DOMAIN(MESH_TOPOLOGY%NODES%NUMBER_OF_NODES),STAT=ERR) + ALLOCATE(DOMAIN%NODE_DOMAIN(MESH_TOPOLOGY%NODES%numberOfNodes),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node domain",ERR,ERROR,*999) DOMAIN%NODE_DOMAIN=-1 - DO node_idx=1,MESH_TOPOLOGY%NODES%NUMBER_OF_NODES + DO node_idx=1,MESH_TOPOLOGY%NODES%numberOfNodes IF(NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%NUMBER_OF_DOMAINS==1) THEN !Internal node domain_no=NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%DOMAIN_NUMBER(1) DOMAIN%NODE_DOMAIN(node_idx)=domain_no LOCAL_NODE_NUMBERS(domain_no)=LOCAL_NODE_NUMBERS(domain_no)+1 NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%LOCAL_NUMBER(1)=LOCAL_NODE_NUMBERS(domain_no) - DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS - ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%DOF_INDEX(version_idx) + DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%numberOfDerivatives + DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%numberOfVersions + ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%dofIndex(version_idx) LOCAL_DOF_NUMBERS(domain_no)=LOCAL_DOF_NUMBERS(domain_no)+1 DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%LOCAL_NUMBER(1)=LOCAL_DOF_NUMBERS(domain_no) ENDDO !version_idx @@ -4529,9 +4549,9 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%LOCAL_NUMBER(1)=LOCAL_NODE_NUMBERS(domain_no) NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%DOMAIN_NUMBER(1)=domain_no NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(node_idx)%LOCAL_TYPE(1)=DOMAIN_LOCAL_BOUNDARY - DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS - ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%DOF_INDEX(version_idx) + DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%numberOfDerivatives + DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%numberOfVersions + ny=MESH_TOPOLOGY%NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%dofIndex(version_idx) LOCAL_DOF_NUMBERS(domain_no)=LOCAL_DOF_NUMBERS(domain_no)+1 DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%NUMBER_OF_DOMAINS=1 DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%LOCAL_NUMBER(1)=LOCAL_DOF_NUMBERS(domain_no) @@ -4569,9 +4589,9 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(ghost_node)%LOCAL_TYPE( & & NODES_MAPPING%GLOBAL_TO_LOCAL_MAP(ghost_node)%NUMBER_OF_DOMAINS)= & & DOMAIN_LOCAL_GHOST - DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(ghost_node)%NUMBER_OF_DERIVATIVES - DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(ghost_node)%DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS - ny=MESH_TOPOLOGY%NODES%NODES(ghost_node)%DERIVATIVES(derivative_idx)%DOF_INDEX(version_idx) + DO derivative_idx=1,MESH_TOPOLOGY%NODES%NODES(ghost_node)%numberOfDerivatives + DO version_idx=1,MESH_TOPOLOGY%NODES%NODES(ghost_node)%DERIVATIVES(derivative_idx)%numberOfVersions + ny=MESH_TOPOLOGY%NODES%NODES(ghost_node)%DERIVATIVES(derivative_idx)%dofIndex(version_idx) LOCAL_DOF_NUMBERS(domain_idx)=LOCAL_DOF_NUMBERS(domain_idx)+1 DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%NUMBER_OF_DOMAINS= & & DOFS_MAPPING%GLOBAL_TO_LOCAL_MAP(ny)%NUMBER_OF_DOMAINS+1 @@ -4593,7 +4613,7 @@ SUBROUTINE DOMAIN_MAPPINGS_NODES_DOFS_CALCULATE(DOMAIN,ERR,ERROR,*) ALLOCATE(NODE_COUNT(0:number_computational_nodes-1),STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node count.",ERR,ERROR,*999) NODE_COUNT=0 - DO node_idx=1,MESH_TOPOLOGY%NODES%NUMBER_OF_NODES + DO node_idx=1,MESH_TOPOLOGY%NODES%numberOfNodes no_computational_node=DOMAIN%NODE_DOMAIN(node_idx) IF(no_computational_node>=0.AND.no_computational_node> DEALLOCATE(MESH) @@ -6509,6 +6529,66 @@ END SUBROUTINE MESH_FINALISE !================================================================================================================================ ! + !>Returns a region nodes pointer corresponding to the mesh global nodes accounting for interfaces. + SUBROUTINE MeshGlobalNodesGet(mesh,nodes,err,error,*) + + !Argument variables + TYPE(MESH_TYPE), POINTER :: mesh !mesh%region + IF(ASSOCIATED(region)) THEN + nodes=>region%nodes + ELSE + INTERFACE=>mesh%INTERFACE + IF(ASSOCIATED(interface)) THEN + nodes=>interface%nodes + ELSE + localError="Mesh number "//TRIM(NumberToVString(mesh%USER_NUMBER,"*",err,error))// & + & " does not have an associated region or interface." + CALL FlagError(localError,err,error,*999) + ENDIF + ENDIF + IF(.NOT.ASSOCIATED(nodes)) THEN + IF(ASSOCIATED(region)) THEN + localError="Mesh number "//TRIM(NumberToVString(mesh%USER_NUMBER,"*",err,error))// & + & " does not have any nodes associated with the mesh region." + ELSE + localError="Mesh number "//TRIM(NumberToVString(mesh%USER_NUMBER,"*",err,error))// & + & " does not have any nodes associated with the mesh interface." + ENDIF + CALL FlagError(localError,err,error,*999) + ENDIF + ENDIF + ELSE + CALL FlagError("Mesh is not associated.",err,error,*999) + ENDIF + + CALL Exits("MeshGlobalNodesGet") + RETURN +999 CALL Errors("MeshGlobalNodesGet",err,error) + CALL Exits("MeshGlobalNodesGet") + RETURN 1 + + END SUBROUTINE MeshGlobalNodesGet + + ! + !================================================================================================================================ + ! + !>Initialises a mesh. SUBROUTINE MESH_INITIALISE(MESH,ERR,ERROR,*) @@ -6624,16 +6704,16 @@ SUBROUTINE MESH_NUMBER_OF_COMPONENTS_SET(MESH,NUMBER_OF_COMPONENTS,ERR,ERROR,*) DO component_idx=MESH%NUMBER_OF_COMPONENTS+1,NUMBER_OF_COMPONENTS ALLOCATE(NEW_TOPOLOGY(component_idx)%PTR,STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate new topology component",ERR,ERROR,*999) - NEW_TOPOLOGY(component_idx)%PTR%MESH=>MESH - NEW_TOPOLOGY(component_idx)%PTR%MESH_COMPONENT_NUMBER=component_idx - NULLIFY(NEW_TOPOLOGY(component_idx)%PTR%ELEMENTS) - NULLIFY(NEW_TOPOLOGY(component_idx)%PTR%NODES) - NULLIFY(NEW_TOPOLOGY(component_idx)%PTR%DOFS) + NEW_TOPOLOGY(component_idx)%PTR%mesh=>mesh + NEW_TOPOLOGY(component_idx)%PTR%meshComponentNumber=component_idx + NULLIFY(NEW_TOPOLOGY(component_idx)%PTR%elements) + NULLIFY(NEW_TOPOLOGY(component_idx)%PTR%nodes) + NULLIFY(NEW_TOPOLOGY(component_idx)%PTR%dofs) NULLIFY(NEW_TOPOLOGY(component_idx)%PTR%dataPoints) !Initialise the topology components CALL MESH_TOPOLOGY_ELEMENTS_INITIALISE(NEW_TOPOLOGY(component_idx)%PTR,ERR,ERROR,*999) - CALL MESH_TOPOLOGY_NODES_INITIALISE(NEW_TOPOLOGY(component_idx)%PTR,ERR,ERROR,*999) - CALL MESH_TOPOLOGY_DOFS_INITIALISE(NEW_TOPOLOGY(component_idx)%PTR,ERR,ERROR,*999) + CALL MeshTopologyNodesInitialise(NEW_TOPOLOGY(component_idx)%PTR,ERR,ERROR,*999) + CALL MeshTopologyDofsInitialise(NEW_TOPOLOGY(component_idx)%PTR,ERR,ERROR,*999) CALL MESH_TOPOLOGY_DATA_POINTS_INITIALISE(NEW_TOPOLOGY(component_idx)%PTR,ERR,ERROR,*999) ENDDO !component_idx ENDIF @@ -6757,6 +6837,63 @@ END SUBROUTINE MESH_NUMBER_OF_ELEMENTS_SET !================================================================================================================================ ! + !>Returns the region for a mesh accounting for regions and interfaces + SUBROUTINE MeshRegionGet(mesh,region,err,error,*) + + !Argument variables + TYPE(MESH_TYPE), POINTER :: mesh !mesh%region + IF(.NOT.ASSOCIATED(region)) THEN + interface=>mesh%interface + IF(ASSOCIATED(interface)) THEN + parentRegion=>interface%PARENT_REGION + IF(ASSOCIATED(parentRegion)) THEN + region=>parentRegion + ELSE + localError="The parent region not associated for mesh number "// & + & TRIM(NumberToVString(mesh%USER_NUMBER,"*",err,error))//" of interface number "// & + & TRIM(NumberToVString(interface%USER_NUMBER,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + ELSE + localError="The region or interface is not associated for mesh number "// & + & TRIM(NumberToVString(mesh%USER_NUMBER,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + ENDIF + ENDIF + ELSE + CALL FlagError("Mesh is not associated.",err,error,*999) + ENDIF + + CALL Exits("MeshRegionGet") + RETURN +999 CALL Errors("MeshRegionGet",err,error) + CALL Exits("MeshRegionGet") + RETURN 1 + + END SUBROUTINE MeshRegionGet + + ! + !================================================================================================================================ + ! + !>Changes/sets the surrounding elements calculate flag. \see OPENCMISS::CMISSMeshSurroundingElementsCalculateSet SUBROUTINE MESH_SURROUNDING_ELEMENTS_CALCULATE_SET(MESH,SURROUNDING_ELEMENTS_CALCULATE_FLAG,ERR,ERROR,*) @@ -6791,273 +6928,273 @@ END SUBROUTINE MESH_SURROUNDING_ELEMENTS_CALCULATE_SET ! !>Calculates the mesh topology. - SUBROUTINE MESH_TOPOLOGY_CALCULATE(TOPOLOGY,SURROUNDING_ELEMENTS_CALCULATE,ERR,ERROR,*) + SUBROUTINE MeshTopologyCalculate(topology,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !Calculates the boundary nodes and elements for a mesh topology. - SUBROUTINE MESH_TOPOLOGY_BOUNDARY_CALCULATE(TOPOLOGY,ERR,ERROR,*) + SUBROUTINE MeshTopologyBoundaryCalculate(topology,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%BASIS - SELECT CASE(BASIS%TYPE) + IF(ASSOCIATED(topology)) THEN + nodes=>topology%nodes + IF(ASSOCIATED(nodes)) THEN + elements=>topology%elements + IF(ASSOCIATED(elements)) THEN + DO elementIdx=1,elements%NUMBER_OF_ELEMENTS + basis=>elements%elements(elementIdx)%basis + SELECT CASE(basis%type) CASE(BASIS_LAGRANGE_HERMITE_TP_TYPE) - DO nic=-BASIS%NUMBER_OF_XI_COORDINATES,BASIS%NUMBER_OF_XI_COORDINATES - IF(nic/=0) THEN - IF(TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%ADJACENT_ELEMENTS(nic)%NUMBER_OF_ADJACENT_ELEMENTS==0) THEN - TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%BOUNDARY_ELEMENT=.TRUE. - IF(nic<0) THEN - XI_DIRECTION=-nic - MATCH_INDEX=1 + DO xiCoordIdx=-basis%NUMBER_OF_XI_COORDINATES,basis%NUMBER_OF_XI_COORDINATES + IF(xiCoordIdx/=0) THEN + IF(elements%elements(elementIdx)%ADJACENT_ELEMENTS(xiCoordIdx)%NUMBER_OF_ADJACENT_ELEMENTS==0) THEN + elements%elements(elementIdx)%BOUNDARY_ELEMENT=.TRUE. + IF(xiCoordIdx<0) THEN + xiDirection=-xiCoordIdx + matchIndex=1 ELSE - XI_DIRECTION=nic - MATCH_INDEX=BASIS%NUMBER_OF_NODES_XIC(nic) + xiDirection=xiCoordIdx + matchIndex=BASIS%NUMBER_OF_NODES_XIC(xiCoordIdx) ENDIF - DO local_node_idx=1,BASIS%NUMBER_OF_NODES - IF(BASIS%NODE_POSITION_INDEX(local_node_idx,XI_DIRECTION)==MATCH_INDEX) THEN - node_idx=TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%MESH_ELEMENT_NODES(local_node_idx) - TOPOLOGY%NODES%NODES(node_idx)%BOUNDARY_NODE=.TRUE. + DO localNodeIdx=1,BASIS%NUMBER_OF_NODES + IF(basis%NODE_POSITION_INDEX(localNodeIdx,XIDIRECTION)==matchIndex) THEN + nodeIdx=elements%elements(elementIdx)%MESH_ELEMENT_NODES(localNodeIdx) + nodes%nodes(nodeIdx)%boundaryNode=.TRUE. ENDIF ENDDO !nn ENDIF ENDIF - ENDDO !nic + ENDDO !xiCoordIdx CASE(BASIS_SIMPLEX_TYPE) - TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%BOUNDARY_ELEMENT=.FALSE. - DO nic=1,BASIS%NUMBER_OF_XI_COORDINATES - TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%BOUNDARY_ELEMENT=TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)% & - & BOUNDARY_ELEMENT.OR.TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%ADJACENT_ELEMENTS(nic)% & - & NUMBER_OF_ADJACENT_ELEMENTS==0 - IF(TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%ADJACENT_ELEMENTS(nic)%NUMBER_OF_ADJACENT_ELEMENTS==0) THEN - DO local_node_idx=1,BASIS%NUMBER_OF_NODES - IF(BASIS%NODE_POSITION_INDEX(local_node_idx,nic)==1) THEN - node_idx=TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%MESH_ELEMENT_NODES(local_node_idx) - TOPOLOGY%NODES%NODES(node_idx)%BOUNDARY_NODE=.TRUE. + elements%elements(elementIdx)%BOUNDARY_ELEMENT=.FALSE. + DO xiCoordIdx=1,basis%NUMBER_OF_XI_COORDINATES + elements%elements(elementIdx)%BOUNDARY_ELEMENT=elements%elements(elementIdx)%BOUNDARY_ELEMENT.OR. & + & elements%elements(elementIdx)%ADJACENT_ELEMENTS(xiCoordIdx)%NUMBER_OF_ADJACENT_ELEMENTS==0 + IF(elements%elements(elementIdx)%ADJACENT_ELEMENTS(xiCoordIdx)%NUMBER_OF_ADJACENT_ELEMENTS==0) THEN + DO localNodeIdx=1,basis%NUMBER_OF_NODES + IF(basis%NODE_POSITION_INDEX(localNodeIdx,xiCoordIdx)==1) THEN + nodeIdx=elements%elements(elementIdx)%MESH_ELEMENT_NODES(localNodeIdx) + nodes%nodes(nodeIdx)%boundaryNode=.TRUE. ENDIF - ENDDO !local_node_idx - + ENDDO !localNodeIdx ENDIF - ENDDO !nic + ENDDO !xiCoordIdx CASE(BASIS_SERENDIPITY_TYPE) - CALL FLAG_ERROR("Not implemented.",ERR,ERROR,*999) + CALL FlagError("Not implemented.",err,error,*999) CASE(BASIS_AUXILLIARY_TYPE) - CALL FLAG_ERROR("Not implemented.",ERR,ERROR,*999) + CALL FlagError("Not implemented.",err,error,*999) CASE(BASIS_B_SPLINE_TP_TYPE) - CALL FLAG_ERROR("Not implemented.",ERR,ERROR,*999) + CALL FlagError("Not implemented.",err,error,*999) CASE(BASIS_FOURIER_LAGRANGE_HERMITE_TP_TYPE) - CALL FLAG_ERROR("Not implemented.",ERR,ERROR,*999) + CALL FlagError("Not implemented.",err,error,*999) CASE(BASIS_EXTENDED_LAGRANGE_TP_TYPE) - CALL FLAG_ERROR("Not implemented.",ERR,ERROR,*999) + CALL FlagError("Not implemented.",err,error,*999) CASE DEFAULT - LOCAL_ERROR="The basis type of "//TRIM(NUMBER_TO_VSTRING(BASIS%TYPE,"*",ERR,ERROR))// & - & " is invalid." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + localError="The basis type of "//TRIM(NumberToVString(basis%TYPE,"*",err,error))//" is invalid." + CALL FlagError(localError,err,error,*999) END SELECT - ENDDO !element_idx + ENDDO !elementIdx ELSE - CALL FLAG_ERROR("Topology elements is not associated.",ERR,ERROR,*999) + CALL FlagError("Topology elements is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Topology nodes is not associated.",ERR,ERROR,*999) + CALL FlagError("Topology nodes is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Topology is not associated.",ERR,ERROR,*999) - ENDIF - - IF(DIAGNOSTICS1) THEN - CALL WRITE_STRING(DIAGNOSTIC_OUTPUT_TYPE,"Boundary elements:",ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,"Number of elements = ",TOPOLOGY%ELEMENTS%NUMBER_OF_ELEMENTS, & - ERR,ERROR,*999) - DO element_idx=1,TOPOLOGY%ELEMENTS%NUMBER_OF_ELEMENTS - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,"Element : ",element_idx,ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Boundary element = ",TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)% & - & BOUNDARY_ELEMENT,ERR,ERROR,*999) - ENDDO !element_idx - CALL WRITE_STRING(DIAGNOSTIC_OUTPUT_TYPE,"Boundary nodes:",ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,"Number of nodes = ",TOPOLOGY%NODES%NUMBER_OF_NODES, & - ERR,ERROR,*999) - DO node_idx=1,TOPOLOGY%NODES%NUMBER_OF_NODES - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,"Node : ",node_idx,ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Boundary node = ",TOPOLOGY%NODES%NODES(node_idx)% & - & BOUNDARY_NODE,ERR,ERROR,*999) - ENDDO !element_idx + CALL FlagError("Topology is not associated.",err,error,*999) + ENDIF + + IF(diagnostics1) THEN + CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"Boundary elements:",err,error,*999) + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"Number of elements = ",elements%NUMBER_OF_ELEMENTS,err,error,*999) + DO elementIdx=1,elements%NUMBER_OF_ELEMENTS + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"Element : ",elementIdx,err,error,*999) + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Boundary element = ",elements%elements(elementIdx)%BOUNDARY_ELEMENT, & + & err,error,*999) + ENDDO !elementIdx + CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"Boundary nodes:",err,error,*999) + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"Number of nodes = ",nodes%numberOfNodes,err,error,*999) + DO nodeIdx=1,nodes%numberOfNodes + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"Node : ",nodeIdx,err,error,*999) + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Boundary node = ",nodes%nodes(nodeIdx)%boundaryNode,err,error,*999) + ENDDO !elementIdx ENDIF - CALL EXITS("MESH_TOPOLOGY_BOUNDARY_CALCULATE") + CALL EXITS("MeshTopologyBoundaryCalculate") RETURN -999 CALL ERRORS("MESH_TOPOLOGY_BOUNDARY_CALCULATE",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_BOUNDARY_CALCULATE") +999 CALL ERRORS("MeshTopologyBoundaryCalculate",err,error) + CALL EXITS("MeshTopologyBoundaryCalculate") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_BOUNDARY_CALCULATE + END SUBROUTINE MeshTopologyBoundaryCalculate ! !=============================================================================================================================== ! !>Calculates the degrees-of-freedom for a mesh topology. - SUBROUTINE MESH_TOPOLOGY_DOFS_CALCULATE(TOPOLOGY,ERR,ERROR,*) + SUBROUTINE MeshTopologyDofsCalculate(topology,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !topology%nodes + IF(ASSOCIATED(nodes)) THEN + dofs=>topology%dofs + IF(ASSOCIATED(dofs)) THEN + numberOfDofs=0 + DO nodeIdx=1,nodes%numberOfNodes + DO derivativeIdx=1,nodes%nodes(nodeIdx)%numberOfDerivatives + ALLOCATE(nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%dofIndex( & + & nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%numberOfVersions),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate mesh topology node derivative version dof index.",err,error,*999) + DO versionIdx=1,nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%numberOfVersions + numberOfDofs=numberOfDofs+1 + nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%dofIndex(versionIdx)=numberOfDofs + ENDDO !versionIdx + ENDDO !derivativeIdx + ENDDO !nodeIdx + dofs%numberOfDofs=numberOfDofs ELSE - CALL FLAG_ERROR("Topology dofs is not assocaited",ERR,ERROR,*999) + CALL FlagError("Topology dofs is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Topology nodes is not assocaited",ERR,ERROR,*999) + CALL FlagError("Topology nodes is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Topology is not associated",ERR,ERROR,*999) + CALL FlagError("Topology is not associated.",err,error,*999) ENDIF - CALL EXITS("MESH_TOPOLOGY_DOFS_CALCULATE") + CALL Exits("MeshTopologyDofsCalculate") RETURN -999 CALL ERRORS("MESH_TOPOLOGY_DOFS_CALCULATE",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_DOFS_CALCULATE") +999 CALL Errors("MeshTopologyDofsCalculate",err,error) + CALL Exits("MeshTopologyDofsCalculate") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_DOFS_CALCULATE + + END SUBROUTINE MeshTopologyDofsCalculate ! !=============================================================================================================================== ! !>Finalises the dof data structures for a mesh topology and deallocates any memory. \todo pass in dofs - SUBROUTINE MESH_TOPOLOGY_DOFS_FINALISE(TOPOLOGY,ERR,ERROR,*) + SUBROUTINE MeshTopologyDofsFinalise(dofs,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !Initialises the dofs in a given mesh topology. \todo finalise on error - SUBROUTINE MESH_TOPOLOGY_DOFS_INITIALISE(TOPOLOGY,ERR,ERROR,*) + !>Initialises the dofs in a given mesh topology. + SUBROUTINE MeshTopologyDofsInitialise(topology,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !TOPOLOGY%MESH + ALLOCATE(topology%dofs,STAT=err) + IF(ERR/=0) CALL FlagError("Could not allocate topology dofs",err,error,*999) + topology%dofs%numberOfDofs=0 + topology%dofs%meshComponentTopology=>topology ENDIF ELSE - CALL FLAG_ERROR("Topology is not associated",ERR,ERROR,*999) + CALL FlagError("Topology is not associated",err,error,*998) ENDIF - CALL EXITS("MESH_TOPOLOGY_DOFS_INITIALISE") + CALL Exits("MeshTopologyDofsInitialise") RETURN -999 CALL ERRORS("MESH_TOPOLOGY_DOFS_INITIALISE",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_DOFS_INITIALISE") +999 CALL MeshTopologyDofsFinalise(topology%dofs,dummyErr,dummyError,*998) +998 CALL Errors("MeshTopologyDofsInitialise",err,error) + CALL Exits("MeshTopologyDofsInitialise") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_DOFS_INITIALISE + + END SUBROUTINE MeshTopologyDofsInitialise ! !================================================================================================================================ @@ -7067,12 +7204,13 @@ END SUBROUTINE MESH_TOPOLOGY_DOFS_INITIALISE SUBROUTINE MESH_TOPOLOGY_ELEMENTS_CREATE_FINISH(ELEMENTS,ERR,ERROR,*) !Argument variables - TYPE(MeshComponentElementsType), POINTER :: ELEMENTS !ELEMENTS%MESH - IF(ASSOCIATED(MESH)) THEN - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,"Number of global elements = ",MESH%NUMBER_OF_ELEMENTS, & - & ERR,ERROR,*999) - DO ne=1,MESH%NUMBER_OF_ELEMENTS - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Element = ",ne,ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Global number = ",ELEMENTS%ELEMENTS(ne)%GLOBAL_NUMBER, & - & ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," User number = ",ELEMENTS%ELEMENTS(ne)%USER_NUMBER, & + meshComponentTopology=>elements%meshComponentTopology + IF(ASSOCIATED(meshComponentTopology)) THEN + mesh=>meshComponentTopology%mesh + IF(ASSOCIATED(MESH)) THEN + CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,"Number of global elements = ",MESH%NUMBER_OF_ELEMENTS, & & ERR,ERROR,*999) - IF(ASSOCIATED(ELEMENTS%ELEMENTS(ne)%BASIS)) THEN - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Basis number = ",ELEMENTS%ELEMENTS(ne)%BASIS%USER_NUMBER, & + DO ne=1,MESH%NUMBER_OF_ELEMENTS + CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Element = ",ne,ERR,ERROR,*999) + CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Global number = ",ELEMENTS%ELEMENTS(ne)%GLOBAL_NUMBER, & & ERR,ERROR,*999) - ELSE - CALL FLAG_ERROR("Basis is not associated.",ERR,ERROR,*999) - ENDIF - IF(ALLOCATED(ELEMENTS%ELEMENTS(ne)%USER_ELEMENT_NODES)) THEN - CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,ELEMENTS%ELEMENTS(ne)% BASIS%NUMBER_OF_NODES,8,8, & - & ELEMENTS%ELEMENTS(ne)%USER_ELEMENT_NODES,'(" User element nodes =",8(X,I6))','(26X,8(X,I6))',ERR,ERROR,*999) - ELSE - CALL FLAG_ERROR("User element nodes are not associated.",ERR,ERROR,*999) - ENDIF - IF(ALLOCATED(ELEMENTS%ELEMENTS(ne)%GLOBAL_ELEMENT_NODES)) THEN - CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,ELEMENTS%ELEMENTS(ne)%BASIS%NUMBER_OF_NODES,8,8, & - & ELEMENTS%ELEMENTS(ne)%GLOBAL_ELEMENT_NODES,'(" Global element nodes =",8(X,I6))','(26X,8(X,I6))',ERR,ERROR,*999) - ELSE - CALL FLAG_ERROR("Global element nodes are not associated.",ERR,ERROR,*999) - ENDIF - ENDDO !ne + CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," User number = ",ELEMENTS%ELEMENTS(ne)%USER_NUMBER, & + & ERR,ERROR,*999) + IF(ASSOCIATED(ELEMENTS%ELEMENTS(ne)%BASIS)) THEN + CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Basis number = ",ELEMENTS%ELEMENTS(ne)%BASIS% & + & USER_NUMBER,ERR,ERROR,*999) + ELSE + CALL FLAG_ERROR("Basis is not associated.",ERR,ERROR,*999) + ENDIF + IF(ALLOCATED(ELEMENTS%ELEMENTS(ne)%USER_ELEMENT_NODES)) THEN + CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,ELEMENTS%ELEMENTS(ne)% BASIS%NUMBER_OF_NODES,8,8, & + & ELEMENTS%ELEMENTS(ne)%USER_ELEMENT_NODES,'(" User element nodes =",8(X,I6))','(26X,8(X,I6))', & + & ERR,ERROR,*999) + ELSE + CALL FLAG_ERROR("User element nodes are not associated.",ERR,ERROR,*999) + ENDIF + IF(ALLOCATED(ELEMENTS%ELEMENTS(ne)%GLOBAL_ELEMENT_NODES)) THEN + CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1,ELEMENTS%ELEMENTS(ne)%BASIS%NUMBER_OF_NODES,8,8, & + & ELEMENTS%ELEMENTS(ne)%GLOBAL_ELEMENT_NODES,'(" Global element nodes =",8(X,I6))','(26X,8(X,I6))', & + & ERR,ERROR,*999) + ELSE + CALL FLAG_ERROR("Global element nodes are not associated.",ERR,ERROR,*999) + ENDIF + ENDDO !ne + ELSE + CALL FLAG_ERROR("Mesh component topology mesh is not associated.",ERR,ERROR,*999) + ENDIF ELSE - CALL FLAG_ERROR("Mesh elements mesh is not associated.",ERR,ERROR,*999) + CALL FLAG_ERROR("Mesh elements mesh component topology is not associated.",ERR,ERROR,*999) ENDIF ENDIF - + CALL EXITS("MESH_TOPOLOGY_ELEMENTS_CREATE_FINISH") RETURN 999 CALL ERRORS("MESH_TOPOLOGY_ELEMENTS_CREATE_FINISH",ERR,ERROR) @@ -7140,7 +7285,7 @@ SUBROUTINE MESH_TOPOLOGY_ELEMENTS_CREATE_START(MESH,MESH_COMPONENT_NUMBER,BASIS, TYPE(MESH_TYPE), POINTER :: MESH !=1.AND.GLOBAL_NUMBER<=ELEMENTS%NUMBER_OF_ELEMENTS) THEN IF(SIZE(USER_ELEMENT_NODES,1)==ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS%NUMBER_OF_NODES) THEN - IF(ASSOCIATED(ELEMENTS%MESH)) THEN - REGION=>ELEMENTS%MESH%REGION - IF(ASSOCIATED(REGION)) THEN - NODES=>REGION%NODES - ELSE - INTERFACE=>ELEMENTS%MESH%INTERFACE - IF(ASSOCIATED(INTERFACE)) THEN - NODES=>INTERFACE%NODES - PARENT_REGION=>INTERFACE%PARENT_REGION - IF(.NOT.ASSOCIATED(PARENT_REGION)) CALL FLAG_ERROR("Mesh interface has no parent region.",ERR,ERROR,*999) + meshComponentTopology=>elements%meshComponentTopology + IF(ASSOCIATED(meshComponentTopology)) THEN + mesh=>meshComponentTopology%mesh + IF(ASSOCIATED(mesh)) THEN + REGION=>MESH%REGION + IF(ASSOCIATED(REGION)) THEN + NODES=>REGION%NODES ELSE - CALL FLAG_ERROR("Elements mesh has no associated region or interface.",ERR,ERROR,*999) - ENDIF - ENDIF - IF(ASSOCIATED(NODES)) THEN - ELEMENT_NODES_OK=.TRUE. - ALLOCATE(GLOBAL_ELEMENT_NODES(ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS%NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate global element nodes.",ERR,ERROR,*999) - ALLOCATE(BAD_NODES(ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS%NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate bad nodes.",ERR,ERROR,*999) - NUMBER_OF_BAD_NODES=0 - DO nn=1,ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS%NUMBER_OF_NODES - CALL NODE_CHECK_EXISTS(NODES,USER_ELEMENT_NODES(nn),NODE_EXISTS,GLOBAL_NODE_NUMBER,ERR,ERROR,*999) - IF(NODE_EXISTS) THEN - GLOBAL_ELEMENT_NODES(nn)=GLOBAL_NODE_NUMBER + INTERFACE=>MESH%INTERFACE + IF(ASSOCIATED(INTERFACE)) THEN + NODES=>INTERFACE%NODES + PARENT_REGION=>INTERFACE%PARENT_REGION + IF(.NOT.ASSOCIATED(PARENT_REGION)) CALL FLAG_ERROR("Mesh interface has no parent region.",ERR,ERROR,*999) ELSE - NUMBER_OF_BAD_NODES=NUMBER_OF_BAD_NODES+1 - BAD_NODES(NUMBER_OF_BAD_NODES)=USER_ELEMENT_NODES(nn) - ELEMENT_NODES_OK=.FALSE. + CALL FLAG_ERROR("Elements mesh has no associated region or interface.",ERR,ERROR,*999) ENDIF - ENDDO !nn - IF(ELEMENT_NODES_OK) THEN - ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%USER_ELEMENT_NODES=USER_ELEMENT_NODES - ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%GLOBAL_ELEMENT_NODES=GLOBAL_ELEMENT_NODES - ELSE - IF(NUMBER_OF_BAD_NODES==1) THEN - IF(ASSOCIATED(REGION)) THEN - LOCAL_ERROR="The element user node number of "//TRIM(NUMBER_TO_VSTRING(BAD_NODES(1),"*",ERR,ERROR))// & - & " is not defined in region "//TRIM(NUMBER_TO_VSTRING(REGION%USER_NUMBER,"*",ERR,ERROR))//"." + ENDIF + IF(ASSOCIATED(NODES)) THEN + ELEMENT_NODES_OK=.TRUE. + ALLOCATE(GLOBAL_ELEMENT_NODES(ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS%NUMBER_OF_NODES),STAT=ERR) + IF(ERR/=0) CALL FLAG_ERROR("Could not allocate global element nodes.",ERR,ERROR,*999) + ALLOCATE(BAD_NODES(ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS%NUMBER_OF_NODES),STAT=ERR) + IF(ERR/=0) CALL FLAG_ERROR("Could not allocate bad nodes.",ERR,ERROR,*999) + NUMBER_OF_BAD_NODES=0 + DO nn=1,ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS%NUMBER_OF_NODES + CALL NODE_CHECK_EXISTS(NODES,USER_ELEMENT_NODES(nn),NODE_EXISTS,GLOBAL_NODE_NUMBER,ERR,ERROR,*999) + IF(NODE_EXISTS) THEN + GLOBAL_ELEMENT_NODES(nn)=GLOBAL_NODE_NUMBER ELSE - LOCAL_ERROR="The element user node number of "//TRIM(NUMBER_TO_VSTRING(BAD_NODES(1),"*",ERR,ERROR))// & - & " is not defined in interface number "//TRIM(NUMBER_TO_VSTRING(INTERFACE%USER_NUMBER,"*",ERR,ERROR))// & - & " of parent region number "//TRIM(NUMBER_TO_VSTRING(PARENT_REGION%USER_NUMBER,"*",ERR,ERROR))//"." + NUMBER_OF_BAD_NODES=NUMBER_OF_BAD_NODES+1 + BAD_NODES(NUMBER_OF_BAD_NODES)=USER_ELEMENT_NODES(nn) + ELEMENT_NODES_OK=.FALSE. ENDIF + ENDDO !nn + IF(ELEMENT_NODES_OK) THEN + ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%USER_ELEMENT_NODES=USER_ELEMENT_NODES + ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%GLOBAL_ELEMENT_NODES=GLOBAL_ELEMENT_NODES ELSE - LOCAL_ERROR="The element user node number of "//TRIM(NUMBER_TO_VSTRING(BAD_NODES(1),"*",ERR,ERROR)) - DO nn=2,NUMBER_OF_BAD_NODES-1 - LOCAL_ERROR=LOCAL_ERROR//","//TRIM(NUMBER_TO_VSTRING(BAD_NODES(nn),"*",ERR,ERROR)) - ENDDO !nn - IF(ASSOCIATED(REGION)) THEN - LOCAL_ERROR=LOCAL_ERROR//" & "//TRIM(NUMBER_TO_VSTRING(BAD_NODES(NUMBER_OF_BAD_NODES),"*",ERR,ERROR))// & - & " are not defined in region number "//TRIM(NUMBER_TO_VSTRING(REGION%USER_NUMBER,"*",ERR,ERROR))//"." + IF(NUMBER_OF_BAD_NODES==1) THEN + IF(ASSOCIATED(REGION)) THEN + LOCAL_ERROR="The element user node number of "//TRIM(NUMBER_TO_VSTRING(BAD_NODES(1),"*",ERR,ERROR))// & + & " is not defined in region "//TRIM(NUMBER_TO_VSTRING(REGION%USER_NUMBER,"*",ERR,ERROR))//"." + ELSE + LOCAL_ERROR="The element user node number of "//TRIM(NUMBER_TO_VSTRING(BAD_NODES(1),"*",ERR,ERROR))// & + & " is not defined in interface number "// & + & TRIM(NUMBER_TO_VSTRING(INTERFACE%USER_NUMBER,"*",ERR,ERROR))// & + & " of parent region number "//TRIM(NUMBER_TO_VSTRING(PARENT_REGION%USER_NUMBER,"*",ERR,ERROR))//"." + ENDIF ELSE - LOCAL_ERROR=LOCAL_ERROR//" & "//TRIM(NUMBER_TO_VSTRING(BAD_NODES(NUMBER_OF_BAD_NODES),"*",ERR,ERROR))// & - & " are not defined in interface number "// & - & TRIM(NUMBER_TO_VSTRING(INTERFACE%USER_NUMBER,"*",ERR,ERROR))//" of parent region number "// & - & TRIM(NUMBER_TO_VSTRING(PARENT_REGION%USER_NUMBER,"*",ERR,ERROR))//"." + LOCAL_ERROR="The element user node number of "//TRIM(NUMBER_TO_VSTRING(BAD_NODES(1),"*",ERR,ERROR)) + DO nn=2,NUMBER_OF_BAD_NODES-1 + LOCAL_ERROR=LOCAL_ERROR//","//TRIM(NUMBER_TO_VSTRING(BAD_NODES(nn),"*",ERR,ERROR)) + ENDDO !nn + IF(ASSOCIATED(REGION)) THEN + LOCAL_ERROR=LOCAL_ERROR//" & "//TRIM(NUMBER_TO_VSTRING(BAD_NODES(NUMBER_OF_BAD_NODES),"*",ERR,ERROR))// & + & " are not defined in region number "//TRIM(NUMBER_TO_VSTRING(REGION%USER_NUMBER,"*",ERR,ERROR))//"." + ELSE + LOCAL_ERROR=LOCAL_ERROR//" & "//TRIM(NUMBER_TO_VSTRING(BAD_NODES(NUMBER_OF_BAD_NODES),"*",ERR,ERROR))// & + & " are not defined in interface number "// & + & TRIM(NUMBER_TO_VSTRING(INTERFACE%USER_NUMBER,"*",ERR,ERROR))//" of parent region number "// & + & TRIM(NUMBER_TO_VSTRING(PARENT_REGION%USER_NUMBER,"*",ERR,ERROR))//"." + ENDIF ENDIF + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) ENDIF - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - ENDIF - ELSE - IF(ASSOCIATED(REGION)) THEN - - CALL FLAG_ERROR("The elements mesh region does not have any associated nodes.",ERR,ERROR,*999) ELSE - CALL FLAG_ERROR("The elements mesh interface does not have any associated nodes.",ERR,ERROR,*999) + IF(ASSOCIATED(REGION)) THEN + CALL FLAG_ERROR("The elements mesh region does not have any associated nodes.",ERR,ERROR,*999) + ELSE + CALL FLAG_ERROR("The elements mesh interface does not have any associated nodes.",ERR,ERROR,*999) + ENDIF ENDIF + ELSE + CALL FLAG_ERROR("The mesh component topology mesh is not associated.",ERR,ERROR,*999) ENDIF ELSE - CALL FLAG_ERROR("The elements do not have an associated mesh.",ERR,ERROR,*999) + CALL FLAG_ERROR("The elements mesh component topology is not associated.",ERR,ERROR,*999) ENDIF ELSE CALL FLAG_ERROR("Number of element nodes does not match number of basis nodes for this element.",ERR,ERROR,*999) @@ -7724,7 +7877,7 @@ SUBROUTINE MESH_TOPOLOGY_ELEMENTS_ELEMENT_NODE_VERSION_SET(GLOBAL_NUMBER,ELEMENT !Argument variables INTEGER(INTG), INTENT(IN) :: GLOBAL_NUMBER !=1.AND.GLOBAL_NUMBER<=ELEMENTS%NUMBER_OF_ELEMENTS) THEN IF(USER_ELEMENT_NODE_INDEX>=1.AND.USER_ELEMENT_NODE_INDEX<=ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS%NUMBER_OF_NODES) THEN - IF(ASSOCIATED(ELEMENTS%MESH)) THEN - REGION=>ELEMENTS%MESH%REGION - IF(ASSOCIATED(REGION)) THEN - NODES=>REGION%NODES - ELSE - INTERFACE=>ELEMENTS%MESH%INTERFACE - IF(ASSOCIATED(INTERFACE)) THEN - NODES=>INTERFACE%NODES - PARENT_REGION=>INTERFACE%PARENT_REGION - IF(.NOT.ASSOCIATED(PARENT_REGION)) CALL FLAG_ERROR("Mesh interface has no parent region.",ERR,ERROR,*999) + meshComponentTopology=>elements%meshComponentTopology + IF(ASSOCIATED(meshComponentTopology)) THEN + mesh=>meshComponentTopology%mesh + IF(ASSOCIATED(mesh)) THEN + REGION=>mesh%REGION + IF(ASSOCIATED(REGION)) THEN + NODES=>REGION%NODES ELSE - CALL FLAG_ERROR("Elements mesh has no associated region or interface.",ERR,ERROR,*999) - ENDIF - ENDIF - IF(DERIVATIVE_NUMBER>=1.AND.DERIVATIVE_NUMBER<=ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS% & - & NUMBER_OF_DERIVATIVES(USER_ELEMENT_NODE_INDEX)) THEN !Check if the specified derivative exists - IF(VERSION_NUMBER>=1) THEN !Check if the specified version is greater than 1 - ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%USER_ELEMENT_NODE_VERSIONS(DERIVATIVE_NUMBER,USER_ELEMENT_NODE_INDEX) & - & = VERSION_NUMBER - !\todo : There is redunancy in USER_ELEMENT_NODE_VERSIONS since it was allocated in MESH_TOPOLOGY_ELEMENTS_CREATE_START based on MAXIMUM_NUMBER_OF_DERIVATIVES for that elements basis:ALLOCATE(ELEMENTS%ELEMENTS(ne)%USER_ELEMENT_NODE_VERSIONS(BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES,BASIS%NUMBER_OF_NODES),STAT=ERR) + INTERFACE=>mesh%INTERFACE + IF(ASSOCIATED(INTERFACE)) THEN + NODES=>INTERFACE%NODES + PARENT_REGION=>INTERFACE%PARENT_REGION + IF(.NOT.ASSOCIATED(PARENT_REGION)) CALL FLAG_ERROR("Mesh interface has no parent region.",ERR,ERROR,*999) + ELSE + CALL FLAG_ERROR("Elements mesh has no associated region or interface.",ERR,ERROR,*999) + ENDIF + ENDIF + IF(DERIVATIVE_NUMBER>=1.AND.DERIVATIVE_NUMBER<=ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS% & + & NUMBER_OF_DERIVATIVES(USER_ELEMENT_NODE_INDEX)) THEN !Check if the specified derivative exists + IF(VERSION_NUMBER>=1) THEN !Check if the specified version is greater than 1 + ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%USER_ELEMENT_NODE_VERSIONS(DERIVATIVE_NUMBER,USER_ELEMENT_NODE_INDEX) & + & = VERSION_NUMBER + !\todo : There is redunancy in USER_ELEMENT_NODE_VERSIONS since it was allocated in MESH_TOPOLOGY_ELEMENTS_CREATE_START based on MAXIMUM_NUMBER_OF_DERIVATIVES for that elements basis:ALLOCATE(ELEMENTS%ELEMENTS(ne)%USER_ELEMENT_NODE_VERSIONS(BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES,BASIS%NUMBER_OF_NODES),STAT=ERR) + ELSE + LOCAL_ERROR="The specified node version number of "//TRIM(NUMBER_TO_VSTRING(VERSION_NUMBER,"*", & + & ERR,ERROR))//" is invalid. The element node index should be greater than 1." + CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + ENDIF ELSE - LOCAL_ERROR="The specified node version number of "//TRIM(NUMBER_TO_VSTRING(VERSION_NUMBER,"*", & - & ERR,ERROR))//" is invalid. The element node index should be greater than 1." + LOCAL_ERROR="The specified node derivative number of "//TRIM(NUMBER_TO_VSTRING(DERIVATIVE_NUMBER,"*", & + & ERR,ERROR))//" is invalid. The element node derivative index should be between 1 and "// & + & TRIM(NUMBER_TO_VSTRING(ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS%NUMBER_OF_DERIVATIVES( & + & USER_ELEMENT_NODE_INDEX),"*",ERR,ERROR))//"." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) ENDIF ELSE - LOCAL_ERROR="The specified node derivative number of "//TRIM(NUMBER_TO_VSTRING(DERIVATIVE_NUMBER,"*", & - & ERR,ERROR))//" is invalid. The element node derivative index should be between 1 and "// & - & TRIM(NUMBER_TO_VSTRING(ELEMENTS%ELEMENTS(GLOBAL_NUMBER)%BASIS%NUMBER_OF_DERIVATIVES(USER_ELEMENT_NODE_INDEX), & - & "*",ERR,ERROR))//"." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + CALL FlagError("The mesh component topology mesh is not associated.",err,error,*999) ENDIF + ELSE + CALL FLAG_ERROR("The elements mesh component topology is not associated.",err,error,*999) ENDIF ELSE LOCAL_ERROR="The specified element node index of "//TRIM(NUMBER_TO_VSTRING(USER_ELEMENT_NODE_INDEX,"*",ERR,ERROR))// & @@ -7837,7 +8000,7 @@ SUBROUTINE MESH_TOPOLOGY_ELEMENTS_ADJACENT_ELEMENTS_CALCULATE(TOPOLOGY,ERR,ERROR IF(ASSOCIATED(TOPOLOGY%ELEMENTS)) THEN !Loop over the global elements in the mesh DO ne=1,TOPOLOGY%ELEMENTS%NUMBER_OF_ELEMENTS - !%%%% first we initialize lists that are required to find the adjacent elements list + !%%%% first we initialize lists that are required to find the adjacent elements list BASIS=>TOPOLOGY%ELEMENTS%ELEMENTS(ne)%BASIS DO nic=-BASIS%NUMBER_OF_XI_COORDINATES,BASIS%NUMBER_OF_XI_COORDINATES NULLIFY(ADJACENT_ELEMENTS_LIST(nic)%PTR) @@ -7935,8 +8098,8 @@ SUBROUTINE MESH_TOPOLOGY_ELEMENTS_ADJACENT_ELEMENTS_CALCULATE(TOPOLOGY,ERR,ERROR IF(NUMBER_NODE_MATCHES>0) THEN !Find list of elements surrounding those nodes np1=NODE_MATCHES(1) - DO nep1=1,TOPOLOGY%NODES%NODES(np1)%NUMBER_OF_SURROUNDING_ELEMENTS - ne1=TOPOLOGY%NODES%NODES(np1)%SURROUNDING_ELEMENTS(nep1) + DO nep1=1,TOPOLOGY%NODES%NODES(np1)%numberOfSurroundingElements + ne1=TOPOLOGY%NODES%NODES(np1)%surroundingElements(nep1) IF(ne1/=ne) THEN !Don't want the current element ! grab the nodes list for current and this surrouding elements ! current face : NODE_MATCHES @@ -7990,8 +8153,8 @@ SUBROUTINE MESH_TOPOLOGY_ELEMENTS_ADJACENT_ELEMENTS_CALCULATE(TOPOLOGY,ERR,ERROR !Find list of elements surrounding those nodes DO node_idx=1,NUMBER_NODE_MATCHES np1=NODE_MATCHES(node_idx) - DO nep1=1,TOPOLOGY%NODES%NODES(np1)%NUMBER_OF_SURROUNDING_ELEMENTS - ne1=TOPOLOGY%NODES%NODES(np1)%SURROUNDING_ELEMENTS(nep1) + DO nep1=1,TOPOLOGY%NODES%NODES(np1)%numberOfSurroundingElements + ne1=TOPOLOGY%NODES%NODES(np1)%surroundingElements(nep1) IF(ne1/=ne) THEN !Don't want the current element ! grab the nodes list for current and this surrouding elements ! current face : NODE_MATCHES @@ -8091,7 +8254,7 @@ END SUBROUTINE MESH_TOPOLOGY_ELEMENTS_ADJACENT_ELEMENTS_CALCULATE SUBROUTINE MESH_TOPOLOGY_ELEMENTS_FINALISE(ELEMENTS,ERR,ERROR,*) !Argument variables - TYPE(MeshComponentElementsType), POINTER :: ELEMENTS !TOPOLOGY%MESH + TOPOLOGY%ELEMENTS%meshComponentTopology=>TOPOLOGY NULLIFY(TOPOLOGY%ELEMENTS%ELEMENTS) NULLIFY(TOPOLOGY%ELEMENTS%ELEMENTS_TREE) ENDIF @@ -8174,7 +8337,7 @@ SUBROUTINE MESH_TOPOLOGY_DATA_POINTS_INITIALISE(TOPOLOGY,ERR,ERROR,*) ALLOCATE(TOPOLOGY%dataPoints,STAT=ERR) IF(ERR/=0) CALL FLAG_ERROR("Could not allocate topology data points",ERR,ERROR,*999) TOPOLOGY%dataPoints%totalNumberOfProjectedData=0 - TOPOLOGY%dataPoints%mesh=>TOPOLOGY%MESH + TOPOLOGY%dataPoints%meshComponentTopology=>topology ENDIF ELSE CALL FLAG_ERROR("Topology is not associated",ERR,ERROR,*999) @@ -8199,7 +8362,7 @@ SUBROUTINE MESH_TOPOLOGY_ELEMENTS_NUMBER_GET(GLOBAL_NUMBER,USER_NUMBER,ELEMENTS, !Argument variables INTEGER(INTG), INTENT(IN) :: GLOBAL_NUMBER !Changes/sets the user numbers for all elements. - SUBROUTINE Mesh_TopologyElementsUserNumbersAllSet(elements,userNumbers,err,error,*) + SUBROUTINE MeshTopologyElementsUserNumbersAllSet(elements,userNumbers,err,error,*) !Argument variables - TYPE(MeshComponentElementsType), POINTER :: elements !Calculates the data points in the given mesh topology. - SUBROUTINE Mesh_TopologyDataPointsCalculateProjection(mesh,dataProjection,err,error,*) + SUBROUTINE MeshTopologyDataPointsCalculateProjection(mesh,dataProjection,err,error,*) !Argument variables TYPE(MESH_TYPE), POINTER :: mesh !Finalises the topology in the given mesh. \todo pass in the mesh topology - SUBROUTINE MESH_TOPOLOGY_FINALISE(MESH,ERR,ERROR,*) + SUBROUTINE MeshTopologyFinalise(mesh,err,error,*) !Argument variables - TYPE(MESH_TYPE), POINTER :: MESH !Finalises the topology in the given mesh. + SUBROUTINE MeshTopologyComponentFinalise(meshComponent,err,error,*) + + !Argument variables + TYPE(MeshComponentTopologyType), POINTER :: meshComponent !Initialises the topology for a given mesh. \todo finalise on error - SUBROUTINE MESH_TOPOLOGY_INITIALISE(MESH,ERR,ERROR,*) + SUBROUTINE MeshTopologyInitialise(mesh,err,error,*) !Argument variables - TYPE(MESH_TYPE), POINTER :: MESH !MESH - NULLIFY(MESH%TOPOLOGY(component_idx)%PTR%ELEMENTS) - NULLIFY(MESH%TOPOLOGY(component_idx)%PTR%NODES) - NULLIFY(MESH%TOPOLOGY(component_idx)%PTR%DOFS) - NULLIFY(MESH%TOPOLOGY(component_idx)%PTR%dataPoints) + ALLOCATE(mesh%topology(mesh%NUMBER_OF_COMPONENTS),STAT=err) + IF(err/=0) CALL FlagError("Mesh topology could not be allocated.",err,error,*999) + DO componentIdx=1,mesh%NUMBER_OF_COMPONENTS + ALLOCATE(mesh%topology(componentIdx)%ptr,STAT=err) + IF(err/=0) CALL FlagError("Mesh topology component could not be allocated.",err,error,*999) + mesh%topology(componentIdx)%ptr%mesh=>mesh + NULLIFY(mesh%topology(componentIdx)%ptr%elements) + NULLIFY(mesh%topology(componentIdx)%ptr%nodes) + NULLIFY(mesh%topology(componentIdx)%ptr%dofs) + NULLIFY(mesh%topology(componentIdx)%ptr%dataPoints) !Initialise the topology components - CALL MESH_TOPOLOGY_ELEMENTS_INITIALISE(MESH%TOPOLOGY(component_idx)%PTR,ERR,ERROR,*999) - CALL MESH_TOPOLOGY_NODES_INITIALISE(MESH%TOPOLOGY(component_idx)%PTR,ERR,ERROR,*999) - CALL MESH_TOPOLOGY_DOFS_INITIALISE(MESH%TOPOLOGY(component_idx)%PTR,ERR,ERROR,*999) - CALL MESH_TOPOLOGY_DATA_POINTS_INITIALISE(MESH%TOPOLOGY(component_idx)%PTR,ERR,ERROR,*999) - ENDDO !component_idx + CALL MESH_TOPOLOGY_ELEMENTS_INITIALISE(mesh%topology(componentIdx)%ptr,err,error,*999) + CALL MeshTopologyNodesInitialise(mesh%topology(componentIdx)%ptr,err,error,*999) + CALL MeshTopologyDofsInitialise(mesh%topology(componentIdx)%ptr,err,error,*999) + CALL MESH_TOPOLOGY_DATA_POINTS_INITIALISE(mesh%topology(componentIdx)%ptr,err,error,*999) + ENDDO !componentIdx ENDIF ELSE - CALL FLAG_ERROR("Mesh is not associated",ERR,ERROR,*999) + CALL FlagError("Mesh is not associated.",err,error,*999) ENDIF - CALL EXITS("MESH_TOPOLOGY_INITIALISE") + CALL EXITS("MeshTopologyInitialise") RETURN -999 CALL ERRORS("MESH_TOPOLOGY_INITIALISE",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_INITIALISE") +999 CALL ERRORS("MeshTopologyInitialise",err,error) + CALL EXITS("MeshTopologyInitialise") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_INITIALISE + END SUBROUTINE MeshTopologyInitialise ! !================================================================================================================================ ! !>Checks that a user element number exists in a mesh component. - SUBROUTINE MESH_TOPOLOGY_ELEMENT_CHECK_EXISTS(MESH,MESH_COMPONENT_NUMBER,USER_ELEMENT_NUMBER,ELEMENT_EXISTS, & - & MESH_GLOBAL_ELEMENT_NUMBER,ERR,ERROR,*) + SUBROUTINE MeshTopologyElementCheckExistsMesh(mesh,meshComponentNumber,userElementNumber,elementExists,globalElementNumber, & + & err,error,*) !Argument variables - TYPE(MESH_TYPE), POINTER :: MESH !=1.AND.MESH_COMPONENT_NUMBER<=MESH%NUMBER_OF_COMPONENTS) THEN - IF(ASSOCIATED(MESH%TOPOLOGY)) THEN - MESH_TOPOLOGY=>MESH%TOPOLOGY(MESH_COMPONENT_NUMBER)%PTR - IF(ASSOCIATED(MESH_TOPOLOGY)) THEN - MESH_ELEMENTS=>MESH_TOPOLOGY%ELEMENTS - IF(ASSOCIATED(MESH_ELEMENTS)) THEN - NULLIFY(TREE_NODE) - CALL TREE_SEARCH(MESH_ELEMENTS%ELEMENTS_TREE,USER_ELEMENT_NUMBER,TREE_NODE,ERR,ERROR,*999) - IF(ASSOCIATED(TREE_NODE)) THEN - CALL TREE_NODE_VALUE_GET(MESH_ELEMENTS%ELEMENTS_TREE,TREE_NODE,MESH_GLOBAL_ELEMENT_NUMBER,ERR,ERROR,*999) - ELEMENT_EXISTS=.TRUE. - ENDIF - ELSE - CALL FLAG_ERROR("Mesh topology elements is not associated.",ERR,ERROR,*999) - ENDIF - ELSE - LOCAL_ERROR="Mesh topology is not associated for component number "// & - & TRIM(NUMBER_TO_VSTRING(MESH_COMPONENT_NUMBER,"*",ERR,ERROR))//"." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - ENDIF - ELSE - CALL FLAG_ERROR("Mesh topology is not associated.",ERR,ERROR,*999) - ENDIF - ELSE - LOCAL_ERROR="The specified mesh component of "//TRIM(NUMBER_TO_VSTRING(MESH_COMPONENT_NUMBER,"*",ERR,ERROR))// & - & " is invalid. Mesh number "//TRIM(NUMBER_TO_VSTRING(MESH%USER_NUMBER,"*",ERR,ERROR))//" has "// & - & TRIM(NUMBER_TO_VSTRING(MESH%NUMBER_OF_COMPONENTS,"*",ERR,ERROR))//"." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - ENDIF + IF(ASSOCIATED(mesh)) THEN + IF(mesh%MESH_FINISHED) THEN + CALL MESH_TOPOLOGY_ELEMENTS_GET(mesh,meshComponentNumber,elements,err,error,*999) + CALL MeshTopologyElementCheckExistsMeshElements(elements,userElementNumber,elementExists,globalElementNumber,err,error,*999) ELSE - CALL FLAG_ERROR("Mesh has not been finished.",ERR,ERROR,*999) + CALL FlagError("Mesh has not been finished.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh is not associated.",ERR,ERROR,*999) + CALL FlagError("Mesh is not associated.",err,error,*999) ENDIF - CALL EXITS("MESH_TOPOLOGY_ELEMENT_CHECK_EXISTS") + CALL Exits("MeshTopologyElementCheckExistsMesh") RETURN -999 CALL ERRORS("MESH_TOPOLOGY_ELEMENT_CHECK_EXISTS",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_ELEMENT_CHECK_EXISTS") +999 CALL Errors("MeshTopologyElementCheckExistsMesh",err,error) + CALL Exits("MeshTopologyElementCheckExistsMesh") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_ELEMENT_CHECK_EXISTS + END SUBROUTINE MeshTopologyElementCheckExistsMesh ! !================================================================================================================================ ! - !>Checks that a user node number exists in a mesh component. - SUBROUTINE MESH_TOPOLOGY_NODE_CHECK_EXISTS(MESH,MESH_COMPONENT_NUMBER,USER_NODE_NUMBER,NODE_EXISTS, & - & MESH_GLOBAL_NODE_NUMBER,ERR,ERROR,*) + !>Checks that a user element number exists in a mesh elements. + SUBROUTINE MeshTopologyElementCheckExistsMeshElements(meshElements,userElementNumber,elementExists,globalElementNumber, & + & err,error,*) !Argument variables - TYPE(MESH_TYPE), POINTER :: MESH !=1.AND.MESH_COMPONENT_NUMBER<=MESH%NUMBER_OF_COMPONENTS) THEN - IF(ASSOCIATED(MESH%TOPOLOGY)) THEN - MESH_TOPOLOGY=>MESH%TOPOLOGY(MESH_COMPONENT_NUMBER)%PTR - IF(ASSOCIATED(MESH_TOPOLOGY)) THEN - MESH_NODES=>MESH_TOPOLOGY%NODES - IF(ASSOCIATED(MESH_NODES)) THEN - NULLIFY(TREE_NODE) - CALL TREE_SEARCH(MESH_NODES%NODES_TREE,USER_NODE_NUMBER,TREE_NODE,ERR,ERROR,*999) - IF(ASSOCIATED(TREE_NODE)) THEN - CALL TREE_NODE_VALUE_GET(MESH_NODES%NODES_TREE,TREE_NODE,MESH_GLOBAL_NODE_NUMBER,ERR,ERROR,*999) - NODE_EXISTS=.TRUE. - ENDIF - ELSE - CALL FLAG_ERROR("Mesh topology nodes is not associated.",ERR,ERROR,*999) + elementExists=.FALSE. + globalElementNumber=0 + IF(ASSOCIATED(meshElements)) THEN + NULLIFY(treeNode) + CALL TREE_SEARCH(meshElements%ELEMENTS_TREE,userElementNumber,treeNode,err,error,*999) + IF(ASSOCIATED(treeNode)) THEN + CALL TREE_NODE_VALUE_GET(meshElements%ELEMENTS_TREE,treeNode,globalElementNumber,err,error,*999) + elementExists=.TRUE. + ENDIF + ELSE + CALL FlagError("Mesh elements is not associated.",err,error,*999) + ENDIF + + CALL Exits("MeshTopologyElementCheckExistsMeshElements") + RETURN +999 CALL Errors("MeshTopologyElementCheckExistsMeshElements",err,error) + CALL Exits("MeshTopologyElementCheckExistsMeshElements") + RETURN 1 + + END SUBROUTINE MeshTopologyElementCheckExistsMeshElements + + ! + !================================================================================================================================ + ! + + !>Checks that a user node number exists in a mesh component. + SUBROUTINE MeshTopologyNodeCheckExistsMesh(mesh,meshComponentNumber,userNodeNumber,nodeExists,meshNodeNumber,err,error,*) + + !Argument variables + TYPE(MESH_TYPE), POINTER :: mesh !region%nodes + IF(ASSOCIATED(nodes)) THEN + CALL NODE_CHECK_EXISTS(nodes,userNodeNumber,nodeExists,globalNodeNumber,err,error,*999) + NULLIFY(treeNode) + CALL TREE_SEARCH(meshNodes%nodesTree,globalNodeNumber,treeNode,err,error,*999) + IF(ASSOCIATED(treeNode)) THEN + CALL TREE_NODE_VALUE_GET(meshNodes%nodesTree,treeNode,meshNodeNumber,err,error,*999) + nodeExists=.TRUE. + ENDIF + ELSE + CALL FlagError("Region nodes is not associated.",err,error,*999) + ENDIF + ELSE + CALL FlagError("Mesh has not been finished.",err,error,*999) + ENDIF + ELSE + CALL FlagError("Mesh is not associated.",err,error,*999) + ENDIF + + CALL Exits("MeshTopologyNodeCheckExistsMesh") + RETURN +999 CALL Errors("MeshTopologyNodeCheckExistsMesh",err,error) + CALL Exits("MeshTopologyNodeCheckExistsMesh") + RETURN 1 + + END SUBROUTINE MeshTopologyNodeCheckExistsMesh + + ! + !================================================================================================================================ + ! + + !>Checks that a user node number exists in a mesh nodes. + SUBROUTINE MeshTopologyNodeCheckExistsMeshNodes(meshNodes,userNodeNumber,nodeExists,meshNodeNumber,err,error,*) + + !Argument variables + TYPE(MeshNodesType), POINTER :: meshNodes !meshNodes%meshComponentTopology + IF(ASSOCIATED(meshComponentTopology)) THEN + mesh=>meshComponentTopology%mesh + IF(ASSOCIATED(mesh)) THEN + IF(mesh%MESH_FINISHED) THEN + CALL MeshRegionGet(mesh,region,err,error,*999) + nodes=>region%nodes + IF(ASSOCIATED(nodes)) THEN + CALL NODE_CHECK_EXISTS(nodes,userNodeNumber,nodeExists,globalNodeNumber,err,error,*999) + NULLIFY(treeNode) + CALL TREE_SEARCH(meshNodes%nodesTree,globalNodeNumber,treeNode,err,error,*999) + IF(ASSOCIATED(treeNode)) THEN + CALL TREE_NODE_VALUE_GET(meshNodes%nodesTree,treeNode,meshNodeNumber,err,error,*999) + nodeExists=.TRUE. ENDIF ELSE - LOCAL_ERROR="Mesh topology is not associated for component number "// & - & TRIM(NUMBER_TO_VSTRING(MESH_COMPONENT_NUMBER,"*",ERR,ERROR))//"." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + CALL FlagError("Region nodes is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology is not associated.",ERR,ERROR,*999) + CALL FlagError("Mesh has not been finished.",err,error,*999) ENDIF ELSE - LOCAL_ERROR="The specified mesh component of "//TRIM(NUMBER_TO_VSTRING(MESH_COMPONENT_NUMBER,"*",ERR,ERROR))// & - & " is invalid. Mesh number "//TRIM(NUMBER_TO_VSTRING(MESH%USER_NUMBER,"*",ERR,ERROR))//" has "// & - & TRIM(NUMBER_TO_VSTRING(MESH%NUMBER_OF_COMPONENTS,"*",ERR,ERROR))//"." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + CALL FlagError("Mesh component topology mesh is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh has not been finished.",ERR,ERROR,*999) + CALL FlagError("Mesh nodes mesh component topology is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh is not associated.",ERR,ERROR,*999) + CALL FLAG_ERROR("Mesh nodes is not associated.",err,error,*999) ENDIF - CALL EXITS("MESH_TOPOLOGY_NODE_CHECK_EXISTS") + CALL Exits("MeshTopologyNodeCheckExistsMeshNodes") RETURN -999 CALL ERRORS("MESH_TOPOLOGY_NODE_CHECK_EXISTS",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_NODE_CHECK_EXISTS") +999 CALL Errors("MeshTopologyNodeCheckExistsMeshNodes",err,error) + CALL Exits("MeshTopologyNodeCheckExistsMeshNodes") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_NODE_CHECK_EXISTS + END SUBROUTINE MeshTopologyNodeCheckExistsMeshNodes ! !================================================================================================================================ ! !>Finalises the given mesh topology node. - SUBROUTINE MESH_TOPOLOGY_NODE_FINALISE(NODE,ERR,ERROR,*) + SUBROUTINE MeshTopologyNodeFinalise(node,err,error,*) !Argument variables - TYPE(MESH_NODE_TYPE) :: NODE !Initialises the given mesh topology node. - SUBROUTINE MESH_TOPOLOGY_NODE_INITIALISE(NODE,ERR,ERROR,*) + SUBROUTINE MeshTopologyNodeInitialise(node,err,error,*) !Argument variables - TYPE(MESH_NODE_TYPE) :: NODE !Calculates the nodes used the mesh identified by a given mesh topology. - SUBROUTINE MESH_TOPOLOGY_NODES_CALCULATE(TOPOLOGY,ERR,ERROR,*) + SUBROUTINE MeshTopologyNodesCalculate(topology,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !TOPOLOGY%ELEMENTS - IF(ASSOCIATED(ELEMENTS)) THEN - MESHNODES=>TOPOLOGY%NODES - IF(ASSOCIATED(MESHNODES)) THEN - MESH=>TOPOLOGY%MESH - IF(ASSOCIATED(MESH)) THEN - NULLIFY(INTERFACE) - REGION=>MESH%REGION - IF(ASSOCIATED(REGION)) THEN - NODES=>REGION%NODES - ELSE - INTERFACE=>MESH%INTERFACE - IF(ASSOCIATED(INTERFACE)) THEN - NODES=>INTERFACE%NODES - ELSE - LOCAL_ERROR="Mesh number "//TRIM(NUMBER_TO_VSTRING(MESH%USER_NUMBER,"*",ERR,ERROR))// & - & " does not have an associated region or interface." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - ENDIF - ENDIF - IF(ASSOCIATED(NODES)) THEN - IF(ASSOCIATED(MESHNODES%NODES)) THEN - LOCAL_ERROR="Mesh number "//TRIM(NUMBER_TO_VSTRING(MESH%USER_NUMBER,"*",ERR,ERROR))// & - & " already has associated mesh topology nodes." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*998) - ELSE - !Work out what nodes are in the mesh - CALL TREE_CREATE_START(MESH_NODES_TREE,ERR,ERROR,*999) - CALL TREE_INSERT_TYPE_SET(MESH_NODES_TREE,TREE_NO_DUPLICATES_ALLOWED,ERR,ERROR,*999) - CALL TREE_CREATE_FINISH(MESH_NODES_TREE,ERR,ERROR,*999) - DO ne=1,ELEMENTS%NUMBER_OF_ELEMENTS - BASIS=>ELEMENTS%ELEMENTS(ne)%BASIS - DO nn=1,BASIS%NUMBER_OF_NODES - np=ELEMENTS%ELEMENTS(ne)%GLOBAL_ELEMENT_NODES(nn) - CALL TREE_ITEM_INSERT(MESH_NODES_TREE,np,np,INSERT_STATUS,ERR,ERROR,*999) - ENDDO !nn - ENDDO !ne - CALL TREE_DETACH_AND_DESTROY(MESH_NODES_TREE,NUMBER_OF_MESH_NODES,MESH_NODES,ERR,ERROR,*999) - !Set up the mesh nodes. - ALLOCATE(MESHNODES%NODES(NUMBER_OF_MESH_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate mesh topology nodes nodes.",ERR,ERROR,*999) - CALL TREE_CREATE_START(MESHNODES%NODES_TREE,ERR,ERROR,*999) - CALL TREE_INSERT_TYPE_SET(MESHNODES%NODES_TREE,TREE_NO_DUPLICATES_ALLOWED,ERR,ERROR,*999) - CALL TREE_CREATE_FINISH(MESHNODES%NODES_TREE,ERR,ERROR,*999) - DO np=1,NUMBER_OF_MESH_NODES - CALL MESH_TOPOLOGY_NODE_INITIALISE(MESHNODES%NODES(np),ERR,ERROR,*999) - MESHNODES%NODES(np)%MESH_NUMBER=np - MESHNODES%NODES(np)%GLOBAL_NUMBER=MESH_NODES(np) - MESHNODES%NODES(np)%USER_NUMBER=NODES%NODES(MESH_NODES(np))%USER_NUMBER - CALL TREE_ITEM_INSERT(MESHNODES%NODES_TREE,MESH_NODES(np),np,INSERT_STATUS,ERR,ERROR,*999) - ENDDO !np - MESHNODES%NUMBER_OF_NODES=NUMBER_OF_MESH_NODES - IF(ASSOCIATED(MESH_NODES)) DEALLOCATE(MESH_NODES) - !Now recalculate the mesh element nodes - DO ne=1,ELEMENTS%NUMBER_OF_ELEMENTS - BASIS=>ELEMENTS%ELEMENTS(ne)%BASIS - ALLOCATE(ELEMENTS%ELEMENTS(ne)%MESH_ELEMENT_NODES(BASIS%NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate mesh topology elements mesh element nodes.",ERR,ERROR,*999) - DO nn=1,BASIS%NUMBER_OF_NODES - np=ELEMENTS%ELEMENTS(ne)%GLOBAL_ELEMENT_NODES(nn) - NULLIFY(TREE_NODE) - CALL TREE_SEARCH(MESHNODES%NODES_TREE,np,TREE_NODE,ERR,ERROR,*999) - IF(ASSOCIATED(TREE_NODE)) THEN - CALL TREE_NODE_VALUE_GET(MESHNODES%NODES_TREE,TREE_NODE,MESH_NUMBER,ERR,ERROR,*999) - ELEMENTS%ELEMENTS(ne)%MESH_ELEMENT_NODES(nn)=MESH_NUMBER - ELSE - LOCAL_ERROR="Could not find global node "//TRIM(NUMBER_TO_VSTRING(np,"*",ERR,ERROR))//" (user node "// & - & TRIM(NUMBER_TO_VSTRING(NODES%NODES(np)%USER_NUMBER,"*",ERR,ERROR))//") in the mesh nodes." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) - ENDIF - ENDDO !nn - ENDDO !ne - ENDIF + INTEGER(INTG) :: dummyErr,elementIdx,insertStatus,localNodeIdx,globalNode,meshNodeIdx,meshNode,numberOfNodes + INTEGER(INTG), POINTER :: globalNodeNumbers(:) + TYPE(BASIS_TYPE), POINTER :: basis + TYPE(MESH_TYPE), POINTER :: mesh + TYPE(MeshElementsType), POINTER :: elements + TYPE(MeshNodesType), POINTER :: meshNodes + TYPE(NODES_TYPE), POINTER :: nodes + TYPE(TREE_TYPE), POINTER :: globalNodesTree + TYPE(TREE_NODE_TYPE), POINTER :: treeNode + TYPE(VARYING_STRING) :: dummyError,localError + + NULLIFY(globalNodeNumbers) + NULLIFY(globalNodesTree) + + CALL ENTERS("MeshTopologyNodesCalculate",err,error,*998) + + IF(ASSOCIATED(topology)) THEN + elements=>topology%elements + IF(ASSOCIATED(elements)) THEN + meshNodes=>topology%nodes + IF(ASSOCIATED(meshNodes)) THEN + mesh=>topology%mesh + IF(ASSOCIATED(mesh)) THEN + NULLIFY(nodes) + CALL MeshGlobalNodesGet(mesh,nodes,err,error,*999) + IF(ALLOCATED(meshNodes%nodes)) THEN + localError="Mesh number "//TRIM(NumberToVString(mesh%USER_NUMBER,"*",err,error))// & + & " already has allocated mesh topology nodes." + CALL FlagError(localError,err,error,*998) ELSE - IF(ASSOCIATED(REGION)) THEN - LOCAL_ERROR="Mesh number "//TRIM(NUMBER_TO_VSTRING(MESH%USER_NUMBER,"*",ERR,ERROR))// & - & " does not have any nodes associated with the mesh region." - ELSE - LOCAL_ERROR="Mesh number "//TRIM(NUMBER_TO_VSTRING(MESH%USER_NUMBER,"*",ERR,ERROR))// & - & " does not have any nodes associated with the mesh interface." - ENDIF - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*998) + !Work out what nodes are in the mesh + CALL TREE_CREATE_START(globalNodesTree,err,error,*999) + CALL TREE_INSERT_TYPE_SET(globalNodesTree,TREE_NO_DUPLICATES_ALLOWED,err,error,*999) + CALL TREE_CREATE_FINISH(globalNodesTree,err,error,*999) + DO elementIdx=1,elements%NUMBER_OF_ELEMENTS + basis=>elements%elements(elementIdx)%basis + DO localNodeIdx=1,basis%NUMBER_OF_NODES + globalNode=elements%elements(elementIdx)%GLOBAL_ELEMENT_NODES(localNodeIdx) + CALL TREE_ITEM_INSERT(globalNodesTree,globalNode,globalNode,insertStatus,err,error,*999) + ENDDO !localNodeIdx + ENDDO !elementIdx + CALL TREE_DETACH_AND_DESTROY(globalNodesTree,numberOfNodes,globalNodeNumbers,err,error,*999) + !Set up the mesh nodes. + ALLOCATE(meshNodes%nodes(numberOfNodes),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate mesh topology nodes nodes.",err,error,*999) + CALL TREE_CREATE_START(meshNodes%nodesTree,err,error,*999) + CALL TREE_INSERT_TYPE_SET(meshNodes%nodesTree,TREE_NO_DUPLICATES_ALLOWED,err,error,*999) + CALL TREE_CREATE_FINISH(meshNodes%nodesTree,err,error,*999) + DO meshNodeIdx=1,numberOfNodes + CALL MeshTopologyNodeInitialise(meshNodes%nodes(meshNodeIdx),err,error,*999) + meshNodes%nodes(meshNodeIdx)%meshNumber=meshNodeIdx + meshNodes%nodes(meshNodeIdx)%globalNumber=globalNodeNumbers(meshNodeIdx) + meshNodes%nodes(meshNodeIdx)%userNumber=nodes%nodes(globalNodeNumbers(meshNodeIdx))%USER_NUMBER + CALL TREE_ITEM_INSERT(meshNodes%nodesTree,globalNodeNumbers(meshNodeIdx),meshNodeIdx,insertStatus,err,error,*999) + ENDDO !nodeIdx + meshNodes%numberOfNodes=numberOfNodes + IF(ASSOCIATED(globalNodeNumbers)) DEALLOCATE(globalNodeNumbers) + !Now recalculate the mesh element nodes + DO elementIdx=1,elements%NUMBER_OF_ELEMENTS + basis=>elements%elements(elementIdx)%basis + ALLOCATE(elements%elements(elementIdx)%MESH_ELEMENT_NODES(basis%NUMBER_OF_NODES),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate mesh topology elements mesh element nodes.",err,error,*999) + DO localNodeIdx=1,basis%NUMBER_OF_NODES + globalNode=elements%elements(elementIdx)%GLOBAL_ELEMENT_NODES(localNodeIdx) + NULLIFY(treeNode) + CALL TREE_SEARCH(meshNodes%nodesTree,globalNode,treeNode,err,error,*999) + IF(ASSOCIATED(treeNode)) THEN + CALL TREE_NODE_VALUE_GET(meshNodes%nodesTree,treeNode,meshNode,err,error,*999) + elements%elements(elementIdx)%MESH_ELEMENT_NODES(localNodeIdx)=meshNode + ELSE + localError="Could not find global node "//TRIM(NumberToVString(globalNode,"*",err,error))//" (user node "// & + & TRIM(NumberToVString(nodes%nodes(globalNode)%USER_NUMBER,"*",err,error))//") in the mesh nodes." + CALL FlagError(localError,err,error,*999) + ENDIF + ENDDO !localNodeIdx + ENDDO !elementIdx ENDIF ELSE - CALL FLAG_ERROR("Mesh topology mesh is not associated.",ERR,ERROR,*998) + CALL FlagError("Mesh topology mesh is not associated.",err,error,*998) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology nodes is not associated.",ERR,ERROR,*998) + CALL FlagError("Mesh topology nodes is not associated.",err,error,*998) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology elements is not associated.",ERR,ERROR,*998) + CALL FlagError("Mesh topology elements is not associated.",err,error,*998) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology is not associated.",ERR,ERROR,*998) + CALL FlagError("Mesh topology is not associated.",err,error,*998) ENDIF IF(DIAGNOSTICS1) THEN - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,"Number of mesh global nodes = ",MESHNODES%NUMBER_OF_NODES,ERR,ERROR,*999) -!Using "NODES%NUMBER_OF_NODES" seems to be a bug, since it exceeds the size of the array "MESHNODES%NODES(np)" -! DO np=1,NODES%NUMBER_OF_NODES - DO np=1,MESHNODES%NUMBER_OF_NODES - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Mesh global node number = ",np,ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Global node number = ",MESHNODES%NODES(np)%GLOBAL_NUMBER, & - & ERR,ERROR,*999) - ENDDO !np + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"Number of mesh nodes = ",meshNodes%numberOfNodes,err,error,*999) + DO meshNodeIdx=1,meshNodes%numberOfNodes + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Mesh node number = ",meshNodeIdx,err,error,*999) + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Global node number = ",meshNodes%nodes(meshNodeIdx)%globalNumber, & + & err,error,*999) + ENDDO !meshNodeIdx ENDIF - CALL EXITS("MESH_TOPOLOGY_NODES_CALCULATE") + CALL Exits("MeshTopologyNodesCalculate") RETURN -999 IF(ASSOCIATED(MESH_NODES)) DEALLOCATE(MESH_NODES) - IF(ASSOCIATED(MESH_NODES_TREE)) CALL TREE_DESTROY(MESH_NODES_TREE,DUMMY_ERR,DUMMY_ERROR,*998) -998 CALL ERRORS("MESH_TOPOLOGY_NODES_CALCULATE",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_NODES_CALCULATE") +999 IF(ASSOCIATED(globalNodeNumbers)) DEALLOCATE(globalNodeNumbers) + IF(ASSOCIATED(globalNodesTree)) CALL TREE_DESTROY(globalNodesTree,dummyErr,dummyError,*998) +998 CALL Errors("MeshTopologyNodesCalculate",err,error) + CALL Exits("MeshTopologyNodesCalculate") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_NODES_CALCULATE + END SUBROUTINE MeshTopologyNodesCalculate + + ! + !================================================================================================================================ + ! + + !>Destroys the nodes in a mesh topology. + SUBROUTINE MeshTopologyNodesDestroy(nodes,err,error,*) + + !Argument variables + TYPE(MeshNodesType), POINTER :: nodes !Returns a pointer to the mesh nodes for a given mesh component. + SUBROUTINE MeshTopologyNodesGet(mesh,meshComponentNumber,nodes,err,error,*) + + !Argument variables + TYPE(MESH_TYPE), POINTER :: mesh !0.AND.meshComponentNumber<=mesh%NUMBER_OF_COMPONENTS) THEN + IF(ASSOCIATED(nodes)) THEN + CALL FlagError("Nodes is already associated.",err,error,*998) + ELSE + IF(ASSOCIATED(mesh%topology(meshComponentNumber)%ptr)) THEN + IF(ASSOCIATED(mesh%topology(meshComponentNumber)%ptr%nodes)) THEN + nodes=>mesh%topology(meshComponentNumber)%ptr%nodes + ELSE + CALL FlagError("Mesh topology nodes is not associated.",err,error,*999) + ENDIF + ELSE + CALL FlagError("Mesh topology is not associated.",err,error,*999) + ENDIF + ENDIF + ELSE + localError="The specified mesh component number of "//TRIM(NumberToVString(meshComponentNumber,"*",err,error))// & + & " is invalid. The component number must be between 1 and "// & + & TRIM(NumberToVString(mesh%NUMBER_OF_COMPONENTS,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + ELSE + CALL FlagError("Mesh is not associated",err,error,*998) + ENDIF + + CALL Exits("MeshTopologyNodesGet") + RETURN +999 NULLIFY(nodes) +998 CALL Errors("MeshTopologyNodesGet",err,error) + CALL Exits("MeshTopologyNodesGet") + RETURN 1 + + END SUBROUTINE MeshTopologyNodesGet ! !================================================================================================================================ ! !>Finalises the given mesh topology node. - SUBROUTINE MESH_TOPOLOGY_NODE_DERIVATIVE_FINALISE(NODE_DERIVATIVE,ERR,ERROR,*) + SUBROUTINE MeshTopologyNodeDerivativeFinalise(nodeDerivative,err,error,*) !Argument variables - TYPE(MESH_NODE_DERIVATIVE_TYPE) :: NODE_DERIVATIVE !Initialises the given mesh topology node. - SUBROUTINE MESH_TOPOLOGY_NODE_DERIVATIVE_INITIALISE(NODE_DERIVATIVE,ERR,ERROR,*) + SUBROUTINE MeshTopologyNodeDerivativeInitialise(nodeDerivative,err,error,*) !Argument variables - TYPE(MESH_NODE_DERIVATIVE_TYPE) :: NODE_DERIVATIVE !Calculates the number of derivatives at each node in a topology. - SUBROUTINE MESH_TOPOLOGY_NODES_DERIVATIVES_CALCULATE(TOPOLOGY,ERR,ERROR,*) + SUBROUTINE MeshTopologyNodesDerivativesCalculate(topology,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !TOPOLOGY%ELEMENTS - IF(ASSOCIATED(TOPOLOGY%NODES)) THEN - NODES=>TOPOLOGY%NODES + IF(ASSOCIATED(topology)) THEN + elements=>topology%elements + IF(ASSOCIATED(elements)) THEN + nodes=>topology%nodes + IF(ASSOCIATED(nodes)) THEN !Loop over the mesh nodes - DO node_idx=1,NODES%NUMBER_OF_NODES - !Calculate the number of derivatives and versions at each node. This needs to be calculated by looking at the mesh elements - !as we may have an adjacent element in another domain with a higher order basis also with versions. - NULLIFY(NODE_DERIVATIVE_LIST) - CALL LIST_CREATE_START(NODE_DERIVATIVE_LIST,ERR,ERROR,*999) - CALL LIST_DATA_TYPE_SET(NODE_DERIVATIVE_LIST,LIST_INTG_TYPE,ERR,ERROR,*999) - CALL LIST_INITIAL_SIZE_SET(NODE_DERIVATIVE_LIST,8,ERR,ERROR,*999) - CALL LIST_CREATE_FINISH(NODE_DERIVATIVE_LIST,ERR,ERROR,*999) - MAX_NUMBER_OF_DERIVATIVES=-1 - DO elem_idx=1,NODES%NODES(node_idx)%NUMBER_OF_SURROUNDING_ELEMENTS - ne=NODES%NODES(node_idx)%SURROUNDING_ELEMENTS(elem_idx) - BASIS=>ELEMENTS%ELEMENTS(ne)%BASIS + DO nodeIdx=1,nodes%numberOfNodes + !Calculate the number of derivatives and versions at each node. This needs to be calculated by looking at the + !mesh elements as we may have an adjacent element in another domain with a higher order basis also with versions. + NULLIFY(nodeDerivativeList) + CALL LIST_CREATE_START(nodeDerivativeList,err,error,*999) + CALL LIST_DATA_TYPE_SET(nodeDerivativeList,LIST_INTG_TYPE,err,error,*999) + CALL LIST_INITIAL_SIZE_SET(nodeDerivativeList,8,err,error,*999) + CALL LIST_CREATE_FINISH(nodeDerivativeList,err,error,*999) + maxNumberOfDerivatives=-1 + DO elementIdx=1,nodes%nodes(nodeIdx)%numberOfSurroundingElements + element=nodes%nodes(nodeIdx)%surroundingElements(elementIdx) + basis=>elements%elements(element)%basis !Find the local node corresponding to this node - FOUND=.FALSE. - DO nn=1,BASIS%NUMBER_OF_NODES - IF(ELEMENTS%ELEMENTS(ne)%MESH_ELEMENT_NODES(nn)==node_idx) THEN - FOUND=.TRUE. + found=.FALSE. + DO localNodeIdx=1,basis%NUMBER_OF_NODES + IF(elements%elements(element)%MESH_ELEMENT_NODES(localNodeIdx)==nodeIdx) THEN + found=.TRUE. EXIT ENDIF ENDDO !nn - IF(FOUND) THEN - DO derivative_idx=1,BASIS%NUMBER_OF_DERIVATIVES(nn) - CALL LIST_ITEM_ADD(NODE_DERIVATIVE_LIST,BASIS%PARTIAL_DERIVATIVE_INDEX(derivative_idx,nn),ERR,ERROR,*999) - ENDDO !derivative_idx - IF(BASIS%NUMBER_OF_DERIVATIVES(nn)>MAX_NUMBER_OF_DERIVATIVES)MAX_NUMBER_OF_DERIVATIVES= & - & BASIS%NUMBER_OF_DERIVATIVES(nn) + IF(found) THEN + DO derivativeIdx=1,basis%NUMBER_OF_DERIVATIVES(localNodeIdx) + CALL LIST_ITEM_ADD(nodeDerivativeList,basis%PARTIAL_DERIVATIVE_INDEX(derivativeIdx,localNodeIdx),err,error,*999) + ENDDO !derivativeIdx + IF(basis%NUMBER_OF_DERIVATIVES(localNodeIdx)>maxNumberOfDerivatives) & + & maxNumberOfDerivatives=basis%NUMBER_OF_DERIVATIVES(localNodeidx) ELSE - CALL FLAG_ERROR("Could not find local node.",ERR,ERROR,*999) + CALL FlagError("Could not find local node.",err,error,*999) ENDIF ENDDO !elem_idx - CALL LIST_REMOVE_DUPLICATES(NODE_DERIVATIVE_LIST,ERR,ERROR,*999) - CALL LIST_DETACH_AND_DESTROY(NODE_DERIVATIVE_LIST,NUMBER_OF_DERIVATIVES,DERIVATIVES,ERR,ERROR,*999) - IF(NUMBER_OF_DERIVATIVES==MAX_NUMBER_OF_DERIVATIVES) THEN + CALL LIST_REMOVE_DUPLICATES(nodeDerivativeList,err,error,*999) + CALL LIST_DETACH_AND_DESTROY(nodeDerivativeList,numberOfDerivatives,derivatives,err,error,*999) + IF(numberOfDerivatives==maxNumberOfDerivatives) THEN !Set up the node derivatives. - ALLOCATE(NODES%NODES(node_idx)%DERIVATIVES(MAX_NUMBER_OF_DERIVATIVES),STAT=ERR) - NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES=MAX_NUMBER_OF_DERIVATIVES - DO derivative_idx=1,NUMBER_OF_DERIVATIVES - CALL MESH_TOPOLOGY_NODE_DERIVATIVE_INITIALISE(NODES%NODES(node_idx)%DERIVATIVES(derivative_idx),ERR,ERROR,*999) - NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%PARTIAL_DERIVATIVE_INDEX = DERIVATIVES(derivative_idx) - global_deriv=PARTIAL_DERIVATIVE_GLOBAL_DERIVATIVE_MAP(DERIVATIVES(derivative_idx)) - IF(global_deriv/=0) THEN - NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%GLOBAL_DERIVATIVE_INDEX=global_deriv + ALLOCATE(nodes%nodes(nodeIdx)%derivatives(maxNumberOfDerivatives),STAT=err) + nodes%nodes(nodeIdx)%numberOfDerivatives=maxNumberOfDerivatives + DO derivativeIdx=1,numberOfDerivatives + CALL MeshTopologyNodeDerivativeInitialise(nodes%nodes(nodeIdx)%derivatives(derivativeIdx),err,error,*999) + nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%partialDerivativeIndex = derivatives(derivativeIdx) + globalDerivative=PARTIAL_DERIVATIVE_GLOBAL_DERIVATIVE_MAP(derivatives(derivativeIdx)) + IF(globalDerivative/=0) THEN + nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%globalDerivativeIndex=globalDerivative ELSE - LOCAL_ERROR="The partial derivative index of "//TRIM(NUMBER_TO_VSTRING(DERIVATIVES(derivative_idx),"*", & - & ERR,ERROR))//" for derivative number "//TRIM(NUMBER_TO_VSTRING(derivative_idx,"*",ERR,ERROR))// & + localError="The partial derivative index of "//TRIM(NumberToVstring(derivatives(derivativeIdx),"*", & + & err,error))//" for derivative number "//TRIM(NumberToVstring(derivativeIdx,"*",err,error))// & & " does not have a corresponding global derivative." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + CALL FlagError(localError,err,error,*999) ENDIF - ENDDO !derivative_idx - DEALLOCATE(DERIVATIVES) + ENDDO !derivativeIdx + DEALLOCATE(derivatives) ELSE - LOCAL_ERROR="Invalid mesh configuration. User node "// & - & TRIM(NUMBER_TO_VSTRING(NODES%NODES(node_idx)%USER_NUMBER,"*",ERR,ERROR))// & + localError="Invalid mesh configuration. User node "// & + & TRIM(NumberToVstring(nodes%nodes(nodeIdx)%userNumber,"*",err,error))// & & " has inconsistent derivative directions." - CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) + CALL FlagError(localError,err,error,*999) ENDIF - ENDDO !node_idx + ENDDO !nodeIdx ELSE - CALL FLAG_ERROR("Mesh topology nodes is not associated.",ERR,ERROR,*999) + CALL FlagError("Mesh topology nodes is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology elements is not associated.",ERR,ERROR,*999) + CALL FlagError("Mesh topology elements is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology is not associated.",ERR,ERROR,*999) + CALL FlagError("Mesh topology is not associated.",err,error,*999) ENDIF - IF(DIAGNOSTICS1) THEN - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,"Number of mesh global nodes = ",NODES%NUMBER_OF_NODES,ERR,ERROR,*999) - DO node_idx=1,NODES%NUMBER_OF_NODES - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Mesh global node number = ",node_idx,ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Number of derivatives = ",NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES, & - & ERR,ERROR,*999) - DO derivative_idx=1,NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - !TODO: change output string below so that it writes out derivative_idx index as well - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Global derivative index(derivative_idx) = ", & - & NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%GLOBAL_DERIVATIVE_INDEX,ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Partial derivative index(derivative_idx) = ", & - & NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%PARTIAL_DERIVATIVE_INDEX,ERR,ERROR,*999) - ENDDO !derivative_idx + IF(diagnostics1) THEN + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"Number of mesh global nodes = ",nodes%numberOfNodes,err,error,*999) + DO nodeIdx=1,nodes%numberOfNodes + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Mesh global node number = ",nodeIdx,err,error,*999) + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Number of derivatives = ",nodes%nodes(nodeIdx)%numberOfDerivatives, & + & err,error,*999) + DO derivativeIdx=1,nodes%nodes(nodeIdx)%numberOfDerivatives + !TODO: change output string below so that it writes out derivativeIdx index as well + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Global derivative index(derivativeIdx) = ", & + & nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%globalDerivativeIndex,err,error,*999) + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Partial derivative index(derivativeIdx) = ", & + & nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%partialDerivativeIndex,err,error,*999) + ENDDO !derivativeIdx ENDDO !node_idx ENDIF - CALL EXITS("MESH_TOPOLOGY_NODES_DERIVATIVES_CALCULATE") + CALL EXITS("MeshTopologyNodesDerivativesCalculate") RETURN -999 IF(ALLOCATED(DERIVATIVES)) DEALLOCATE(DERIVATIVES) - IF(ASSOCIATED(NODE_DERIVATIVE_LIST)) CALL LIST_DESTROY(NODE_DERIVATIVE_LIST,ERR,ERROR,*998) -998 CALL ERRORS("MESH_TOPOLOGY_NODES_DERIVATIVES_CALCULATE",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_NODES_DERIVATIVES_CALCULATE") +999 IF(ALLOCATED(derivatives)) DEALLOCATE(derivatives) + IF(ASSOCIATED(nodeDerivativeList)) CALL LIST_DESTROY(nodeDerivativeList,err,error,*998) +998 CALL errorS("MeshTopologyNodesDerivativesCalculate",err,error) + CALL EXITS("MeshTopologyNodesDerivativesCalculate") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_NODES_DERIVATIVES_CALCULATE + END SUBROUTINE MeshTopologyNodesDerivativesCalculate + + ! + !================================================================================================================================ + ! + + !>Returns the number of derivatives for a node in a mesh + SUBROUTINE MeshTopologyNodeNumberOfDerivativesGet(meshNodes,userNumber,numberOfDerivatives,err,error,*) + + !Argument variables + TYPE(MeshNodesType), POINTER :: meshNodes !meshNodes%meshComponentTopology + IF(ASSOCIATED(meshComponentTopology)) THEN + mesh=>meshComponentTopology%mesh + IF(ASSOCIATED(mesh)) THEN + meshComponentNumber=meshComponentTopology%meshComponentNumber + localError="The user node number "//TRIM(NumberToVString(userNumber,"*",err,error))// & + & " does not exist in mesh component number "//TRIM(NumberToVString(meshComponentNumber,"*",err,error))// & + & " of mesh number "//TRIM(NumberToVString(mesh%USER_NUMBER,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ELSE + CALL FlagError("Mesh component topology mesh is not associated.",err,error,*999) + ENDIF + ELSE + CALL FlagError("Mesh nodes mesh component topology is not associated.",err,error,*999) + ENDIF + ENDIF + ELSE + CALL FlagError("Mesh nodes is not associated.",err,error,*999) + ENDIF + + CALL Exits("MeshTopologyNodeNumberOfDerivativesGet") + RETURN +999 CALL Errors("MeshTopologyNodeNumberOfDerivativesGet",err,error) + CALL Exits("MeshTopologyNodeNumberOfDerivativesGet") + RETURN 1 + + END SUBROUTINE MeshTopologyNodeNumberOfDerivativesGet + + ! + !================================================================================================================================ + ! + + !>Returns the global derivative numbers for a node in mesh nodes + SUBROUTINE MeshTopologyNodeDerivativesGet(meshNodes,userNumber,derivatives,err,error,*) + + !Argument variables + TYPE(MeshNodesType), POINTER :: meshNodes !=numberOfDerivatives) THEN + DO derivativeIdx=1,numberOfDerivatives + derivatives(derivativeIdx)=meshNodes%nodes(meshNumber)%derivatives(derivativeIdx)%globalDerivativeIndex + ENDDO !derivativeIdx + ELSE + localError="The size of the supplied derivatives array of "// & + & TRIM(NumberToVString(SIZE(derivatives,1),"*",err,error))// & + & " is too small. The size should be >= "// & + & TRIM(NumberToVString(numberOfDerivatives,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + ELSE + meshComponentTopology=>meshNodes%meshComponentTopology + IF(ASSOCIATED(meshComponentTopology)) THEN + mesh=>meshComponentTopology%mesh + IF(ASSOCIATED(mesh)) THEN + meshComponentNumber=meshComponentTopology%meshComponentNumber + localError="The user node number "//TRIM(NumberToVString(userNumber,"*",err,error))// & + & " does not exist in mesh component number "//TRIM(NumberToVString(meshComponentNumber,"*",err,error))// & + & " of mesh number "//TRIM(NumberToVString(mesh%USER_NUMBER,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ELSE + CALL FlagError("Mesh component topology mesh is not associated.",err,error,*999) + ENDIF + ELSE + CALL FlagError("Mesh nodes mesh component topology is not associated.",err,error,*999) + ENDIF + ENDIF + ELSE + CALL FlagError("Mesh nodes is not associated.",err,error,*999) + ENDIF + + CALL Exits("MeshTopologyNodeDerivativesGet") + RETURN +999 CALL Errors("MeshTopologyNodeDerivativesGet",err,error) + CALL Exits("MeshTopologyNodeDerivativesGet") + RETURN 1 + + END SUBROUTINE MeshTopologyNodeDerivativesGet + + ! + !================================================================================================================================ + ! + + !>Returns the number of versions for a derivative of a node in mesh nodes + SUBROUTINE MeshTopologyNodeNumberOfVersionsGet(meshNodes,derivativeNumber,userNumber,numberOfVersions,err,error,*) + + !Argument variables + TYPE(MeshNodesType), POINTER :: meshNodes !=1.AND.derivativeNumber<=meshNodes%nodes(meshNumber)%numberOfDerivatives) THEN + numberOfVersions=meshNodes%nodes(meshNumber)%derivatives(derivativeNumber)%numberOfVersions + ELSE + localError="The specified derivative index of "// & + & TRIM(NumberToVString(derivativeNumber,"*",err,error))// & + & " is invalid. The derivative index must be >= 1 and <= "// & + & TRIM(NumberToVString(meshNodes%nodes(meshNumber)%numberOfDerivatives,"*",err,error))// & + & " for user node number "//TRIM(NumberToVString(userNumber,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + ELSE + meshComponentTopology=>meshNodes%meshComponentTopology + IF(ASSOCIATED(meshComponentTopology)) THEN + mesh=>meshComponentTopology%mesh + IF(ASSOCIATED(mesh)) THEN + meshComponentNumber=meshComponentTopology%meshComponentNumber + localError="The user node number "//TRIM(NumberToVString(userNumber,"*",err,error))// & + & " does not exist in mesh component number "//TRIM(NumberToVString(meshComponentNumber,"*",err,error))// & + & " of mesh number "//TRIM(NumberToVString(mesh%USER_NUMBER,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ELSE + CALL FlagError("Mesh component topology mesh is not associated.",err,error,*999) + ENDIF + ENDIF + ENDIF + ELSE + CALL FlagError("Mesh nodes is not associated.",err,error,*999) + ENDIF + + CALL Exits("MeshTopologyNodeNumberOfVersionsGet") + RETURN +999 CALL Errors("MeshTopologyNodeNumberOfVersionsGet",err,error) + CALL Exits("MeshTopologyNodeNumberOfVersionsGet") + RETURN 1 + + END SUBROUTINE MeshTopologyNodeNumberOfVersionsGet + ! !================================================================================================================================ ! !>Calculates the number of versions at each node in a topology. - SUBROUTINE MESH_TOPOLOGY_NODES_VERSIONS_CALCULATE(TOPOLOGY,ERR,ERROR,*) + SUBROUTINE MeshTopologyNodesVersionCalculate(topology,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !TOPOLOGY%ELEMENTS - IF(ASSOCIATED(TOPOLOGY%NODES)) THEN - NODES=>TOPOLOGY%NODES + IF(ASSOCIATED(topology)) THEN + elements=>topology%elements + IF(ASSOCIATED(elements)) THEN + nodes=>topology%nodes + IF(ASSOCIATED(nodes)) THEN !Loop over the mesh elements !Calculate the number of versions at each node. This needs to be calculated by looking at all the mesh elements - !as we may have an adjacent elements in another domain with a higher order basis along with different versions being assigned to its derivatives. + !as we may have an adjacent elements in another domain with a higher order basis along with different versions + !being assigned to its derivatives. !\todo : See if there are any constraints that can be applied to restrict the amount of memory being allocated here - ALLOCATE(NODE_VERSION_LIST(MAXIMUM_GLOBAL_DERIV_NUMBER,NODES%NUMBER_OF_NODES),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node version list.",ERR,ERROR,*999) - DO node_idx=1,NODES%NUMBER_OF_NODES - DO derivative_idx=1,NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - NULLIFY(NODE_VERSION_LIST(derivative_idx,node_idx)%PTR) - CALL LIST_CREATE_START(NODE_VERSION_LIST(derivative_idx,node_idx)%PTR,ERR,ERROR,*999) - CALL LIST_DATA_TYPE_SET(NODE_VERSION_LIST(derivative_idx,node_idx)%PTR,LIST_INTG_TYPE,ERR,ERROR,*999) - CALL LIST_INITIAL_SIZE_SET(NODE_VERSION_LIST(derivative_idx,node_idx)%PTR,8,ERR,ERROR,*999) - CALL LIST_CREATE_FINISH(NODE_VERSION_LIST(derivative_idx,node_idx)%PTR,ERR,ERROR,*999) - ENDDO!derivative_idx - ENDDO!node_idx - DO ne=1,ELEMENTS%NUMBER_OF_ELEMENTS - BASIS=>ELEMENTS%ELEMENTS(ne)%BASIS - DO nn=1,BASIS%NUMBER_OF_NODES - DO derivative_idx=1,BASIS%NUMBER_OF_DERIVATIVES(nn) - CALL LIST_ITEM_ADD(NODE_VERSION_LIST(derivative_idx,ELEMENTS%ELEMENTS(ne)%MESH_ELEMENT_NODES(nn))%PTR, & - & ELEMENTS%ELEMENTS(ne)%USER_ELEMENT_NODE_VERSIONS(derivative_idx,nn),ERR,ERROR,*999) - ENDDO!derivative_idx - ENDDO!nn - ENDDO!ne - DO node_idx=1,NODES%NUMBER_OF_NODES - DO derivative_idx=1,NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - CALL LIST_REMOVE_DUPLICATES(NODE_VERSION_LIST(derivative_idx,node_idx)%PTR,ERR,ERROR,*999) - CALL LIST_DETACH_AND_DESTROY(NODE_VERSION_LIST(derivative_idx,node_idx)%PTR,NUMBER_OF_VERSIONS,VERSIONS, & - & ERR,ERROR,*999) - NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS = MAXVAL(VERSIONS(1:NUMBER_OF_VERSIONS)) - ALLOCATE(NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%USER_VERSION_NUMBERS(NODES%NODES(node_idx)% & - & DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate node global derivative index.",ERR,ERROR,*999) - DO version_idx=1,NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS - NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%USER_VERSION_NUMBERS(version_idx) = version_idx - ENDDO !version_idx - DEALLOCATE(VERSIONS) - ENDDO!derivative_idx - ENDDO!node_idx - DEALLOCATE(NODE_VERSION_LIST) - NULLIFY(NODE_VERSION_LIST) + ALLOCATE(nodeVersionList(MAXIMUM_GLOBAL_DERIV_NUMBER,nodes%numberOfNodes),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate node version list.",err,error,*999) + DO nodeIdx=1,nodes%numberOfNodes + DO derivativeIdx=1,nodes%nodes(nodeIdx)%numberOfDerivatives + NULLIFY(nodeVersionList(derivativeIdx,nodeIdx)%ptr) + CALL LIST_CREATE_START(nodeVersionList(derivativeIdx,nodeIdx)%ptr,err,error,*999) + CALL LIST_DATA_TYPE_SET(nodeVersionList(derivativeIdx,nodeIdx)%ptr,LIST_INTG_TYPE,err,error,*999) + CALL LIST_INITIAL_SIZE_SET(nodeVersionList(derivativeIdx,nodeIdx)%ptr,8,err,error,*999) + CALL LIST_CREATE_FINISH(nodeVersionList(derivativeIdx,nodeIdx)%ptr,err,error,*999) + ENDDO!derivativeIdx + ENDDO!nodeIdx + DO element=1,elements%NUMBER_OF_ELEMENTS + basis=>elements%elements(element)%basis + DO localNodeIdx=1,BASIS%NUMBER_OF_NODES + DO derivativeIdx=1,basis%NUMBER_OF_DERIVATIVES(localNodeIdx) + CALL LIST_ITEM_ADD(nodeVersionList(derivativeIdx,elements%elements(element)% & + & MESH_ELEMENT_NODES(localNodeIdx))%ptr,elements%elements(element)%USER_ELEMENT_NODE_VERSIONS( & + & derivativeIdx,localNodeIdx),err,error,*999) + ENDDO!derivativeIdx + ENDDO!localNodeIdx + ENDDO!element + DO nodeIdx=1,nodes%numberOfNodes + DO derivativeIdx=1,nodes%nodes(nodeIdx)%numberOfDerivatives + CALL LIST_REMOVE_DUPLICATES(nodeVersionList(derivativeIdx,nodeIdx)%ptr,err,error,*999) + CALL LIST_DETACH_AND_DESTROY(nodeVersionList(derivativeIdx,nodeIdx)%ptr,numberOfVersions,versions, & + & err,error,*999) + nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%numberOfVersions = MAXVAL(versions(1:numberOfVersions)) + ALLOCATE(nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%userVersionNumbers(nodes%nodes(nodeIdx)% & + & derivatives(derivativeIdx)%numberOfVersions),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate node global derivative index.",err,error,*999) + DO versionIdx=1,nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%numberOfVersions + nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%userVersionNumbers(versionIdx) = versionIdx + ENDDO !versionIdx + DEALLOCATE(versions) + ENDDO!derivativeIdx + ENDDO!nodeIdx + DEALLOCATE(nodeVersionList) + NULLIFY(nodeVersionList) ELSE - CALL FLAG_ERROR("Mesh topology nodes is not associated.",ERR,ERROR,*999) + CALL FlagError("Mesh topology nodes is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology elements is not associated.",ERR,ERROR,*999) + CALL FlagError("Mesh topology elements is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology is not associated.",ERR,ERROR,*999) + CALL FlagError("Mesh topology is not associated.",err,error,*999) ENDIF - IF(DIAGNOSTICS1) THEN - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE,"Number of mesh global nodes = ",NODES%NUMBER_OF_NODES,ERR,ERROR,*999) - DO node_idx=1,NODES%NUMBER_OF_NODES - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Mesh global node number = ",node_idx,ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Number of derivatives = ", & - & NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES,ERR,ERROR,*999) - DO derivative_idx=1,NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - !\todo : change output string below so that it writes out derivative_idx index as well - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Global derivative index(derivative_idx) = ", & - & NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%GLOBAL_DERIVATIVE_INDEX,ERR,ERROR,*999) - CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," Partial derivative index(derivative_idx) = ", & - & NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%PARTIAL_DERIVATIVE_INDEX,ERR,ERROR,*999) + IF(diagnostics1) THEN + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE,"Number of mesh global nodes = ",nodes%numberOfNodes,err,error,*999) + DO nodeIdx=1,nodes%numberOfNodes + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Mesh global node number = ",nodeIdx,err,error,*999) + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Number of derivatives = ", & + & nodes%nodes(nodeIdx)%numberOfDerivatives,err,error,*999) + DO derivativeIdx=1,nodes%nodes(nodeIdx)%numberOfDerivatives + !\todo : change output string below so that it writes out derivativeIdx index as well + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Global derivative index(derivativeIdx) = ", & + & nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%globalDerivativeIndex,err,error,*999) + CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Partial derivative index(derivativeIdx) = ", & + & nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%partialDerivativeIndex,err,error,*999) CALL WRITE_STRING_VECTOR(DIAGNOSTIC_OUTPUT_TYPE,1,1, & - & NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%NUMBER_OF_VERSIONS,8,8, & - & NODES%NODES(node_idx)%DERIVATIVES(derivative_idx)%USER_VERSION_NUMBERS, & - & '(" User Version index(derivative_idx,:) :",8(X,I2))','(36X,8(X,I2))',ERR,ERROR,*999) - ENDDO!derivative_idx - ENDDO !node_idx + & nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%numberOfVersions,8,8, & + & nodes%nodes(nodeIdx)%derivatives(derivativeIdx)%userVersionNumbers, & + & '(" User Version index(derivativeIdx,:) :",8(X,I2))','(36X,8(X,I2))',err,error,*999) + ENDDO!derivativeIdx + ENDDO !nodeIdx ENDIF - CALL EXITS("MESH_TOPOLOGY_NODES_VERSIONS_CALCULATE") + CALL EXITS("MeshTopologyNodesVersionCalculate") RETURN -999 IF(ALLOCATED(VERSIONS)) DEALLOCATE(VERSIONS) - IF(ASSOCIATED(NODE_VERSION_LIST)) THEN - DO node_idx=1,NODES%NUMBER_OF_NODES - DO derivative_idx=1,NODES%NODES(node_idx)%NUMBER_OF_DERIVATIVES - CALL LIST_DESTROY(NODE_VERSION_LIST(derivative_idx,node_idx)%PTR,ERR,ERROR,*998) - ENDDO !derivative_idx - ENDDO !node_idx - DEALLOCATE(NODE_VERSION_LIST) - NULLIFY(NODE_VERSION_LIST) +999 IF(ALLOCATED(versions)) DEALLOCATE(versions) + IF(ASSOCIATED(nodeVersionList)) THEN + DO nodeIdx=1,nodes%numberOfNodes + DO derivativeIdx=1,nodes%nodes(nodeIdx)%numberOfDerivatives + CALL LIST_DESTROY(nodeVersionList(derivativeIdx,nodeIdx)%ptr,err,error,*998) + ENDDO !derivativeIdx + ENDDO !nodeIdx + DEALLOCATE(nodeVersionList) ENDIF -998 CALL ERRORS("MESH_TOPOLOGY_NODES_VERSIONS_CALCULATE",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_NODES_VERSIONS_CALCULATE") +998 CALL ERRORS("MeshTopologyNodesVersionCalculate",err,error) + CALL EXITS("MeshTopologyNodesVersionCalculate") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_NODES_VERSIONS_CALCULATE + END SUBROUTINE MeshTopologyNodesVersionCalculate ! !================================================================================================================================ ! !>Calculates the element numbers surrounding a node for a mesh. - SUBROUTINE MESH_TOPOLOGY_NODES_SURROUNDING_ELEMENTS_CALCULATE(TOPOLOGY,ERR,ERROR,*) + SUBROUTINE MeshTopologySurroundingElementsCalculate(topology,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !TOPOLOGY%ELEMENTS%ELEMENTS(ne)%BASIS - DO nn=1,BASIS%NUMBER_OF_NODES - np=TOPOLOGY%ELEMENTS%ELEMENTS(ne)%MESH_ELEMENT_NODES(nn) - FOUND_ELEMENT=.FALSE. - element_no=1 - insert_position=1 - DO WHILE(element_no<=TOPOLOGY%NODES%NODES(np)%NUMBER_OF_SURROUNDING_ELEMENTS.AND..NOT.FOUND_ELEMENT) - surrounding_elem_no=TOPOLOGY%NODES%NODES(np)%SURROUNDING_ELEMENTS(element_no) - IF(surrounding_elem_no==ne) THEN - FOUND_ELEMENT=.TRUE. + INTEGER(INTG) :: element,elementIdx,insertPosition,localNodeIdx,node,surroundingElementNumber + INTEGER(INTG), POINTER :: newSurroundingElements(:) + LOGICAL :: foundElement + TYPE(BASIS_TYPE), POINTER :: basis + TYPE(MeshElementsType), POINTER :: elements + TYPE(MeshNodesType), POINTER :: nodes + + NULLIFY(newSurroundingElements) + + CALL ENTERS("MeshTopologySurroundingElementsCalculate",err,error,*999) + + IF(ASSOCIATED(topology)) THEN + elements=>topology%elements + IF(ASSOCIATED(elements)) THEN + nodes=>topology%nodes + IF(ASSOCIATED(nodes)) THEN + IF(ALLOCATED(nodes%nodes)) THEN + DO elementIdx=1,elements%NUMBER_OF_ELEMENTS + basis=>elements%elements(elementIdx)%basis + DO localNodeIdx=1,basis%NUMBER_OF_NODES + node=elements%elements(elementIdx)%MESH_ELEMENT_NODES(localNodeIdx) + foundElement=.FALSE. + element=1 + insertPosition=1 + DO WHILE(element<=nodes%nodes(node)%numberOfSurroundingElements.AND..NOT.foundElement) + surroundingElementNumber=nodes%nodes(node)%surroundingElements(element) + IF(surroundingElementNumber==elementIdx) THEN + foundElement=.TRUE. ENDIF - element_no=element_no+1 - IF(ne>=surrounding_elem_no) THEN - insert_position=element_no + element=element+1 + IF(elementIdx>=surroundingElementNumber) THEN + insertPosition=element ENDIF ENDDO - IF(.NOT.FOUND_ELEMENT) THEN + IF(.NOT.foundElement) THEN !Insert element into surrounding elements - ALLOCATE(NEW_SURROUNDING_ELEMENTS(TOPOLOGY%NODES%NODES(np)%NUMBER_OF_SURROUNDING_ELEMENTS+1),STAT=ERR) - IF(ERR/=0) CALL FLAG_ERROR("Could not allocate new surrounding elements",ERR,ERROR,*999) - IF(ASSOCIATED(TOPOLOGY%NODES%NODES(np)%SURROUNDING_ELEMENTS)) THEN - NEW_SURROUNDING_ELEMENTS(1:insert_position-1)=TOPOLOGY%NODES%NODES(np)% & - & SURROUNDING_ELEMENTS(1:insert_position-1) - NEW_SURROUNDING_ELEMENTS(insert_position)=ne - NEW_SURROUNDING_ELEMENTS(insert_position+1:TOPOLOGY%NODES%NODES(np)%NUMBER_OF_SURROUNDING_ELEMENTS+1)= & - & TOPOLOGY%NODES%NODES(np)%SURROUNDING_ELEMENTS(insert_position:TOPOLOGY%NODES%NODES(np)% & - & NUMBER_OF_SURROUNDING_ELEMENTS) - DEALLOCATE(TOPOLOGY%NODES%NODES(np)%SURROUNDING_ELEMENTS) + ALLOCATE(newSurroundingElements(nodes%nodes(node)%numberOfSurroundingElements+1),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate new surrounding elements.",err,error,*999) + IF(ASSOCIATED(nodes%nodes(node)%surroundingElements)) THEN + newSurroundingElements(1:insertPosition-1)=nodes%nodes(node)%surroundingElements(1:insertPosition-1) + newSurroundingElements(insertPosition)=elementIdx + newSurroundingElements(insertPosition+1:nodes%nodes(node)%numberOfSurroundingElements+1)= & + & nodes%nodes(node)%surroundingElements(insertPosition:nodes%nodes(node)%numberOfSurroundingElements) + DEALLOCATE(nodes%nodes(node)%surroundingElements) ELSE - NEW_SURROUNDING_ELEMENTS(1)=ne + newSurroundingElements(1)=elementIdx ENDIF - TOPOLOGY%NODES%NODES(np)%SURROUNDING_ELEMENTS=>NEW_SURROUNDING_ELEMENTS - TOPOLOGY%NODES%NODES(np)%NUMBER_OF_SURROUNDING_ELEMENTS= & - & TOPOLOGY%NODES%NODES(np)%NUMBER_OF_SURROUNDING_ELEMENTS+1 + nodes%nodes(node)%surroundingElements=>newSurroundingElements + nodes%nodes(node)%numberOfSurroundingElements=nodes%nodes(node)%numberOfSurroundingElements+1 ENDIF - ENDDO !nn - ENDDO !ne + ENDDO !localNodeIdx + ENDDO !elementIdx ELSE - CALL FLAG_ERROR("Mesh topology nodes nodes are not associated",ERR,ERROR,*999) + CALL FlagError("Mesh topology nodes nodes have not been allocated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology nodes are not associated",ERR,ERROR,*999) + CALL FlagError("Mesh topology nodes are not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology elements is not associated",ERR,ERROR,*999) + CALL FlagError("Mesh topology elements is not associated.",err,error,*999) ENDIF ELSE - CALL FLAG_ERROR("Mesh topology not associated",ERR,ERROR,*999) + CALL FlagError("Mesh topology not associated.",err,error,*999) ENDIF - CALL EXITS("MESH_TOPOLOGY_NODES_SURROUNDING_ELEMENTS_CALCULATE") + CALL EXITS("MeshTopologySurroundingElementsCalculate") RETURN -999 IF(ASSOCIATED(NEW_SURROUNDING_ELEMENTS)) DEALLOCATE(NEW_SURROUNDING_ELEMENTS) - CALL ERRORS("MESH_TOPOLOGY_NODES_SURROUNDING_ELEMENTS_CALCULATE",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_NODES_SURROUNDING_ELEMENTS_CALCULATE") +999 IF(ASSOCIATED(newSurroundingElements)) DEALLOCATE(newSurroundingElements) + CALL ERRORS("MeshTopologySurroundingElementsCalculate",err,error) + CALL EXITS("MeshTopologySurroundingElementsCalculate") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_NODES_SURROUNDING_ELEMENTS_CALCULATE + END SUBROUTINE MeshTopologySurroundingElementsCalculate ! !=============================================================================================================================== ! - !>Finalises the nodes data structures for a mesh topology and deallocates any memory. \todo pass in nodes - SUBROUTINE MESH_TOPOLOGY_NODES_FINALISE(TOPOLOGY,ERR,ERROR,*) + !>Finalises the nodes data structures for a mesh topology and deallocates any memory. + SUBROUTINE MeshTopologyNodesFinalise(nodes,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !Initialises the nodes in a given mesh topology. \todo finalise on errors - SUBROUTINE MESH_TOPOLOGY_NODES_INITIALISE(TOPOLOGY,ERR,ERROR,*) + SUBROUTINE MeshTopologyNodesInitialise(topology,err,error,*) !Argument variables - TYPE(MeshComponentTopologyType), POINTER :: TOPOLOGY !TOPOLOGY%MESH - NULLIFY(TOPOLOGY%NODES%NODES) - NULLIFY(TOPOLOGY%NODES%NODES_TREE) + ALLOCATE(topology%nodes,STAT=err) + IF(err/=0) CALL FlagError("Could not allocate topology nodes.",err,error,*999) + topology%nodes%numberOfNodes=0 + topology%nodes%meshComponentTopology=>topology + NULLIFY(topology%nodes%nodesTree) ENDIF ELSE - CALL FLAG_ERROR("Topology is not associated",ERR,ERROR,*999) + CALL FlagError("Topology is not associated.",err,error,*999) ENDIF - CALL EXITS("MESH_TOPOLOGY_NODES_INITIALISE") + CALL EXITS("MeshTopologyNodesInitialise") RETURN -999 CALL ERRORS("MESH_TOPOLOGY_NODES_INITIALISE",ERR,ERROR) - CALL EXITS("MESH_TOPOLOGY_NODES_INITIALISE") +999 CALL ERRORS("MeshTopologyNodesInitialise",err,error) + CALL EXITS("MeshTopologyNodesInitialise") RETURN 1 - END SUBROUTINE MESH_TOPOLOGY_NODES_INITIALISE + END SUBROUTINE MeshTopologyNodesInitialise ! !================================================================================================================================ @@ -9580,7 +10082,7 @@ SUBROUTINE MESHES_INITIALISE_INTERFACE(INTERFACE,ERR,ERROR,*) IF(ASSOCIATED(INTERFACE)) THEN IF(ASSOCIATED(INTERFACE%MESHES)) THEN - LOCAL_ERROR="Interface number "//TRIM(NUMBER_TO_VSTRING(INTERFACE%USER_NUMBER,"*",ERR,ERROR))// & + LOCAL_ERROR="Interface number "//TRIM(NumberToVString(INTERFACE%USER_NUMBER,"*",ERR,ERROR))// & & " already has a mesh associated" CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) ELSE @@ -9616,7 +10118,7 @@ SUBROUTINE MESHES_INITIALISE_REGION(REGION,ERR,ERROR,*) IF(ASSOCIATED(REGION)) THEN IF(ASSOCIATED(REGION%MESHES)) THEN - LOCAL_ERROR="Region number "//TRIM(NUMBER_TO_VSTRING(REGION%USER_NUMBER,"*",ERR,ERROR))// & + LOCAL_ERROR="Region number "//TRIM(NumberToVString(REGION%USER_NUMBER,"*",ERR,ERROR))// & & " already has a mesh associated" CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) ELSE @@ -9654,7 +10156,7 @@ SUBROUTINE DECOMPOSITION_NODE_DOMAIN_GET(DECOMPOSITION,USER_NODE_NUMBER,MESH_COM TYPE(VARYING_STRING) :: LOCAL_ERROR INTEGER(INTG) :: GLOBAL_NODE_NUMBER TYPE(TREE_NODE_TYPE), POINTER :: TREE_NODE - TYPE(MESH_NODES_TYPE), POINTER :: MESH_NODES + TYPE(MeshNodesType), POINTER :: MESH_NODES TYPE(DOMAIN_TYPE), POINTER :: MESH_DOMAIN CALL ENTERS("DECOMPOSITION_NODE_DOMAIN_GET",ERR,ERROR,*999) @@ -9667,49 +10169,49 @@ SUBROUTINE DECOMPOSITION_NODE_DOMAIN_GET(DECOMPOSITION,USER_NODE_NUMBER,MESH_COM IF(ASSOCIATED(MESH)) THEN MESH_TOPOLOGY=>MESH%TOPOLOGY(MESH_COMPONENT_NUMBER)%PTR IF(ASSOCIATED(MESH_TOPOLOGY)) THEN - MESH_NODES=>MESH_TOPOLOGY%NODES + MESH_NODES=>MESH_TOPOLOGY%nodes IF(ASSOCIATED(MESH_NODES)) THEN NULLIFY(TREE_NODE) - CALL TREE_SEARCH(MESH_NODES%NODES_TREE,USER_NODE_NUMBER,TREE_NODE,ERR,ERROR,*999) + CALL TREE_SEARCH(MESH_NODES%nodesTree,USER_NODE_NUMBER,TREE_NODE,ERR,ERROR,*999) IF(ASSOCIATED(TREE_NODE)) THEN - CALL TREE_NODE_VALUE_GET(MESH_NODES%NODES_TREE,TREE_NODE,GLOBAL_NODE_NUMBER,ERR,ERROR,*999) - IF(GLOBAL_NODE_NUMBER>0.AND.GLOBAL_NODE_NUMBER<=MESH_TOPOLOGY%NODES%NUMBER_OF_NODES) THEN + CALL TREE_NODE_VALUE_GET(MESH_NODES%nodesTree,TREE_NODE,GLOBAL_NODE_NUMBER,ERR,ERROR,*999) + IF(GLOBAL_NODE_NUMBER>0.AND.GLOBAL_NODE_NUMBER<=MESH_TOPOLOGY%NODES%numberOfNodes) THEN IF(MESH_COMPONENT_NUMBER>0.AND.MESH_COMPONENT_NUMBER<=MESH%NUMBER_OF_COMPONENTS) THEN MESH_DOMAIN=>DECOMPOSITION%DOMAIN(MESH_COMPONENT_NUMBER)%PTR IF(ASSOCIATED(MESH_DOMAIN)) THEN DOMAIN_NUMBER=MESH_DOMAIN%NODE_DOMAIN(GLOBAL_NODE_NUMBER) ELSE - CALL FLAG_ERROR("Decomposition domain is not associated",ERR,ERROR,*999) + CALL FLAG_ERROR("Decomposition domain is not associated.",ERR,ERROR,*999) ENDIF ELSE - LOCAL_ERROR="Mesh Component number "//TRIM(NUMBER_TO_VSTRING(MESH_COMPONENT_NUMBER,"*",ERR,ERROR))// & + LOCAL_ERROR="Mesh Component number "//TRIM(NumberToVString(MESH_COMPONENT_NUMBER,"*",ERR,ERROR))// & & " is invalid. The limits are 1 to "// & - & TRIM(NUMBER_TO_VSTRING(MESH%NUMBER_OF_COMPONENTS,"*",ERR,ERROR)) + & TRIM(NumberToVString(MESH%NUMBER_OF_COMPONENTS,"*",ERR,ERROR))//"." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) ENDIF ELSE - LOCAL_ERROR="Global element number found "//TRIM(NUMBER_TO_VSTRING(GLOBAL_NODE_NUMBER,"*",ERR,ERROR))// & + LOCAL_ERROR="Global element number found "//TRIM(NumberToVString(GLOBAL_NODE_NUMBER,"*",ERR,ERROR))// & & " is invalid. The limits are 1 to "// & - & TRIM(NUMBER_TO_VSTRING(MESH_TOPOLOGY%NODES%NUMBER_OF_NODES,"*",ERR,ERROR)) + & TRIM(NumberToVString(MESH_TOPOLOGY%NODES%numberOfNodes,"*",ERR,ERROR))//"." CALL FLAG_ERROR(LOCAL_ERROR,ERR,ERROR,*999) ENDIF ELSE - CALL FLAG_ERROR("Decomposition mesh node corresponding to user number not found",ERR,ERROR,*999) + CALL FLAG_ERROR("Decomposition mesh node corresponding to user number not found.",ERR,ERROR,*999) ENDIF ELSE - CALL FLAG_ERROR("Decomposition mesh nodes are not associated",ERR,ERROR,*999) + CALL FLAG_ERROR("Decomposition mesh nodes are not associated.",ERR,ERROR,*999) ENDIF ELSE - CALL FLAG_ERROR("Decomposition mesh topology is not associated",ERR,ERROR,*999) + CALL FLAG_ERROR("Decomposition mesh topology is not associated.",ERR,ERROR,*999) ENDIF ELSE - CALL FLAG_ERROR("Decomposition mesh is not associated",ERR,ERROR,*999) + CALL FLAG_ERROR("Decomposition mesh is not associated.",ERR,ERROR,*999) ENDIF ELSE - CALL FLAG_ERROR("Decomposition has not been finished",ERR,ERROR,*999) + CALL FLAG_ERROR("Decomposition has not been finished.",ERR,ERROR,*999) ENDIF ELSE - CALL FLAG_ERROR("Decomposition is not associated",ERR,ERROR,*999) + CALL FLAG_ERROR("Decomposition is not associated.",ERR,ERROR,*999) ENDIF CALL EXITS("DECOMPOSITION_NODE_DOMAIN_GET") @@ -9890,8 +10392,8 @@ SUBROUTINE MESH_USER_NUMBER_TO_MESH( USER_NUMBER, REGION, MESH, ERR, ERROR, * ) NULLIFY( MESH ) CALL MESH_USER_NUMBER_FIND( USER_NUMBER, REGION, MESH, ERR, ERROR, *999 ) IF( .NOT.ASSOCIATED( MESH ) ) THEN - LOCAL_ERROR = "A mesh with an user number of "//TRIM(NUMBER_TO_VSTRING( USER_NUMBER, "*", ERR, ERROR ))// & - & " does not exist on region number "//TRIM(NUMBER_TO_VSTRING( REGION%USER_NUMBER, "*", ERR, ERROR ))//"." + LOCAL_ERROR = "A mesh with an user number of "//TRIM(NumberToVString( USER_NUMBER, "*", ERR, ERROR ))// & + & " does not exist on region number "//TRIM(NumberToVString( REGION%USER_NUMBER, "*", ERR, ERROR ))//"." CALL FLAG_ERROR( LOCAL_ERROR, ERR, ERROR, *999 ) ENDIF @@ -9924,8 +10426,8 @@ SUBROUTINE DECOMPOSITION_USER_NUMBER_TO_DECOMPOSITION( USER_NUMBER, MESH, DECOMP NULLIFY( DECOMPOSITION ) CALL DECOMPOSITION_USER_NUMBER_FIND( USER_NUMBER, MESH, DECOMPOSITION, ERR, ERROR, *999 ) IF( .NOT.ASSOCIATED( DECOMPOSITION ) ) THEN - LOCAL_ERROR = "A decomposition with an user number of "//TRIM(NUMBER_TO_VSTRING( USER_NUMBER, "*", ERR, ERROR ))// & - & " does not exist on mesh number "//TRIM(NUMBER_TO_VSTRING( MESH%USER_NUMBER, "*", ERR, ERROR ))//"." + LOCAL_ERROR = "A decomposition with an user number of "//TRIM(NumberToVString( USER_NUMBER, "*", ERR, ERROR ))// & + & " does not exist on mesh number "//TRIM(NumberToVString( MESH%USER_NUMBER, "*", ERR, ERROR ))//"." CALL FLAG_ERROR( LOCAL_ERROR, ERR, ERROR, *999 ) ENDIF diff --git a/src/node_routines.f90 b/src/node_routines.f90 index c2e219db..91d21e00 100644 --- a/src/node_routines.f90 +++ b/src/node_routines.f90 @@ -98,7 +98,7 @@ MODULE NODE_ROUTINES PUBLIC NODES_USER_NUMBER_GET,NODES_USER_NUMBER_SET - PUBLIC Nodes_UserNumbersAllSet + PUBLIC NodesUserNumbersAllSet !PUBLIC NODES_NUMBER_OF_VERSIONS_SET @@ -858,7 +858,7 @@ END SUBROUTINE NODES_USER_NUMBER_SET ! !>Changes/sets the user numbers for all nodes. \see OPENCMISS::CMISSNodes_UserNumbersAllSet - SUBROUTINE Nodes_UserNumbersAllSet(nodes,userNumbers,err,error,*) + SUBROUTINE NodesUserNumbersAllSet(nodes,userNumbers,err,error,*) !Argument variables TYPE(NODES_TYPE), POINTER :: nodes !Contains information on a mesh elements defined in a mesh TYPE CMISSMeshElementsType PRIVATE - TYPE(MeshComponentElementsType), POINTER :: MESH_ELEMENTS + TYPE(MeshElementsType), POINTER :: MESH_ELEMENTS END TYPE CMISSMeshElementsType !>Contains information on an embedded mesh @@ -259,6 +259,12 @@ MODULE OPENCMISS TYPE(MESH_EMBEDDING_TYPE), POINTER :: MESH_EMBEDDING END TYPE CMISSMeshEmbeddingType + !>Contains information on a mesh nodes defined in a mesh + TYPE CMISSMeshNodesType + PRIVATE + TYPE(MeshNodesType), POINTER :: meshNodes + END TYPE CMISSMeshNodesType + !>Contains information on the nodes defined on a region. TYPE CMISSNodesType PRIVATE @@ -376,6 +382,8 @@ MODULE OPENCMISS PUBLIC CMISSMeshElementsType,CMISSMeshElements_Finalise,CMISSMeshElements_Initialise + PUBLIC CMISSMeshNodesType,CMISSMeshNodes_Finalise,CMISSMeshNodes_Initialise + PUBLIC CMISSNodesType,CMISSNodes_Finalise,CMISSNodes_Initialise PUBLIC CMISSProblemType,CMISSProblem_Finalise,CMISSProblem_Initialise @@ -4874,7 +4882,6 @@ MODULE OPENCMISS MODULE PROCEDURE CMISSMeshElements_LocalElementNodeVersionSetObj END INTERFACE !CMISSMeshElements_LocalElementNodeVersionSet - !>Returns the element user number for an element in a mesh. INTERFACE CMISSMeshElements_UserNumberGet MODULE PROCEDURE CMISSMeshElements_UserNumberGetNumber @@ -4905,7 +4912,31 @@ MODULE OPENCMISS MODULE PROCEDURE CMISSMesh_ElementExistsObj END INTERFACE !CMISSMesh_ElementExists - !>Returns the domain for a given element in a decomposition of a mesh. + !>Get the mesh nodes belonging to a mesh component. + INTERFACE CMISSMesh_NodesGet + MODULE PROCEDURE CMISSMesh_NodesGetNumber + MODULE PROCEDURE CMISSMesh_NodesGetObj + END INTERFACE CMISSMesh_NodesGet + + !>Returns the number of derivatives for a node in a mesh. + INTERFACE CMISSMeshNodes_NumberOfDerivativesGet + MODULE PROCEDURE CMISSMeshNodes_NumberOfDerivativesGetNumber + MODULE PROCEDURE CMISSMeshNodes_NumberOfDerivativesGetObj + END INTERFACE CMISSMeshNodes_NumberOfDerivativesGet + + !>Returns the derivatives for a node in a mesh. + INTERFACE CMISSMeshNodes_DerivativesGet + MODULE PROCEDURE CMISSMeshNodes_DerivativesGetNumber + MODULE PROCEDURE CMISSMeshNodes_DerivativesGetObj + END INTERFACE CMISSMeshNodes_DerivativesGet + + !>Returns the number of versions for a derivative at a node in a mesh. + INTERFACE CMISSMeshNodes_NumberOfVersionsGet + MODULE PROCEDURE CMISSMeshNodes_NumberOfVersionsGetNumber + MODULE PROCEDURE CMISSMeshNodes_NumberOfVersionsGetObj + END INTERFACE CMISSMeshNodes_NumberOfVersionsGet + + !>Returns the domain for a given element in a decomposition of a mesh. INTERFACE CMISSDecomposition_NodeDomainGet MODULE PROCEDURE CMISSDecomposition_NodeDomainGetNumber MODULE PROCEDURE CMISSDecomposition_NodeDomainGetObj @@ -4979,8 +5010,14 @@ MODULE OPENCMISS PUBLIC CMISSMeshElements_UserNumbersAllSet + PUBLIC CMISSMeshNodes_NumberOfDerivativesGet,CMISSMeshNodes_DerivativesGet + + PUBLIC CMISSMeshNodes_NumberOfVersionsGet + PUBLIC CMISSMesh_ElementsGet + PUBLIC CMISSMesh_NodesGet + PUBLIC CMISSMesh_NodeExists,CMISSMesh_ElementExists PUBLIC CMISSMesh_SurroundingElementsCalculateSet @@ -8219,6 +8256,57 @@ END SUBROUTINE CMISSMeshElements_Initialise !================================================================================================================================ ! + !>Finalises a CMISSMeshNodesType object. + SUBROUTINE CMISSMeshNodes_Finalise(CMISSMeshNodes,err) + + !Argument variables + TYPE(CMISSMeshNodesType), INTENT(OUT) :: CMISSMeshNodes !Initialises a CMISSMeshNodesType object. + SUBROUTINE CMISSMeshNodes_Initialise(CMISSMeshNodes,err) + + !Argument variables + TYPE(CMISSMeshNodesType), INTENT(OUT) :: CMISSMeshNodes !Finalises a CMISSNodesType object. SUBROUTINE CMISSNodes_Finalise(CMISSNodes,err) @@ -41012,7 +41100,7 @@ SUBROUTINE CMISSDecomposition_TopologyDataProjectionCalculateObj(decomposition,e CALL ENTERS("CMISSDecomposition_TopologyDataProjectionCalculateObj",err,error,*999) - CALL DecompositionTopology_DataProjectionCalculate(decomposition%DECOMPOSITION%TOPOLOGY,err,error,*999) + CALL DecompositionTopologyDataProjectionCalculate(decomposition%DECOMPOSITION%TOPOLOGY,err,error,*999) #ifdef TAUPROF CALL TAU_STATIC_PHASE_STOP('CMISSDecomposition_TopologyDataProjectionCalculateObj',err,error,*999) @@ -41045,7 +41133,7 @@ SUBROUTINE CMISSDecomposition_TopologyElementDataPointLocalNumberGetObj(decompos CALL ENTERS("CMISSDecomposition_TopologyElementDataPointLocalNumberGetObj",err,error,*999) - CALL DecompositionTopology_ElementDataPointLocalNumberGet(decomposition%DECOMPOSITION%TOPOLOGY,elementNumber,dataPointIndex, & + CALL DecompositionTopologyElementDataPointLocalNumberGet(decomposition%DECOMPOSITION%TOPOLOGY,elementNumber,dataPointIndex, & & dataPointLocalNumber,err,error,*999) #ifdef TAUPROF @@ -41079,7 +41167,7 @@ SUBROUTINE CMISSDecomposition_TopologyElementDataPointUserNumberGetObj(decomposi CALL ENTERS("CMISSDecomposition_TopologyElementDataPointUserNumberGetObj",err,error,*999) - CALL DecompositionTopology_ElementDataPointUserNumberGet(decomposition%DECOMPOSITION%TOPOLOGY,elementNumber,dataPointIndex, & + CALL DecompositionTopologyElementDataPointUserNumberGet(decomposition%DECOMPOSITION%TOPOLOGY,elementNumber,dataPointIndex, & & dataPointUserNumber,err,error,*999) #ifdef TAUPROF @@ -41111,7 +41199,7 @@ SUBROUTINE CMISSDecomposition_TopologyNumberOfElementDataPointsGetObj(decomposit CALL ENTERS("CMISSDecomposition_TopologyNumberOfElementDataPointsGetObj",err,error,*999) - CALL DecompositionTopology_NumberOfElementDataPointsGet(decomposition%DECOMPOSITION%TOPOLOGY,elementNumber, & + CALL DecompositionTopologyNumberOfElementDataPointsGet(decomposition%DECOMPOSITION%TOPOLOGY,elementNumber, & & numberOfDataPoints,err,error,*999) #ifdef TAUPROF @@ -42930,7 +43018,7 @@ SUBROUTINE CMISSMesh_TopologyDataPointsCalculateProjectionRegionNumber(regionUse IF(ASSOCIATED(REGION)) THEN CALL MESH_USER_NUMBER_FIND(MeshUserNumber,REGION,MESH,Err,ERROR,*999) IF(ASSOCIATED(MESH)) THEN - CALL Mesh_TopologyDataPointsCalculateProjection(MESH,DataProjection%DATA_PROJECTION,Err,ERROR,*999) + CALL MeshTopologyDataPointsCalculateProjection(MESH,DataProjection%DATA_PROJECTION,Err,ERROR,*999) ELSE LOCAL_ERROR="A mesh with an user number of "//TRIM(NUMBER_TO_VSTRING(MeshUserNumber,"*",Err,ERROR))// & & " does not exist on the region with an user number of "//TRIM(NUMBER_TO_VSTRING(regionUserNumber,"*",Err,ERROR))//"." @@ -42981,7 +43069,7 @@ SUBROUTINE CMISSMesh_TopologyDataPointsCalculateProjectionInterfaceNumber(parent IF(ASSOCIATED(INTERFACE)) THEN CALL MESH_USER_NUMBER_FIND(MeshUserNumber,INTERFACE,MESH,Err,ERROR,*999) IF(ASSOCIATED(MESH)) THEN - CALL Mesh_TopologyDataPointsCalculateProjection(MESH,DataProjection%DATA_PROJECTION,Err,ERROR,*999) + CALL MeshTopologyDataPointsCalculateProjection(MESH,DataProjection%DATA_PROJECTION,Err,ERROR,*999) ELSE LOCAL_ERROR="A mesh with an user number of "//TRIM(NUMBER_TO_VSTRING(MeshUserNumber,"*",Err,ERROR))// & & " does not exist on the region with an user number of "//TRIM(NUMBER_TO_VSTRING(parentregionUserNumber, & @@ -43023,7 +43111,7 @@ SUBROUTINE CMISSMesh_TopologyDataPointsCalculateProjectionObj(Mesh,DataProjectio CALL ENTERS("CMISSMesh_TopologyDataPointsCalculateProjectionObj",Err,ERROR,*999) - CALL Mesh_TopologyDataPointsCalculateProjection(Mesh%MESH,DataProjection%DATA_PROJECTION,Err,ERROR,*999) + CALL MeshTopologyDataPointsCalculateProjection(Mesh%MESH,DataProjection%DATA_PROJECTION,Err,ERROR,*999) CALL EXITS("CMISSMesh_TopologyDataPointsCalculateProjectionObj") RETURN @@ -43048,7 +43136,7 @@ SUBROUTINE CMISSMeshElements_CreateFinishNumber(regionUserNumber,meshUserNumber, INTEGER(INTG), INTENT(OUT) :: err !Returns the mesh elements for a mesh component on a mesh identified by an - !user number. + !>Returns the mesh elements for a mesh component on a mesh identified by an user number. SUBROUTINE CMISSMesh_ElementsGetNumber(regionUserNumber,meshUserNumber,meshComponentNumber,meshElements,err) !Argument variables @@ -43291,7 +43378,7 @@ SUBROUTINE CMISSMeshElements_BasisGetNumber(regionUserNumber,meshUserNumber,mesh !Local variables TYPE(BASIS_TYPE), POINTER :: BASIS TYPE(MESH_TYPE), POINTER :: MESH - TYPE(MeshComponentElementsType), POINTER :: MESH_ELEMENTS + TYPE(MeshElementsType), POINTER :: MESH_ELEMENTS TYPE(REGION_TYPE), POINTER :: REGION TYPE(VARYING_STRING) :: LOCAL_ERROR @@ -43382,7 +43469,7 @@ SUBROUTINE CMISSMeshElements_BasisSetNumber(regionUserNumber,meshUserNumber,mesh !Local variables TYPE(BASIS_TYPE), POINTER :: BASIS TYPE(MESH_TYPE), POINTER :: MESH - TYPE(MeshComponentElementsType), POINTER :: MESH_ELEMENTS + TYPE(MeshElementsType), POINTER :: MESH_ELEMENTS TYPE(REGION_TYPE), POINTER :: REGION TYPE(VARYING_STRING) :: LOCAL_ERROR @@ -43470,7 +43557,7 @@ SUBROUTINE CMISSMeshElements_AdjacentElementGetNumber(regionUserNumber,meshUserN INTEGER(INTG), INTENT(OUT) :: err !Returns the mesh nodes for a mesh component on a mesh identified by an user number. + SUBROUTINE CMISSMesh_NodesGetNumber(regionUserNumber,meshUserNumber,meshComponentNumber,meshNodes,err) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: regionUserNumber !Returns the mesh nodes for a mesh component on a mesh identified by an object. + SUBROUTINE CMISSMesh_NodesGetObj(mesh,meshComponentNumber,meshNodes,err) + + !Argument variables + TYPE(CMISSMeshType), INTENT(IN) :: mesh !Returns the number of derivatives at a node in a mesh identified by an user number. + SUBROUTINE CMISSMeshNodes_NumberOfDerivativesGetNumber(regionUserNumber,meshUserNumber,meshComponentNumber,userNodeNumber, & + & numberOfDerivatives,err) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: regionUserNumber !Returns the number of derivatives for a node in a mesh identified by an object. + SUBROUTINE CMISSMeshNodes_NumberOfDerivativesGetObj(meshNodes,userNodeNumber,numberOfDerivatives,err) + + !Argument variables + TYPE(CMISSMeshNodesType), INTENT(IN) :: meshNodes!Returns the derivatives at a node in a mesh identified by an user number. + SUBROUTINE CMISSMeshNodes_DerivativesGetNumber(regionUserNumber,meshUserNumber,meshComponentNumber,userNodeNumber, & + & derivatives,err) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: regionUserNumber !Returns the derivatives for a node in a mesh identified by an object. + SUBROUTINE CMISSMeshNodes_DerivativesGetObj(meshNodes,userNodeNumber,derivatives,err) + + !Argument variables + TYPE(CMISSMeshNodesType), INTENT(IN) :: meshNodes!Returns the number of version at a derivative for a node in a mesh identified by an user number. + SUBROUTINE CMISSMeshNodes_NumberOfVersionsGetNumber(regionUserNumber,meshUserNumber,meshComponentNumber,derivativeNumber, & + & userNodeNumber,numberOfVersions,err) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: regionUserNumber !Returns the number of versions for an node in a mesh identified by an object. + SUBROUTINE CMISSMeshNodes_NumberOfVersionsGetObj(meshNodes,derivativeNumber,userNodeNumber,numberOfVersions,err) + + !Argument variables + TYPE(CMISSMeshNodesType), INTENT(IN) :: meshNodes !Contains information on the dofs for a mesh. - TYPE MESH_DOFS_TYPE - TYPE(MESH_TYPE), POINTER :: MESH !Contains information on the mesh adjacent elements for a xi coordinate TYPE MESH_ADJACENT_ELEMENT_TYPE @@ -398,42 +397,42 @@ MODULE TYPES END TYPE MESH_ELEMENT_TYPE !>Contains the information for the elements of a mesh. - TYPE MeshComponentElementsType - TYPE(MESH_TYPE), POINTER :: MESH !Contains the information for a node derivative of a mesh. - TYPE MESH_NODE_DERIVATIVE_TYPE - INTEGER(INTG) :: NUMBER_OF_VERSIONS !The number of global versions at the node for the mesh. - INTEGER(INTG), ALLOCATABLE :: USER_VERSION_NUMBERS(:) !The user version index of the nk'th global derivative for the node. - INTEGER(INTG), ALLOCATABLE :: DOF_INDEX(:) !The global dof version index (nv) in the domain of the nk'th global derivative for the node. - INTEGER(INTG) :: GLOBAL_DERIVATIVE_INDEX !The global derivative index of the nk'th global derivative for the node. - INTEGER(INTG) :: PARTIAL_DERIVATIVE_INDEX !The partial derivative index (nu) of the nk'th global derivative for the node. Old CMISS name NUNK(nk,nj,np). - END TYPE MESH_NODE_DERIVATIVE_TYPE + TYPE MeshNodeDerivativeType + INTEGER(INTG) :: numberOfVersions !The number of global versions at the node for the mesh. + INTEGER(INTG), ALLOCATABLE :: userVersionNumbers(:) !userVersionNumbers(versionIdx). The user version numbers for the versionIdx'th version for the node. + INTEGER(INTG), ALLOCATABLE :: dofIndex(:) !The global dof version index (nv) in the domain of the nk'th global derivative for the node. + INTEGER(INTG) :: globalDerivativeIndex !The global derivative index of the nk'th global derivative for the node. + INTEGER(INTG) :: partialDerivativeIndex !The partial derivative index (nu) of the nk'th global derivative for the node. Old CMISS name NUNK(nk,nj,np). + END TYPE MeshNodeDerivativeType !>Contains the topology information for a global node of a mesh. - TYPE MESH_NODE_TYPE - INTEGER(INTG) :: MESH_NUMBER !Contains the information for the nodes of a mesh. - TYPE MESH_NODES_TYPE - TYPE(MESH_TYPE), POINTER :: MESH !Contains information on the (global) topology of a mesh. TYPE MeshComponentTopologyType - TYPE(MESH_TYPE), POINTER :: MESH !A buffer type to allow for an array of pointers to a MeshComponentTopologyType. TYPE MeshComponentTopologyPtrType - TYPE(MeshComponentTopologyType), POINTER :: PTR !Embedded mesh types @@ -513,7 +512,7 @@ MODULE TYPES INTEGER(INTG) :: NUMBER_OF_EMBEDDED_MESHES !