From 85a3134e1dcc512f71aefa992f9ba15ede28309a Mon Sep 17 00:00:00 2001 From: Chris Bradley Date: Fri, 17 May 2019 14:43:56 +1200 Subject: [PATCH 01/17] Give the ability to set the projected xi location for data projections. --- src/data_projection_access_routines.F90 | 13 +- src/data_projection_routines.F90 | 2105 +++++++++++++---------- src/interface_routines.F90 | 2 +- src/mesh_routines.F90 | 4 +- src/opencmiss_iron.F90 | 254 ++- src/types.F90 | 4 +- 6 files changed, 1462 insertions(+), 920 deletions(-) diff --git a/src/data_projection_access_routines.F90 b/src/data_projection_access_routines.F90 index 1d831ea5..04d94b93 100644 --- a/src/data_projection_access_routines.F90 +++ b/src/data_projection_access_routines.F90 @@ -62,10 +62,11 @@ MODULE DataProjectionAccessRoutines !> \see DataProjectionRoutines,OPENCMISS_DataProjectionExitTags !>@{ INTEGER(INTG), PARAMETER :: DATA_PROJECTION_CANCELLED=0 !@} !Module types @@ -74,8 +75,8 @@ MODULE DataProjectionAccessRoutines !Interfaces - PUBLIC DATA_PROJECTION_CANCELLED,DATA_PROJECTION_EXIT_TAG_CONVERGED,DATA_PROJECTION_EXIT_TAG_BOUNDS, & - & DATA_PROJECTION_EXIT_TAG_MAX_ITERATION,DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + PUBLIC DATA_PROJECTION_CANCELLED,DATA_PROJECTION_USER_SPECIFIED,DATA_PROJECTION_EXIT_TAG_CONVERGED, & + & DATA_PROJECTION_EXIT_TAG_BOUNDS,DATA_PROJECTION_EXIT_TAG_MAX_ITERATION,DATA_PROJECTION_EXIT_TAG_NO_ELEMENT PUBLIC DataProjection_DataPointsGet diff --git a/src/data_projection_routines.F90 b/src/data_projection_routines.F90 index 8792f40f..b3082bc4 100644 --- a/src/data_projection_routines.F90 +++ b/src/data_projection_routines.F90 @@ -182,7 +182,11 @@ MODULE DataProjectionRoutines PUBLIC DataProjection_ResultDistanceGet - PUBLIC DataProjection_ResultElementNumberGet,DataProjection_ResultElementFaceNumberGet,DataProjection_ResultElementLineNumberGet + PUBLIC DataProjection_ResultElementNumberGet,DataProjection_ResultElementNumberSet + + PUBLIC DataProjection_ResultElementFaceNumberGet,DataProjection_ResultElementFaceNumberSet + + PUBLIC DataProjection_ResultElementLineNumberGet,DataProjection_ResultElementLineNumberSet PUBLIC DataProjection_ResultExitTagGet @@ -893,6 +897,35 @@ END SUBROUTINE DataProjection_DataProjectionCandidatesInitialise !================================================================================================================================ ! + !>Cancels a data projection result + SUBROUTINE DataProjection_DataProjectionResultCancel(dataProjectionResult,err,error,*) + + !Argument variables + TYPE(DataProjectionResultType) :: dataProjectionResult !Finalises the data projection result and deallocates all memory SUBROUTINE DataProjection_DataProjectionResultFinalise(dataProjectionResult,err,error,*) @@ -906,7 +939,8 @@ SUBROUTINE DataProjection_DataProjectionResultFinalise(dataProjectionResult,err, dataProjectionResult%userNumber=0 dataProjectionResult%distance=0.0 - dataProjectionResult%elementNumber=0 + dataProjectionResult%elementLocalNumber=0 + dataProjectionResult%elementGlobalNumber=0 dataProjectionResult%elementLineFaceNumber=0 dataProjectionResult%exitTag=0 IF(ALLOCATED(dataProjectionResult%xi)) DEALLOCATE(dataProjectionResult%xi) @@ -938,7 +972,8 @@ SUBROUTINE DataProjection_DataProjectionResultInitialise(dataProjectionResult,er dataProjectionResult%userNumber=0 dataProjectionResult%distance=0.0_DP - dataProjectionResult%elementNumber=0 + dataProjectionResult%elementLocalNumber=0 + dataProjectionResult%elementGlobalNumber=0 dataProjectionResult%elementLineFaceNumber=0 dataProjectionResult%exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT @@ -1372,7 +1407,7 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection INTEGER(INTG), ALLOCATABLE :: globalNumberOfProjectedPoints(:) INTEGER(INTG) :: MPIClosestDistances,dataProjectionGlobalNumber INTEGER(INTG) :: MPIIError - INTEGER(INTG), ALLOCATABLE :: projectedElement(:),projectedLineFace(:),projectionExitTag(:) + INTEGER(INTG), ALLOCATABLE :: projectedGlobalElement(:),projectedLineFace(:),projectedLocalElement(:),projectionExitTag(:) REAL(DP), ALLOCATABLE :: closestDistances(:,:),globalClosestDistances(:,:) REAL(DP), ALLOCATABLE :: projectedDistance(:,:),projectedXi(:,:),projectionVectors(:,:) INTEGER(INTG) :: elementIdx,elementNumber,lineFaceIdx,lineFaceNumber,computationalNodeIdx,xiIdx,localElementNumber, & @@ -1673,8 +1708,10 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection IF(err/=0) CALL FlagError("Could not allocate all number of projected points.",err,error,*999) ALLOCATE(projectionExitTag(numberOfDataPoints),STAT=err) IF(err/=0) CALL FlagError("Could not allocate projected.",err,error,*999) - ALLOCATE(projectedElement(numberOfDataPoints),STAT=err) - IF(err/=0) CALL FlagError("Could not allocate projected element.",err,error,*999) + ALLOCATE(projectedLocalElement(numberOfDataPoints),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate projected local element.",err,error,*999) + ALLOCATE(projectedGlobalElement(numberOfDataPoints),STAT=err) + IF(err/=0) CALL FlagError("Could not allocate projected global element.",err,error,*999) IF(boundaryProjection) THEN ALLOCATE(projectedLineFace(numberOfDataPoints),STAT=err) IF(err/=0) CALL FlagError("Could not allocate projected sub element.",err,error,*999) @@ -1724,6 +1761,8 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection projectedDistance(2,:)=myComputationalNode !Find the globally closest distance in the current domain DO dataPointIdx=1,numberOfDataPoints + projectionExitTag(dataPointIdx)=dataProjection%dataProjectionResults(dataPointIdx)%exitTag + projectedLocalElement(dataPointIdx)=dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber CALL Sorting_BubbleIndexSort(globalClosestDistances(dataPointIdx,:),sortingIndices1,err,error,*999) sortingIndices1(1:totalNumberOfClosestCandidates)=sortingIndices1(1:totalNumberOfClosestCandidates)- & & globalMPIDisplacements(myComputationalNode+1) !shift the index to current computational node @@ -1747,11 +1786,11 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonLinesEvaluate(dataProjection,interpolatedPoint, & & dataPoints%dataPoints(dataPointIdx)%position,closestElements( & & dataPointIdx,1:numberOfClosestCandidates),closestLinesFaces(dataPointIdx,1: & - & numberOfClosestCandidates),projectionExitTag(dataPointIdx),projectedElement(dataPointIdx), & + & numberOfClosestCandidates),projectionExitTag(dataPointIdx),projectedLocalElement(dataPointIdx), & & projectedLineFace(dataPointIdx),projectedDistance(1,dataPointIdx),projectedXi(:,dataPointIdx), & & projectionVectors(:,dataPointIdx),err,error,*999) !Map the element number to global number - projectedElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedElement(dataPointIdx)) + projectedGlobalElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedLocalElement(dataPointIdx)) ENDIF ENDDO !dataPointIdx CASE(DATA_PROJECTION_BOUNDARY_FACES_PROJECTION_TYPE) @@ -1762,11 +1801,11 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonFacesEvaluate(dataProjection,interpolatedPoint, & & dataPoints%dataPoints(dataPointIdx)%position,closestElements( & & dataPointIdx,1:numberOfClosestCandidates),closestLinesFaces(dataPointIdx, & - & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx),projectedElement(dataPointIdx), & + & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx),projectedLocalElement(dataPointIdx), & & projectedLineFace(dataPointIdx),projectedDistance(1,dataPointIdx),projectedXi(:,dataPointIdx), & & projectionVectors(:,dataPointIdx),err,error,*999) !Map the element number to global number - projectedElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedElement(dataPointIdx)) + projectedGlobalElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedLocalElement(dataPointIdx)) ENDIF ENDDO !dataPointIdx CASE(DATA_PROJECTION_ALL_ELEMENTS_PROJECTION_TYPE) @@ -1779,10 +1818,10 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx, & & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx), & - & projectedElement(dataPointIdx),projectedDistance(1,dataPointIdx), & + & projectedLocalElement(dataPointIdx),projectedDistance(1,dataPointIdx), & & projectedXi(:,dataPointIdx),projectionVectors(:,dataPointIdx),err,error,*999) !Map the element number to global number - projectedElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedElement(dataPointIdx)) + projectedGlobalElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedLocalElement(dataPointIdx)) ENDIF ENDDO !dataPointIdx CASE(2) !2D element @@ -1792,11 +1831,11 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonElementsEvaluate_2(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx, & & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx), & - & projectedElement(dataPointIdx),projectedDistance(1,dataPointIdx), & + & projectedLocalElement(dataPointIdx),projectedDistance(1,dataPointIdx), & & projectedXi(:,dataPointIdx),projectionVectors(:,dataPointIdx), & & err,error,*999) !Map the element number to global number - projectedElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedElement(dataPointIdx)) + projectedGlobalElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedLocalElement(dataPointIdx)) ENDIF ENDDO !dataPointIdx CASE(3) !3D element @@ -1806,10 +1845,10 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonElementsEvaluate_3(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx, & & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx), & - & projectedElement(dataPointIdx),projectedDistance(1,dataPointIdx), & + & projectedLocalElement(dataPointIdx),projectedDistance(1,dataPointIdx), & & projectedXi(:,dataPointIdx),projectionVectors(:,dataPointIdx),err,error,*999) !Map the element number to global number - projectedElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedElement(dataPointIdx)) + projectedGlobalElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedLocalElement(dataPointIdx)) ENDIF ENDDO !dataPointIdx CASE DEFAULT @@ -1841,8 +1880,12 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection & globalNumberOfProjectedPoints(computationalNodeIdx) ENDDO !computationalNodeIdx !Shares minimum projection information between all domains - CALL MPI_ALLGATHERV(projectedElement(sortingIndices2(startIdx:finishIdx)),globalNumberOfProjectedPoints( & - & myComputationalNode+1),MPI_INTEGER,projectedElement,globalNumberOfProjectedPoints, & + CALL MPI_ALLGATHERV(projectedLocalElement(sortingIndices2(startIdx:finishIdx)),globalNumberOfProjectedPoints( & + & myComputationalNode+1),MPI_INTEGER,projectedLocalElement,globalNumberOfProjectedPoints, & + & globalMPIDisplacements,MPI_INTEGER,computationalEnvironment%mpiCommunicator,MPIIError) !projectedElement + CALL MPI_ERROR_CHECK("MPI_ALLGATHERV",MPIIError,err,error,*999) + CALL MPI_ALLGATHERV(projectedGlobalElement(sortingIndices2(startIdx:finishIdx)),globalNumberOfProjectedPoints( & + & myComputationalNode+1),MPI_INTEGER,projectedGlobalElement,globalNumberOfProjectedPoints, & & globalMPIDisplacements,MPI_INTEGER,computationalEnvironment%mpiCommunicator,MPIIError) !projectedElement CALL MPI_ERROR_CHECK("MPI_ALLGATHERV",MPIIError,err,error,*999) IF(boundaryProjection) THEN @@ -1871,7 +1914,10 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection !Assign projection information to projected points DO dataPointIdx=1,numberOfDataPoints dataProjection%dataProjectionResults(sortingIndices2(dataPointIdx))%exitTag=projectionExitTag(dataPointIdx) - dataProjection%dataProjectionResults(sortingIndices2(dataPointIdx))%elementNumber=projectedElement(dataPointIdx) + dataProjection%dataProjectionResults(sortingIndices2(dataPointIdx))%elementLocalNumber= & + & projectedLocalElement(dataPointIdx) + dataProjection%dataProjectionResults(sortingIndices2(dataPointIdx))%elementGlobalNumber= & + & projectedGlobalElement(dataPointIdx) dataProjection%dataProjectionResults(dataPointIdx)%DISTANCE=projectedDistance(1,dataPointIdx) dataProjection%dataProjectionResults(sortingIndices2(dataPointIdx))%xi(1:dataProjection%numberOfXi)= & & projectedXi(1:dataProjection%numberOfXi,dataPointIdx) @@ -1880,7 +1926,8 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection ENDDO !dataPointIdx projectedXi(:,sortingIndices2)=projectedXi projectionVectors(:, sortingIndices2)=projectionVectors - projectedElement(sortingIndices2)=projectedElement + projectedLocalElement(sortingIndices2)=projectedLocalElement + projectedGlobalElement(sortingIndices2)=projectedGlobalElement IF(dataProjection%projectionType==DATA_PROJECTION_BOUNDARY_LINES_PROJECTION_TYPE) THEN DO dataPointIdx=1,numberOfDataPoints dataProjection%dataProjectionResults(sortingIndices2(dataPointIdx))%elementLineFaceNumber=projectedLineFace(dataPointIdx) @@ -1899,10 +1946,12 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonLinesEvaluate(dataProjection,interpolatedPoint,dataPoints%dataPoints( & & dataPointIdx)%position,closestElements(dataPointIdx,:),closestLinesFaces(dataPointIdx,:), & & dataProjection%dataProjectionResults(dataPointIdx)%exitTag,dataProjection% & - & dataProjectionResults(dataPointIdx)%elementNumber,dataProjection% & + & dataProjectionResults(dataPointIdx)%elementLocalNumber,dataProjection% & & dataProjectionResults(dataPointIdx)%elementLineFaceNumber,dataProjection%dataProjectionResults( & - & dataPointIdx)%DISTANCE,dataProjection%dataProjectionResults(dataPointIdx)%xi, & + & dataPointIdx)%distance,dataProjection%dataProjectionResults(dataPointIdx)%xi, & & dataProjection%dataProjectionResults(dataPointIdx)%projectionVector,err,error,*999) + dataProjection%dataProjectionResults(dataPointIdx)%elementGlobalNumber= & + & dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber ENDDO !dataPointIdx CASE(DATA_PROJECTION_BOUNDARY_FACES_PROJECTION_TYPE) !Newton project to closest faces, and find miminum projection @@ -1910,10 +1959,12 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonFacesEvaluate(dataProjection,interpolatedPoint,dataPoints%dataPoints( & & dataPointIdx)%position,closestElements(dataPointIdx,:),closestLinesFaces(dataPointIdx,:), & & dataProjection%dataProjectionResults(dataPointIdx)%exitTag,dataProjection% & - & dataProjectionResults(dataPointIdx)%elementNumber,dataProjection% & + & dataProjectionResults(dataPointIdx)%elementLocalNumber,dataProjection% & & dataProjectionResults(dataPointIdx)%elementLineFaceNumber,dataProjection%dataProjectionResults( & & dataPointIdx)%DISTANCE,dataProjection%dataProjectionResults(dataPointIdx)%xi, & & dataProjection%dataProjectionResults(dataPointIdx)%projectionVector,err,error,*999) + dataProjection%dataProjectionResults(dataPointIdx)%elementGlobalNumber= & + & dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber ENDDO !dataPointIdx CASE(DATA_PROJECTION_ALL_ELEMENTS_PROJECTION_TYPE) !Newton project to closest elements, and find miminum projection @@ -1923,30 +1974,36 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx,:),dataProjection% & & dataProjectionResults(dataPointIdx)%exitTag,dataProjection%dataProjectionResults( & - & dataPointIdx)%elementNumber,dataProjection%dataProjectionResults(dataPointIdx)%distance, & + & dataPointIdx)%elementLocalNumber,dataProjection%dataProjectionResults(dataPointIdx)%distance, & & dataProjection%dataProjectionResults(dataPointIdx)%xi, & & dataProjection%dataProjectionResults(dataPointIdx)%projectionVector,& & err,error,*999) - ENDDO !dataPointIdx + dataProjection%dataProjectionResults(dataPointIdx)%elementGlobalNumber= & + & dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber + ENDDO !dataPointIdx CASE(2) !2D mesh DO dataPointIdx=1,numberOfDataPoints CALL DataProjection_NewtonElementsEvaluate_2(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx,:),dataProjection% & & dataProjectionResults(dataPointIdx)%exitTag,dataProjection%dataProjectionResults( & - & dataPointIdx)%elementNumber,dataProjection%dataProjectionResults(dataPointIdx)%distance, & + & dataPointIdx)%elementLocalNumber,dataProjection%dataProjectionResults(dataPointIdx)%distance, & & dataProjection%dataProjectionResults(dataPointIdx)%xi, & & dataProjection%dataProjectionResults(dataPointIdx)%projectionVector, & & err,error,*999) + dataProjection%dataProjectionResults(dataPointIdx)%elementGlobalNumber= & + & dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber ENDDO !dataPointIdx CASE(3) !3D mesh DO dataPointIdx=1,numberOfDataPoints CALL DataProjection_NewtonElementsEvaluate_3(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx,:),dataProjection% & & dataProjectionResults(dataPointIdx)%exitTag,dataProjection%dataProjectionResults( & - & dataPointIdx)%elementNumber,dataProjection%dataProjectionResults(dataPointIdx)%distance, & + & dataPointIdx)%elementLocalNumber,dataProjection%dataProjectionResults(dataPointIdx)%distance, & & dataProjection%dataProjectionResults(dataPointIdx)%xi, & & dataProjection%dataProjectionResults(dataPointIdx)%projectionVector, & & err,error,*999) + dataProjection%dataProjectionResults(dataPointIdx)%elementGlobalNumber= & + & dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber ENDDO !dataPointIdx CASE DEFAULT localError="The data projection number of xi of "//TRIM(NumberToVString(dataProjection%numberOfXi,"*",err,error))// & @@ -1967,7 +2024,7 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection ENDDO !dataPointIdx ELSE DO dataPointIdx=1,numberOfDataPoints - elementNumber=dataProjection%dataProjectionResults(dataPointIdx)%elementNumber + elementNumber=dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber localLineFaceNumber=dataProjection%dataProjectionResults(dataPointIdx)%elementLineFaceNumber basis=>domainElements%elements(elementNumber)%basis CALL Basis_BoundaryXiToXi(basis,localLineFaceNumber,dataProjection% & @@ -1991,7 +2048,8 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection IF(ALLOCATED(globalMPIDisplacements)) DEALLOCATE(globalMPIDisplacements) IF(ALLOCATED(globalNumberOfProjectedPoints)) DEALLOCATE(globalNumberOfProjectedPoints) IF(ALLOCATED(projectionExitTag)) DEALLOCATE(projectionExitTag) - IF(ALLOCATED(projectedElement)) DEALLOCATE(projectedElement) + IF(ALLOCATED(projectedLocalElement)) DEALLOCATE(projectedLocalElement) + IF(ALLOCATED(projectedGlobalElement)) DEALLOCATE(projectedGlobalElement) IF(ALLOCATED(projectedLineFace)) DEALLOCATE(projectedLineFace) IF(ALLOCATED(projectedDistance)) DEALLOCATE(projectedDistance) IF(ALLOCATED(projectedXi)) DEALLOCATE(projectedXi) @@ -2012,7 +2070,8 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection IF(ALLOCATED(globalMPIDisplacements)) DEALLOCATE(globalMPIDisplacements) IF(ALLOCATED(globalNumberOfProjectedPoints)) DEALLOCATE(globalNumberOfProjectedPoints) IF(ALLOCATED(projectionExitTag)) DEALLOCATE(projectionExitTag) - IF(ALLOCATED(projectedElement)) DEALLOCATE(projectedElement) + IF(ALLOCATED(projectedLocalElement)) DEALLOCATE(projectedLocalElement) + IF(ALLOCATED(projectedGlobalElement)) DEALLOCATE(projectedGlobalElement) IF(ALLOCATED(projectedLineFace)) DEALLOCATE(projectedLineFace) IF(ALLOCATED(projectedDistance)) DEALLOCATE(projectedDistance) IF(ALLOCATED(projectedXi)) DEALLOCATE(projectedXi) @@ -2072,7 +2131,7 @@ SUBROUTINE DataProjection_DataPointsPositionEvaluate(dataProjection,field,fieldV interpolatedPoint=>interpolatedPoints(fieldVariableType)%ptr !Loop through data points DO dataPointIdx=1,dataPoints%numberOfDataPoints - elementNumber=dataProjection%dataProjectionResults(dataPointIdx)%elementNumber + elementNumber=dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber CALL Field_InterpolationParametersElementGet(fieldParameterSetType,elementNumber, & & interpolationParameters(fieldVariableType)%ptr,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) CALL Field_InterpolateXi(NO_PART_DERIV,dataProjection%dataProjectionResults(dataPointIdx)%xi, & @@ -2253,118 +2312,128 @@ SUBROUTINE DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPo IF(.NOT.ASSOCIATED(dataProjection)) CALL FlagError("Data projection is not associated.",err,error,*999) IF(.NOT.dataProjection%dataProjectionFinished) CALL FlagError("Data projection has not been finished.",err,error,*999) - projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT numberOfCoordinates=dataProjection%numberOfCoordinates - relativeTolerance=dataProjection%relativeTolerance - absoluteTolerance=dataProjection%absoluteTolerance - maximumDelta=dataProjection%maximumIterationUpdate - minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small - !Project on each candidate elements - DO elementIdx=1,SIZE(candidateElements,1) - elementNumber=candidateElements(elementIdx) - IF(elementNumber>0) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT - converged=.FALSE. - !Start at half the maximumDelta as we do not know if quadratic model is a good approximation yet - delta=0.5_DP*maximumDelta - CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,elementNumber,interpolatedPoint% & - & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - xi=dataProjection%startingXi - CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector(1:numberOfCoordinates)=dataPointLocation-interpolatedPoint%values(:,NO_PART_DERIV) - functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations !(outer loop) - !Check for bounds [0,1] - IF(ABS(xi(1))absoluteTolerance) THEN !positive: minimum exists - updateXi(1)=-functionGradient/functionHessian - updateXiNorm=DABS(updateXi(1)) - insideRegion=updateXiNorm<=delta - ENDIF !positive - IF(.NOT.insideRegion) THEN !minimum not in the region - updateXi(1)=-DSIGN(delta,functionGradient) - updateXiNorm=delta - ENDIF - IF((bound/=0).AND.(bound>0.EQV.updateXi(1)>0.0_DP)) THEN !projection go out of element bound - exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS - EXIT main_loop - ENDIF - converged=updateXiNorm1.0_DP) THEN - newXi(1)=1.0_DP + IF(projectionExitTag==DATA_PROJECTION_USER_SPECIFIED) THEN + IF(projectionElementNumber==0) CALL FlagError("The projection element number has not been set.",err,error,*999) + CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,projectionElementNumber, & + & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + CALL Field_InterpolateXi(NO_PART_DERIV,projectionXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + projectionVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + projectionDistance=DSQRT(DOT_PRODUCT(projectionVector(1:numberOfCoordinates),projectionVector(1:numberOfCoordinates))) + ELSE + projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + relativeTolerance=dataProjection%relativeTolerance + absoluteTolerance=dataProjection%absoluteTolerance + maximumDelta=dataProjection%maximumIterationUpdate + minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small + !Project on each candidate elements + DO elementIdx=1,SIZE(candidateElements,1) + elementNumber=candidateElements(elementIdx) + IF(elementNumber>0) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + converged=.FALSE. + !Start at half the maximumDelta as we do not know if quadratic model is a good approximation yet + delta=0.5_DP*maximumDelta + CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,elementNumber,interpolatedPoint% & + & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + xi=dataProjection%startingXi + CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation-interpolatedPoint%values(:,NO_PART_DERIV) + functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations !(outer loop) + !Check for bounds [0,1] + IF(ABS(xi(1))absoluteTolerance) THEN !bad model: reduce step size - IF(delta<=minimumDelta) THEN - !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + !functionGradient + functionGradient=-2.0_DP* & + & (DOT_PRODUCT(distanceVector(1:numberOfCoordinates),interpolatedPoint%values(:,FIRST_PART_DERIV))) + !functionHessian + functionHessian=-2.0_DP* & + & (DOT_PRODUCT(distanceVector(1:numberOfCoordinates),interpolatedPoint%values(:,SECOND_PART_DERIV))- & + & DOT_PRODUCT(interpolatedPoint%values(:,FIRST_PART_DERIV),interpolatedPoint%values(:,FIRST_PART_DERIV))) + !A model trust region approach, directly taken from CMISS CLOS22: V = -(H + EIGEN_SHIFT*I)g + !The calculation of EIGEN_SHIFT are only approximated as opposed to the common trust region approach + !(inner loop: adjust region size) usually EXIT at 1 or 2 iterations + DO iterationIdx2=1,dataProjection%maximumNumberOfIterations + insideRegion=.FALSE. + IF(functionHessian>absoluteTolerance) THEN !positive: minimum exists + updateXi(1)=-functionGradient/functionHessian + updateXiNorm=DABS(updateXi(1)) + insideRegion=updateXiNorm<=delta + ENDIF !positive + IF(.NOT.insideRegion) THEN !minimum not in the region + updateXi(1)=-DSIGN(delta,functionGradient) + updateXiNorm=delta + ENDIF + IF((bound/=0).AND.(bound>0.EQV.updateXi(1)>0.0_DP)) THEN !projection go out of element bound + exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.25_DP*delta) - ELSE - predictedReduction=updateXi(1)*(functionGradient+0.5_DP*functionHessian*updateXi(1)) - predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction - IF(predictionAccuracy<0.01_DP) THEN !bad model: reduce region size + converged=updateXiNorm1.0_DP) THEN + newXi(1)=1.0_DP + ENDIF + CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector=dataPointLocation-interpolatedPoint%values(:,NO_PART_DERIV) + newFunctionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !second half of the convergence test + converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN !bad model: reduce step size IF(delta<=minimumDelta) THEN !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.5_DP*delta) - ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN - !good model: increase region size - delta=DMIN1(maximumDelta,2.0_DP*delta) - EXIT + delta=DMAX1(minimumDelta,0.25_DP*delta) ELSE - !ok model: keep the current region size - EXIT + predictedReduction=updateXi(1)*(functionGradient+0.5_DP*functionHessian*updateXi(1)) + predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction + IF(predictionAccuracy<0.01_DP) THEN !bad model: reduce region size + IF(delta<=minimumDelta) THEN + !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + EXIT main_loop + ENDIF + delta=DMAX1(minimumDelta,0.5_DP*delta) + ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN + !good model: increase region size + delta=DMIN1(maximumDelta,2.0_DP*delta) + EXIT + ELSE + !ok model: keep the current region size + EXIT + ENDIF ENDIF + ENDDO !iterationIdx2 (inner loop: adjust region size) + functionValue=newFunctionValue + xi=newXi + IF(converged) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED + EXIT ENDIF - ENDDO !iterationIdx2 (inner loop: adjust region size) - functionValue=newFunctionValue - xi=newXi - IF(converged) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED - EXIT + ENDDO main_loop !iterationIdx1 (outer loop) + IF(exitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT.AND.iterationIdx1>=dataProjection%maximumNumberOfIterations) & + & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION + IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)=dataProjection%maximumNumberOfIterations) & - & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION - IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%domain(meshComponentNumber)%ptr% & & mappings%elements numberOfCoordinates=dataProjection%numberOfCoordinates - relativeTolerance=dataProjection%relativeTolerance - absoluteTolerance=dataProjection%absoluteTolerance - maximumDelta=dataProjection%maximumIterationUpdate - minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small - DO elementIdx=1,SIZE(candidateElements,1) !project on each candidate elements - elementNumber=candidateElements(elementIdx) - IF(elementNumber>0) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT - converged=.FALSE. - !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet - delta=0.5_DP*maximumDelta - CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,elementNumber, & - & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - xi=dataProjection%startingXi - CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & - & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - !Outer loop - main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations - !Check for bounds [0,1] - DO xiIdx=1,2 - IF(ABS(xi(xiIdx))=temp1).AND.(minEigen>absoluteTolerance) - IF(insideRegion) THEN - determinant=minEigen*maxEigen !det(H) - hessianDiagonal(1)=functionHessian(1,1) - hessianDiagonal(2)=functionHessian(2,2) - ELSE - eigenShift=MAX(temp1-minEigen,absoluteTolerance) !shift towards steepest decent - determinant=temp1*(maxEigen+eigenShift) !det(H) - hessianDiagonal(1)=functionHessian(1,1)+eigenShift - hessianDiagonal(2)=functionHessian(2,2)+eigenShift - ENDIF - updateXi(1)=-(hessianDiagonal(2)*functionGradient(1)-functionHessian(1,2)*functionGradient(2))/determinant - updateXi(2)=(functionHessian(1,2)*functionGradient(1)-hessianDiagonal(1)*functionGradient(2))/determinant - updateXiNorm=DSQRT(DOT_PRODUCT(updateXi,updateXi)) - free=.TRUE. - DO xiIdx=1,2 - IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN - !Projection will go out of element bounds - IF(.NOT.free) THEN !both xi are fixed - exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS - EXIT main_loop - ENDIF - free=.FALSE. - fixedXiIdx=xiIdx - ENDIF - ENDDO !xiIdx - IF(free) THEN - !Both xi are free - IF(.NOT.insideRegion) THEN - IF(updateXiNorm>0.0_DP) THEN - updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound - ENDIF - ENDIF - ELSE - !xi are not free - updateXi(fixedXiIdx)=0.0_DP - xiIdx=3-fixedXiIdx - insideRegion=.FALSE. - IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN - !Positive: minimum exists in the unbounded direction - updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) - updateXiNorm=DABS(updateXi(xiIdx)) - insideRegion=updateXiNorm<=delta - ENDIF - IF(.NOT.insideRegion) THEN - !Minimum not in the region - updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) - updateXiNorm=delta + IF(projectionExitTag==DATA_PROJECTION_USER_SPECIFIED) THEN + IF(projectionElementNumber==0) CALL FlagError("The projection element number has not been set.",err,error,*999) + CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,projectionElementNumber, & + & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + CALL Field_InterpolateXi(NO_PART_DERIV,projectionXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + projectionVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + projectionDistance=DSQRT(DOT_PRODUCT(projectionVector(1:numberOfCoordinates),projectionVector(1:numberOfCoordinates))) + ELSE + projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + relativeTolerance=dataProjection%relativeTolerance + absoluteTolerance=dataProjection%absoluteTolerance + maximumDelta=dataProjection%maximumIterationUpdate + minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small + DO elementIdx=1,SIZE(candidateElements,1) !project on each candidate elements + elementNumber=candidateElements(elementIdx) + IF(elementNumber>0) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + converged=.FALSE. + !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet + delta=0.5_DP*maximumDelta + CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,elementNumber, & + & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + xi=dataProjection%startingXi + CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Outer loop + main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations + !Check for bounds [0,1] + DO xiIdx=1,2 + IF(ABS(xi(xiIdx))1.0_DP) THEN - newXi(xiIdx)=1.0_DP - newXi(3-xiIdx)=xi(3-xiIdx)+updateXi(3-xiIdx)*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) + ENDDO !xiIdx + !functionGradient + functionGradient(1)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1))) + functionGradient(2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + !functionHessian + functionHessian(1,1)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1_S1))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1))) + functionHessian(1,2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1_S2))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + functionHessian(2,2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2_S2))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + !A model trust region approach, a Newton step is taken if the minimum lies inside the trust region (delta), + !if not, shift the step towards the steepest descent + temp1=0.5_DP*(functionHessian(1,1)+functionHessian(2,2)) + temp2=DSQRT((0.5_DP*(functionHessian(1,1)-functionHessian(2,2)))**2+functionHessian(1,2)**2) + minEigen=temp1-temp2 + maxEigen=temp1+temp2 + functionGradientNorm=DSQRT(DOT_PRODUCT(functionGradient,functionGradient)) + !Inner loop: adjust region size, usually EXIT at 1 or 2 iterations + DO iterationIdx2=1,dataProjection%maximumNumberOfIterations + temp1=functionGradientNorm/delta + !Estimate if the solution is inside the trust region without calculating a Newton step, this also guarantees + !the Hessian matrix is positive definite + insideRegion=(minEigen>=temp1).AND.(minEigen>absoluteTolerance) + IF(insideRegion) THEN + determinant=minEigen*maxEigen !det(H) + hessianDiagonal(1)=functionHessian(1,1) + hessianDiagonal(2)=functionHessian(2,2) + ELSE + eigenShift=MAX(temp1-minEigen,absoluteTolerance) !shift towards steepest decent + determinant=temp1*(maxEigen+eigenShift) !det(H) + hessianDiagonal(1)=functionHessian(1,1)+eigenShift + hessianDiagonal(2)=functionHessian(2,2)+eigenShift ENDIF - ENDDO - CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector(1:numberOfCoordinates)=dataPointLocation-interpolatedPoint%values(:,NO_PART_DERIV) - newFunctionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - !Second half of the convergence test (before collision detection) - converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN - !Bad model: reduce step size - IF(delta<=minimumDelta) THEN - !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! - EXIT main_loop + updateXi(1)=-(hessianDiagonal(2)*functionGradient(1)-functionHessian(1,2)*functionGradient(2))/determinant + updateXi(2)=(functionHessian(1,2)*functionGradient(1)-hessianDiagonal(1)*functionGradient(2))/determinant + updateXiNorm=DSQRT(DOT_PRODUCT(updateXi,updateXi)) + free=.TRUE. + DO xiIdx=1,2 + IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN + !Projection will go out of element bounds + IF(.NOT.free) THEN !both xi are fixed + exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS + EXIT main_loop + ENDIF + free=.FALSE. + fixedXiIdx=xiIdx + ENDIF + ENDDO !xiIdx + IF(free) THEN + !Both xi are free + IF(.NOT.insideRegion) THEN + IF(updateXiNorm>0.0_DP) THEN + updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound + ENDIF + ENDIF + ELSE + !xi are not free + updateXi(fixedXiIdx)=0.0_DP + xiIdx=3-fixedXiIdx + insideRegion=.FALSE. + IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN + !Positive: minimum exists in the unbounded direction + updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) + updateXiNorm=DABS(updateXi(xiIdx)) + insideRegion=updateXiNorm<=delta + ENDIF + IF(.NOT.insideRegion) THEN + !Minimum not in the region + updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) + updateXiNorm=delta + ENDIF ENDIF - delta=DMAX1(minimumDelta,0.25_DP*delta) - ELSE - predictedReduction=DOT_PRODUCT(functionGradient,updateXi)+ & - & 0.5_DP*(updateXi(1)*(updateXi(1)*functionHessian(1,1)+2.0_DP*updateXi(2)*functionHessian(1,2))+ & - & updateXi(2)**2*functionHessian(2,2)) - predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction - IF(predictionAccuracy<0.01_DP) THEN - !Bad model: reduce region size + !First half of the convergence test + converged=updateXiNorm1.0_DP) THEN + newXi(xiIdx)=1.0_DP + newXi(3-xiIdx)=xi(3-xiIdx)+updateXi(3-xiIdx)*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) + ENDIF + ENDDO + CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation-interpolatedPoint%values(:,NO_PART_DERIV) + newFunctionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Second half of the convergence test (before collision detection) + converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN + !Bad model: reduce step size IF(delta<=minimumDelta) THEN !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.5_DP*delta) - ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN - !Good model: increase region size - delta=DMIN1(maximumDelta,2.0_DP*delta) - EXIT + delta=DMAX1(minimumDelta,0.25_DP*delta) ELSE - !OK model: keep the current region size - EXIT + predictedReduction=DOT_PRODUCT(functionGradient,updateXi)+ & + & 0.5_DP*(updateXi(1)*(updateXi(1)*functionHessian(1,1)+2.0_DP*updateXi(2)*functionHessian(1,2))+ & + & updateXi(2)**2*functionHessian(2,2)) + predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction + IF(predictionAccuracy<0.01_DP) THEN + !Bad model: reduce region size + IF(delta<=minimumDelta) THEN + !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + EXIT main_loop + ENDIF + delta=DMAX1(minimumDelta,0.5_DP*delta) + ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN + !Good model: increase region size + delta=DMIN1(maximumDelta,2.0_DP*delta) + EXIT + ELSE + !OK model: keep the current region size + EXIT + ENDIF ENDIF + ENDDO !iterationIdx2 (inner loop: adjust region size) + functionValue=newFunctionValue + xi=newXi + IF(converged) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED + EXIT ENDIF - ENDDO !iterationIdx2 (inner loop: adjust region size) - functionValue=newFunctionValue - xi=newXi - IF(converged) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED - EXIT + ENDDO main_loop !iterationIdx1 (outer loop) + IF(exitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT.AND.iterationIdx1>=dataProjection%maximumNumberOfIterations) & + & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION + IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)=dataProjection%maximumNumberOfIterations) & - & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION - IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%domain(meshComponentNumber)%ptr% & & mappings%elements numberOfCoordinates=dataProjection%numberOfCoordinates - relativeTolerance=dataProjection%relativeTolerance - absoluteTolerance=dataProjection%absoluteTolerance - maximumDelta=dataProjection%maximumIterationUpdate - minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small - !Project on each candidate elements - DO elementIdx=1,SIZE(candidateElements,1) - elementNumber=candidateElements(elementIdx) - IF(elementNumber>0) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT - converged=.FALSE. - !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet - delta=0.5_DP*maximumDelta - CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,elementNumber, & - & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - xi=dataProjection%startingXi - CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & - interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - functionValue=DOT_PRODUCT(distanceVector,distanceVector) - !Outer loop - main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations - !Check for bounds [0,1] - DO xiIdx=1,3 - IF(ABS(xi(xiIdx))-1.0E-5_DP) THEN !include some negatives for numerical errors - minEigen=-temp1 !all eigenvalues are the same - ELSE - temp3=DSQRT(-temp3) - temp4=(determinant+3.0_DP*(temp1*temp2)-2.0_DP*temp1**3)/(2.0_DP*temp3**3) - minEigen=2.0_DP*temp3*DCOS((DACOS(temp4)+TWOPI)/3.0_DP)-temp1 - ENDIF - functionGradientNorm=DSQRT(DOT_PRODUCT(functionGradient,functionGradient)) - !Inner loop: adjust region size, usually EXIT at 1 or 2 iterations - DO iterationIdx2=1,dataProjection%maximumNumberOfIterations - temp1=functionGradientNorm/delta - !Estimate if the solution is inside the trust region without calculating a Newton step, this also guarantees - !the Hessian matrix is positive definite - insideRegion=(minEigen>=temp1).AND.(minEigen>absoluteTolerance) - IF(insideRegion) THEN - hessianDiagonal(1)=functionHessian(1,1) - hessianDiagonal(2)=functionHessian(2,2) - hessianDiagonal(3)=functionHessian(3,3) - ELSE - eigenShift=MAX(temp1-minEigen,absoluteTolerance) !shift towards steepest decent - determinant=determinant+eigenShift*(trace2+eigenShift*(trace+eigenShift)) !shift the determinant - hessianDiagonal(1)=functionHessian(1,1)+eigenShift - hessianDiagonal(2)=functionHessian(2,2)+eigenShift - hessianDiagonal(3)=functionHessian(3,3)+eigenShift - ENDIF - temp2=functionHessian(1,3)*functionHessian(2,3)-functionHessian(1,2)*hessianDiagonal(3) - temp3=functionHessian(1,2)*functionHessian(2,3)-functionHessian(1,3)*hessianDiagonal(2) - temp4=functionHessian(1,2)*functionHessian(1,3)-functionHessian(2,3)*hessianDiagonal(1) - updateXi(1)=((functionHessian(2,3)**2-hessianDiagonal(2)*hessianDiagonal(3))*functionGradient(1)- & - & temp2*functionGradient(2)-temp3*functionGradient(3))/determinant - updateXi(2)=((functionHessian(1,3)**2-hessianDiagonal(1)*hessianDiagonal(3))*functionGradient(2)- & - & temp2*functionGradient(1)-temp4*functionGradient(3))/determinant - updateXi(3)=((functionHessian(1,2)**2-hessianDiagonal(1)*hessianDiagonal(2))*functionGradient(3)- & - & temp3*functionGradient(1)-temp4*functionGradient(2))/determinant - updateXiNorm=DSQRT(DOT_PRODUCT(updateXi,updateXi)) - free=.TRUE. - nBound=0 + IF(projectionExitTag==DATA_PROJECTION_USER_SPECIFIED) THEN + IF(projectionElementNumber==0) CALL FlagError("The projection element number has not been set.",err,error,*999) + CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,projectionElementNumber, & + & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + CALL Field_InterpolateXi(NO_PART_DERIV,projectionXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + projectionVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + projectionDistance=DSQRT(DOT_PRODUCT(projectionVector(1:numberOfCoordinates),projectionVector(1:numberOfCoordinates))) + ELSE + projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + relativeTolerance=dataProjection%relativeTolerance + absoluteTolerance=dataProjection%absoluteTolerance + maximumDelta=dataProjection%maximumIterationUpdate + minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small + !Project on each candidate elements + DO elementIdx=1,SIZE(candidateElements,1) + elementNumber=candidateElements(elementIdx) + IF(elementNumber>0) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + converged=.FALSE. + !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet + delta=0.5_DP*maximumDelta + CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,elementNumber, & + & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + xi=dataProjection%startingXi + CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + functionValue=DOT_PRODUCT(distanceVector,distanceVector) + !Outer loop + main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations + !Check for bounds [0,1] DO xiIdx=1,3 - IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN - !Projection will go out of element bounds - nBound=nBound+1 - free=.FALSE. - IF(nBound<=2) THEN - fixedXiIdx2(nBound)=xiIdx - ELSE - !All xi are fixed - exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS - EXIT main_loop - ENDIF - ENDIF - ENDDO !xiIdx - IF(free) THEN - !All xi are free - IF(.NOT.insideRegion) THEN - IF(updateXiNorm>0.0_DP) THEN - updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound - ENDIF + IF(ABS(xi(xiIdx))-1.0E-5_DP) THEN !include some negatives for numerical errors + minEigen=-temp1 !all eigenvalues are the same ELSE - !At least one of the xi are not free - !Try 2D projection - free=.TRUE. - fixedXiIdx=fixedXiIdx2(1) - IF(nBound==2) THEN - !only fix the direction that is most strongly suggesting leaving the element - IF(updateXi(fixedXiIdx2(2))>updateXi(fixedXiIdx2(1))) fixedXiIdx=fixedXiIdx2(2) - ENDIF - updateXi(fixedXiIdx)=0.0_DP - faceXiIdxs(1)=1+MOD(fixedXiIdx,3) - faceXiIdxs(2)=1+MOD(fixedXiIdx+1,3) - !functionGradient2 - functionGradient2(1)=functionGradient(faceXiIdxs(1)) - functionGradient2(2)=functionGradient(faceXiIdxs(2)) - !functionHessian2 - functionHessian2(1,1)=functionHessian(faceXiIdxs(1),faceXiIdxs(1)) - functionHessian2(1,2)=functionHessian(faceXiIdxs(1),faceXiIdxs(2)) - functionHessian2(2,2)=functionHessian(faceXiIdxs(2),faceXiIdxs(2)) - !re-estimate the trust solution in 2D - temp1=0.5_DP*(functionHessian2(1,1)+functionHessian2(2,2)) - temp2=DSQRT((0.5_DP*(functionHessian2(1,1)-functionHessian2(2,2)))**2+functionHessian2(1,2)**2) - minEigen=temp1-temp2 - maxEigen=temp1+temp2 - temp3=DSQRT(DOT_PRODUCT(functionGradient2,functionGradient2))/delta + temp3=DSQRT(-temp3) + temp4=(determinant+3.0_DP*(temp1*temp2)-2.0_DP*temp1**3)/(2.0_DP*temp3**3) + minEigen=2.0_DP*temp3*DCOS((DACOS(temp4)+TWOPI)/3.0_DP)-temp1 + ENDIF + functionGradientNorm=DSQRT(DOT_PRODUCT(functionGradient,functionGradient)) + !Inner loop: adjust region size, usually EXIT at 1 or 2 iterations + DO iterationIdx2=1,dataProjection%maximumNumberOfIterations + temp1=functionGradientNorm/delta !Estimate if the solution is inside the trust region without calculating a Newton step, this also guarantees !the Hessian matrix is positive definite - insideRegion=(minEigen>=temp3).AND.(minEigen>absoluteTolerance) + insideRegion=(minEigen>=temp1).AND.(minEigen>absoluteTolerance) IF(insideRegion) THEN - determinant=minEigen*maxEigen !determinant of functionHessian - hessianDiagonal2(1)=functionHessian2(1,1) - hessianDiagonal2(2)=functionHessian2(2,2) + hessianDiagonal(1)=functionHessian(1,1) + hessianDiagonal(2)=functionHessian(2,2) + hessianDiagonal(3)=functionHessian(3,3) ELSE - eigenShift=MAX(temp3-minEigen,absoluteTolerance) !shift towards steepest decent - determinant=temp3*(maxEigen+eigenShift) !determinant of shifted functionHessian - hessianDiagonal2(1)=functionHessian2(1,1)+eigenShift - hessianDiagonal2(2)=functionHessian2(2,2)+eigenShift + eigenShift=MAX(temp1-minEigen,absoluteTolerance) !shift towards steepest decent + determinant=determinant+eigenShift*(trace2+eigenShift*(trace+eigenShift)) !shift the determinant + hessianDiagonal(1)=functionHessian(1,1)+eigenShift + hessianDiagonal(2)=functionHessian(2,2)+eigenShift + hessianDiagonal(3)=functionHessian(3,3)+eigenShift ENDIF - updateXi(faceXiIdxs(1))=-(hessianDiagonal2(2)*functionGradient2(1)- & - & functionHessian2(1,2)*functionGradient2(2))/determinant - updateXi(faceXiIdxs(2))=(functionHessian2(1,2)*functionGradient2(1)- & - & hessianDiagonal2(1)*functionGradient2(2))/determinant + temp2=functionHessian(1,3)*functionHessian(2,3)-functionHessian(1,2)*hessianDiagonal(3) + temp3=functionHessian(1,2)*functionHessian(2,3)-functionHessian(1,3)*hessianDiagonal(2) + temp4=functionHessian(1,2)*functionHessian(1,3)-functionHessian(2,3)*hessianDiagonal(1) + updateXi(1)=((functionHessian(2,3)**2-hessianDiagonal(2)*hessianDiagonal(3))*functionGradient(1)- & + & temp2*functionGradient(2)-temp3*functionGradient(3))/determinant + updateXi(2)=((functionHessian(1,3)**2-hessianDiagonal(1)*hessianDiagonal(3))*functionGradient(2)- & + & temp2*functionGradient(1)-temp4*functionGradient(3))/determinant + updateXi(3)=((functionHessian(1,2)**2-hessianDiagonal(1)*hessianDiagonal(2))*functionGradient(3)- & + & temp3*functionGradient(1)-temp4*functionGradient(2))/determinant updateXiNorm=DSQRT(DOT_PRODUCT(updateXi,updateXi)) - !Check again for bounds - DO boundXiIdx=1,2 - IF((bound(faceXiIdxs(boundXiIdx))/=0).AND.(bound(faceXiIdxs(boundXiIdx))>0.EQV. & - & updateXi(faceXiIdxs(boundXiIdx))>0.0_DP)) THEN + free=.TRUE. + nBound=0 + DO xiIdx=1,3 + IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN !Projection will go out of element bounds - IF(.NOT.free) THEN - !Both xi are fixed + nBound=nBound+1 + free=.FALSE. + IF(nBound<=2) THEN + fixedXiIdx2(nBound)=xiIdx + ELSE + !All xi are fixed exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS EXIT main_loop ENDIF - free=.FALSE. - fixedBoundXiIdx=boundXiIdx ENDIF ENDDO !xiIdx IF(free) THEN - !Both xi are free + !All xi are free IF(.NOT.insideRegion) THEN IF(updateXiNorm>0.0_DP) THEN updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound ENDIF ENDIF ELSE - !xi are not free - updateXi(faceXiIdxs(fixedBoundXiIdx))=0.0_DP - xiIdx=faceXiIdxs(3-fixedBoundXiIdx) - insideRegion=.FALSE. - IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN - !positive: minimum exists in the unbounded direction - updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) - updateXiNorm=DABS(updateXi(xiIdx)) - insideRegion=updateXiNorm<=delta + !At least one of the xi are not free + !Try 2D projection + free=.TRUE. + fixedXiIdx=fixedXiIdx2(1) + IF(nBound==2) THEN + !only fix the direction that is most strongly suggesting leaving the element + IF(updateXi(fixedXiIdx2(2))>updateXi(fixedXiIdx2(1))) fixedXiIdx=fixedXiIdx2(2) ENDIF - IF(.NOT.insideRegion) THEN - !minimum not in the region - updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) - updateXiNorm=delta + updateXi(fixedXiIdx)=0.0_DP + faceXiIdxs(1)=1+MOD(fixedXiIdx,3) + faceXiIdxs(2)=1+MOD(fixedXiIdx+1,3) + !functionGradient2 + functionGradient2(1)=functionGradient(faceXiIdxs(1)) + functionGradient2(2)=functionGradient(faceXiIdxs(2)) + !functionHessian2 + functionHessian2(1,1)=functionHessian(faceXiIdxs(1),faceXiIdxs(1)) + functionHessian2(1,2)=functionHessian(faceXiIdxs(1),faceXiIdxs(2)) + functionHessian2(2,2)=functionHessian(faceXiIdxs(2),faceXiIdxs(2)) + !re-estimate the trust solution in 2D + temp1=0.5_DP*(functionHessian2(1,1)+functionHessian2(2,2)) + temp2=DSQRT((0.5_DP*(functionHessian2(1,1)-functionHessian2(2,2)))**2+functionHessian2(1,2)**2) + minEigen=temp1-temp2 + maxEigen=temp1+temp2 + temp3=DSQRT(DOT_PRODUCT(functionGradient2,functionGradient2))/delta + !Estimate if the solution is inside the trust region without calculating a Newton step, this also guarantees + !the Hessian matrix is positive definite + insideRegion=(minEigen>=temp3).AND.(minEigen>absoluteTolerance) + IF(insideRegion) THEN + determinant=minEigen*maxEigen !determinant of functionHessian + hessianDiagonal2(1)=functionHessian2(1,1) + hessianDiagonal2(2)=functionHessian2(2,2) + ELSE + eigenShift=MAX(temp3-minEigen,absoluteTolerance) !shift towards steepest decent + determinant=temp3*(maxEigen+eigenShift) !determinant of shifted functionHessian + hessianDiagonal2(1)=functionHessian2(1,1)+eigenShift + hessianDiagonal2(2)=functionHessian2(2,2)+eigenShift ENDIF - ENDIF !if xi are free (2D) - ENDIF !if xi are free (3D) - !First half of the convergence test - converged=updateXiNorm0.EQV. & + & updateXi(faceXiIdxs(boundXiIdx))>0.0_DP)) THEN + !Projection will go out of element bounds + IF(.NOT.free) THEN + !Both xi are fixed + exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS + EXIT main_loop ENDIF - ENDDO - ELSEIF(newXi(xiIdx)>1.0_DP) THEN - newXi(xiIdx)=1.0_DP - DO xiIdx2 = 1,3 - IF(xiIdx2 /= xiIdx) THEN - newXi(xiIdx2)=xi(xiIdx2)+updateXi(xiIdx2) + free=.FALSE. + fixedBoundXiIdx=boundXiIdx + ENDIF + ENDDO !xiIdx + IF(free) THEN + !Both xi are free + IF(.NOT.insideRegion) THEN + IF(updateXiNorm>0.0_DP) THEN + updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound ENDIF - ENDDO + ENDIF + ELSE + !xi are not free + updateXi(faceXiIdxs(fixedBoundXiIdx))=0.0_DP + xiIdx=faceXiIdxs(3-fixedBoundXiIdx) + insideRegion=.FALSE. + IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN + !positive: minimum exists in the unbounded direction + updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) + updateXiNorm=DABS(updateXi(xiIdx)) + insideRegion=updateXiNorm<=delta + ENDIF + IF(.NOT.insideRegion) THEN + !minimum not in the region + updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) + updateXiNorm=delta + ENDIF + ENDIF !if xi are free (2D) + ENDIF !if xi are free (3D) + !First half of the convergence test + converged=updateXiNorm1.0_DP) THEN + newXi(xiIdx)=1.0_DP + DO xiIdx2 = 1,3 + IF(xiIdx2 /= xiIdx) THEN + newXi(xiIdx2)=xi(xiIdx2)+updateXi(xiIdx2) + ENDIF + ENDDO + ENDIF + ELSE IF(newXi(xiIdx)<0.0_DP) THEN + !boundary collision check + newXi(xiIdx)=0.0_DP + newXi=xi-updateXi*xi(xiIdx)/updateXi(xiIdx) + ELSE IF(newXi(xiIdx)>1.0_DP) THEN + newXi(xiIdx)=1.0_DP + newXi=xi+updateXi*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) ENDIF - ELSE IF(newXi(xiIdx)<0.0_DP) THEN - !boundary collision check - newXi(xiIdx)=0.0_DP - newXi=xi-updateXi*xi(xiIdx)/updateXi(xiIdx) - ELSE IF(newXi(xiIdx)>1.0_DP) THEN - newXi(xiIdx)=1.0_DP - newXi=xi+updateXi*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) - ENDIF - ENDDO !xiIdx - CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector=dataPointLocation-interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - newFunctionValue=DOT_PRODUCT(distanceVector,distanceVector) - !Second half of the convergence test (before collision detection) - converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN - !bad model: reduce step size - IF(delta<=minimumDelta) THEN - !something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION ! \it will get stucked!! - EXIT main_loop - ENDIF - delta=DMAX1(minimumDelta,0.25_DP*delta) - ELSE - predictedReduction=DOT_PRODUCT(functionGradient,updateXi)+ & - & 0.5_DP*(updateXi(1)*(updateXi(1)*functionHessian(1,1)+2.0_DP*updateXi(2)*functionHessian(1,2)+ & - & 2.0_DP*updateXi(3)*functionHessian(1,3))+updateXi(2)*(updateXi(2)*functionHessian(2,2)+ & - & 2.0_DP*updateXi(3)*functionHessian(2,3))+updateXi(2)**2*functionHessian(2,2)) - IF (ABS(predictedReduction) < ZERO_TOLERANCE) THEN - converged = .TRUE. !prediction reduction converged - EXIT - ENDIF - predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction - IF(predictionAccuracy<0.01_DP) THEN - !bad model: reduce region size + ENDDO !xiIdx + CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector=dataPointLocation-interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + newFunctionValue=DOT_PRODUCT(distanceVector,distanceVector) + !Second half of the convergence test (before collision detection) + converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN + !bad model: reduce step size IF(delta<=minimumDelta) THEN !something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION ! \it will get stucked!! EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.5_DP*delta) - ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN - !good model: increase region size - delta=DMIN1(maximumDelta,2.0_DP*delta) - EXIT + delta=DMAX1(minimumDelta,0.25_DP*delta) ELSE - !ok model: keep the current region size - EXIT + predictedReduction=DOT_PRODUCT(functionGradient,updateXi)+ & + & 0.5_DP*(updateXi(1)*(updateXi(1)*functionHessian(1,1)+2.0_DP*updateXi(2)*functionHessian(1,2)+ & + & 2.0_DP*updateXi(3)*functionHessian(1,3))+updateXi(2)*(updateXi(2)*functionHessian(2,2)+ & + & 2.0_DP*updateXi(3)*functionHessian(2,3))+updateXi(2)**2*functionHessian(2,2)) + IF (ABS(predictedReduction) < ZERO_TOLERANCE) THEN + converged = .TRUE. !prediction reduction converged + EXIT + ENDIF + predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction + IF(predictionAccuracy<0.01_DP) THEN + !bad model: reduce region size + IF(delta<=minimumDelta) THEN + !something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + EXIT main_loop + ENDIF + delta=DMAX1(minimumDelta,0.5_DP*delta) + ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN + !good model: increase region size + delta=DMIN1(maximumDelta,2.0_DP*delta) + EXIT + ELSE + !ok model: keep the current region size + EXIT + ENDIF ENDIF + ENDDO !iterationIdx2 (inner loop: adjust region size) + functionValue=newFunctionValue + xi=newXi + IF(converged) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED + EXIT ENDIF - ENDDO !iterationIdx2 (inner loop: adjust region size) - functionValue=newFunctionValue - xi=newXi - IF(converged) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED - EXIT + ENDDO main_loop !iterationIdx1 (outer loop) + IF(exitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT.AND.iterationIdx1>=dataProjection%maximumNumberOfIterations) & + & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION + IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)=dataProjection%maximumNumberOfIterations) & - & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION - IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)interpolatedPoint%INTERPOLATION_PARAMETERS%FIELD%decomposition%domain(meshComponentNumber)% & & PTR%MAPPINGS%ELEMENTS numberOfCoordinates=dataProjection%numberOfCoordinates - relativeTolerance=dataProjection%relativeTolerance - absoluteTolerance=dataProjection%absoluteTolerance - maximumDelta=dataProjection%maximumIterationUpdate - minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small - !Project on each candidate elements - DO elementIdx=1,SIZE(candidateElements,1) - elementNumber=candidateElements(elementIdx) - IF(elementNumber>0) THEN - elementFaceNumber=candidateElementFaces(elementIdx) - faceNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & - & elements(elementNumber)%ELEMENT_FACES(elementFaceNumber) - exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT - converged=.FALSE. - !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet - delta=0.5_DP*maximumDelta - CALL Field_InterpolationParametersFaceGet(dataProjection%projectionSetType,faceNumber,interpolatedPoint% & - & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - xi=dataProjection%startingXi - CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & - & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - !Outer loop - main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations - !Check for bounds [0,1] - DO xiIdx=1,2 - IF(ABS(xi(xiIdx))=temp1).AND.(minEigen>absoluteTolerance) - IF(insideRegion) THEN - determinant=minEigen*maxEigen !det(H) - hessianDiagonal(1)=functionHessian(1,1) - hessianDiagonal(2)=functionHessian(2,2) - ELSE - eigenShift=MAX(temp1-minEigen,absoluteTolerance) !shift towards steepest decent - determinant=temp1*(maxEigen+eigenShift) !det(H) - hessianDiagonal(1)=functionHessian(1,1)+eigenShift - hessianDiagonal(2)=functionHessian(2,2)+eigenShift - ENDIF - updateXi(1)=-(hessianDiagonal(2)*functionGradient(1)-functionHessian(1,2)*functionGradient(2))/determinant - updateXi(2)=(functionHessian(1,2)*functionGradient(1)-hessianDiagonal(1)*functionGradient(2))/determinant - updateXiNorm=DSQRT(DOT_PRODUCT(updateXi,updateXi)) - free=.TRUE. - DO xiIdx=1,2 - IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN - !Projection will go out of element bounds - IF(.NOT.free) THEN !both xi are fixed - exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS - EXIT main_loop - ENDIF - free=.FALSE. - fixedXiIdx=xiIdx - ENDIF - ENDDO !xiIdx - IF(free) THEN - !Both xi are free - IF(.NOT.insideRegion) THEN - IF(updateXiNorm>0.0_DP) THEN - updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound - ENDIF - ENDIF - ELSE - !xi are not free - updateXi(fixedXiIdx)=0.0_DP - xiIdx=3-fixedXiIdx - insideRegion=.FALSE. - IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN - !Positive: minimum exists in the unbounded direction - updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) - updateXiNorm=DABS(updateXi(xiIdx)) - insideRegion=updateXiNorm<=delta - ENDIF - IF(.NOT.insideRegion) THEN - !Minimum not in the region - updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) - updateXiNorm=delta + IF(projectionExitTag==DATA_PROJECTION_USER_SPECIFIED) THEN + IF(projectionElementNumber==0) CALL FlagError("The projection element number has not been set.",err,error,*999) + IF(projectionElementFaceNumber==0) CALL FlagError("The projection element face number has not been set.",err,error,*999) + faceNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & + & elements(projectionElementNumber)%ELEMENT_FACES(projectionElementFaceNumber) + CALL Field_InterpolationParametersFaceGet(dataProjection%projectionSetType,faceNumber,interpolatedPoint% & + & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + CALL Field_InterpolateXi(NO_PART_DERIV,projectionXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + projectionVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + projectionDistance=DSQRT(DOT_PRODUCT(projectionVector(1:numberOfCoordinates),projectionVector(1:numberOfCoordinates))) + ELSE + projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + relativeTolerance=dataProjection%relativeTolerance + absoluteTolerance=dataProjection%absoluteTolerance + maximumDelta=dataProjection%maximumIterationUpdate + minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small + !Project on each candidate elements + DO elementIdx=1,SIZE(candidateElements,1) + elementNumber=candidateElements(elementIdx) + IF(elementNumber>0) THEN + elementFaceNumber=candidateElementFaces(elementIdx) + faceNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & + & elements(elementNumber)%ELEMENT_FACES(elementFaceNumber) + exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + converged=.FALSE. + !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet + delta=0.5_DP*maximumDelta + CALL Field_InterpolationParametersFaceGet(dataProjection%projectionSetType,faceNumber,interpolatedPoint% & + & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + xi=dataProjection%startingXi + CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Outer loop + main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations + !Check for bounds [0,1] + DO xiIdx=1,2 + IF(ABS(xi(xiIdx))1.0_DP) THEN - newXi(xiIdx)=1.0_DP - newXi(3-xiIdx)=xi(3-xiIdx)+updateXi(3-xiIdx)*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) + ENDDO !xiIdx + !functionGradient + functionGradient(1)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1))) + functionGradient(2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + !functionHessian + functionHessian(1,1)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1_S1))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1))) + functionHessian(1,2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1_S2))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + functionHessian(2,2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2_S2))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + !A model trust region approach, a Newton step is taken if the minimum lies inside the trust region (delta), + !if not, shift the step towards the steepest descent + temp1=0.5_DP*(functionHessian(1,1)+functionHessian(2,2)) + temp2=DSQRT((0.5_DP*(functionHessian(1,1)-functionHessian(2,2)))**2+functionHessian(1,2)**2) + minEigen=temp1-temp2 + maxEigen=temp1+temp2 + functionGradientNorm=DSQRT(DOT_PRODUCT(functionGradient,functionGradient)) + !Inner loop: adjust region size, usually EXIT at 1 or 2 iterations + DO iterationIdx2=1,dataProjection%maximumNumberOfIterations + temp1=functionGradientNorm/delta + !Estimate if the solution is inside the trust region without calculating a Newton step, this also guarantees + !the Hessian matrix is positive definite + insideRegion=(minEigen>=temp1).AND.(minEigen>absoluteTolerance) + IF(insideRegion) THEN + determinant=minEigen*maxEigen !det(H) + hessianDiagonal(1)=functionHessian(1,1) + hessianDiagonal(2)=functionHessian(2,2) + ELSE + eigenShift=MAX(temp1-minEigen,absoluteTolerance) !shift towards steepest decent + determinant=temp1*(maxEigen+eigenShift) !det(H) + hessianDiagonal(1)=functionHessian(1,1)+eigenShift + hessianDiagonal(2)=functionHessian(2,2)+eigenShift ENDIF - ENDDO - CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector=dataPointLocation-interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - newFunctionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - !Second half of the convergence test (before collision detection) - converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN - !Bad model: reduce step size - IF(delta<=minimumDelta) THEN - !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + updateXi(1)=-(hessianDiagonal(2)*functionGradient(1)-functionHessian(1,2)*functionGradient(2))/determinant + updateXi(2)=(functionHessian(1,2)*functionGradient(1)-hessianDiagonal(1)*functionGradient(2))/determinant + updateXiNorm=DSQRT(DOT_PRODUCT(updateXi,updateXi)) + free=.TRUE. + DO xiIdx=1,2 + IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN + !Projection will go out of element bounds + IF(.NOT.free) THEN !both xi are fixed + exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS + EXIT main_loop + ENDIF + free=.FALSE. + fixedXiIdx=xiIdx + ENDIF + ENDDO !xiIdx + IF(free) THEN + !Both xi are free + IF(.NOT.insideRegion) THEN + IF(updateXiNorm>0.0_DP) THEN + updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound + ENDIF + ENDIF + ELSE + !xi are not free + updateXi(fixedXiIdx)=0.0_DP + xiIdx=3-fixedXiIdx + insideRegion=.FALSE. + IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN + !Positive: minimum exists in the unbounded direction + updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) + updateXiNorm=DABS(updateXi(xiIdx)) + insideRegion=updateXiNorm<=delta + ENDIF + IF(.NOT.insideRegion) THEN + !Minimum not in the region + updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) + updateXiNorm=delta + ENDIF + ENDIF !if xi is free + !First half of the convergence test + converged=updateXiNorm1.0_DP) THEN + newXi(xiIdx)=1.0_DP + newXi(3-xiIdx)=xi(3-xiIdx)+updateXi(3-xiIdx)*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) + ENDIF + ENDDO + CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector=dataPointLocation-interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + newFunctionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Second half of the convergence test (before collision detection) + converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN + !Bad model: reduce step size IF(delta<=minimumDelta) THEN !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! functionValue=newFunctionValue EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.5_DP*delta) - ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN - !Good model: increase region size - delta=DMIN1(maximumDelta,2.0_DP*delta) - functionValue=newFunctionValue - EXIT + delta=DMAX1(minimumDelta,0.25_DP*delta) ELSE - !OK model: keep the current region size - functionValue=newFunctionValue - EXIT + predictedReduction=DOT_PRODUCT(functionGradient,updateXi)+ & + & 0.5_DP*(updateXi(1)*(updateXi(1)*functionHessian(1,1)+2.0_DP*updateXi(2)*functionHessian(1,2))+ & + & updateXi(2)**2*functionHessian(2,2)) + predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction + IF(predictionAccuracy<0.01_DP) THEN + !Bad model: reduce region size + IF(delta<=minimumDelta) THEN + !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + functionValue=newFunctionValue + EXIT main_loop + ENDIF + delta=DMAX1(minimumDelta,0.5_DP*delta) + ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN + !Good model: increase region size + delta=DMIN1(maximumDelta,2.0_DP*delta) + functionValue=newFunctionValue + EXIT + ELSE + !OK model: keep the current region size + functionValue=newFunctionValue + EXIT + ENDIF ENDIF + ENDDO !iterationIdx2 (inner loop: adjust region size) + functionValue=newFunctionValue + xi=newXi + IF(converged) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED + EXIT ENDIF - ENDDO !iterationIdx2 (inner loop: adjust region size) - functionValue=newFunctionValue - xi=newXi - IF(converged) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED - EXIT + ENDDO main_loop !iterationIdx1 (outer loop) + IF(exitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT.AND.iterationIdx1>=dataProjection%maximumNumberOfIterations) & + & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION + IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)=dataProjection%maximumNumberOfIterations) & - & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION - IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)0) THEN - elementLineNumber=candidateElementLines(elementIdx) - lineNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & - & elements(elementNumber)%ELEMENT_LINES(elementLineNumber) - exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT - converged=.FALSE. - !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet - delta=0.5_DP*maximumDelta - CALL Field_InterpolationParametersLineGet(dataProjection%projectionSetType,lineNumber,interpolatedPoint% & - & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - xi=dataProjection%startingXi - CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & - & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - !Outer loop - main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations - !Check for bounds [0,1] - IF(ABS(xi(1))absoluteTolerance) THEN - !Positive: minimum exists - updateXi(1)=-functionGradient/functionHessian - updateXiNorm=DABS(updateXi(1)) - insideRegion=updateXiNorm<=delta - ENDIF - IF(.NOT.insideRegion) THEN - !minimum not in the region - updateXi(1)=-DSIGN(delta,functionGradient) - updateXiNorm=delta - ENDIF - IF((bound/=0).AND.(bound>0.EQV.updateXi(1)>0.0_DP)) THEN - !Projection will go out of element bounds - exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS - EXIT main_loop - ENDIF - !First half of the convergence test (before collision detection) - converged=updateXiNorm1.0_DP) THEN - newXi(1)=1.0_DP + IF(projectionExitTag==DATA_PROJECTION_USER_SPECIFIED) THEN + IF(projectionElementNumber==0) CALL FlagError("The projection element number has not been set.",err,error,*999) + IF(projectionElementLineNumber==0) CALL FlagError("The projection element line number has not been set.",err,error,*999) + lineNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & + & elements(projectionElementNumber)%ELEMENT_LINES(projectionElementLineNumber) + CALL Field_InterpolationParametersLineGet(dataProjection%projectionSetType,lineNumber,interpolatedPoint% & + & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + CALL Field_InterpolateXi(NO_PART_DERIV,projectionXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + projectionVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + projectionDistance=DSQRT(DOT_PRODUCT(projectionVector(1:numberOfCoordinates),projectionVector(1:numberOfCoordinates))) + ELSE + projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + relativeTolerance=dataProjection%relativeTolerance + absoluteTolerance=dataProjection%absoluteTolerance + maximumDelta=dataProjection%maximumIterationUpdate + minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small + !Project on each candidate elements + DO elementIdx=1,SIZE(candidateElements,1) + elementNumber=candidateElements(elementIdx) + IF(elementNumber>0) THEN + elementLineNumber=candidateElementLines(elementIdx) + lineNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & + & elements(elementNumber)%ELEMENT_LINES(elementLineNumber) + exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + converged=.FALSE. + !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet + delta=0.5_DP*maximumDelta + CALL Field_InterpolationParametersLineGet(dataProjection%projectionSetType,lineNumber,interpolatedPoint% & + & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + xi=dataProjection%startingXi + CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Outer loop + main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations + !Check for bounds [0,1] + IF(ABS(xi(1))absoluteTolerance) THEN - !Bad model: reduce step size - IF(delta<=minimumDelta) THEN - !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + !functionGradient + functionGradient=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,FIRST_PART_DERIV))) + !functionHessian + functionHessian=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,SECOND_PART_DERIV))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,FIRST_PART_DERIV), & + & interpolatedPoint%values(1:numberOfCoordinates,FIRST_PART_DERIV))) + !A model trust region approach, directly taken from CMISS CLOS22: V = -(H + eigenShift*I)g + !The calculation of eigenShift are only approximated as oppose to the common trust region approach + !Inner loop: adjust region size, usually EXIT at 1 or 2 iterations + DO iterationIdx2=1,dataProjection%maximumNumberOfIterations + insideRegion=.FALSE. + IF(functionHessian>absoluteTolerance) THEN + !Positive: minimum exists + updateXi(1)=-functionGradient/functionHessian + updateXiNorm=DABS(updateXi(1)) + insideRegion=updateXiNorm<=delta + ENDIF + IF(.NOT.insideRegion) THEN + !minimum not in the region + updateXi(1)=-DSIGN(delta,functionGradient) + updateXiNorm=delta + ENDIF + IF((bound/=0).AND.(bound>0.EQV.updateXi(1)>0.0_DP)) THEN + !Projection will go out of element bounds + exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.25_DP*delta) - ELSE - predictedReduction=updateXi(1)*(functionGradient+0.5_DP*functionHessian*updateXi(1)) - predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction - IF(predictionAccuracy<0.01_DP) THEN - !Bad model: reduce region size + !First half of the convergence test (before collision detection) + converged=updateXiNorm1.0_DP) THEN + newXi(1)=1.0_DP + ENDIF + CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + newFunctionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Second half of the convergence test + converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN + !Bad model: reduce step size IF(delta<=minimumDelta) THEN !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.5_DP*delta) - ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN - !Good model: increase region size - delta=DMIN1(maximumDelta,2.0_DP*delta) - EXIT + delta=DMAX1(minimumDelta,0.25_DP*delta) ELSE - !OK model: keep the current region size - EXIT + predictedReduction=updateXi(1)*(functionGradient+0.5_DP*functionHessian*updateXi(1)) + predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction + IF(predictionAccuracy<0.01_DP) THEN + !Bad model: reduce region size + IF(delta<=minimumDelta) THEN + !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + EXIT main_loop + ENDIF + delta=DMAX1(minimumDelta,0.5_DP*delta) + ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN + !Good model: increase region size + delta=DMIN1(maximumDelta,2.0_DP*delta) + EXIT + ELSE + !OK model: keep the current region size + EXIT + ENDIF ENDIF + ENDDO !iterationIdx2 (inner loop: adjust region size) + functionValue=newFunctionValue + xi=newXi + IF(converged) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED + EXIT ENDIF - ENDDO !iterationIdx2 (inner loop: adjust region size) - functionValue=newFunctionValue - xi=newXi - IF(converged) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED - EXIT + ENDDO main_loop !iterationIdx1 (outer loop) + IF(exitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT.AND.iterationIdx1>=dataProjection%maximumNumberOfIterations) & + & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION + IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)=dataProjection%maximumNumberOfIterations) & - & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION - IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)distance) THEN - dataProjection%dataProjectionResults(dataPointIdx)%exitTag=DATA_PROJECTION_CANCELLED - dataProjection%dataProjectionResults(dataPointIdx)%elementNumber=0 - dataProjection%dataProjectionResults(dataPointIdx)%elementLineFaceNumber=0 - dataProjection%dataProjectionResults(dataPointIdx)%xi=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%elementXi=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%distance=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%projectionVector=0.0_DP + CALL DataProjection_DataProjectionResultCancel(dataProjection%dataProjectionResults(dataPointIdx),err,error,*999) ENDIF ENDDO !dataPointIdx CASE(DATA_PROJECTION_DISTANCE_GREATER_EQUAL) DO dataPointIdx=1,dataPoints%numberOfDataPoints IF(dataProjection%dataProjectionResults(dataPointIdx)%distance>=distance) THEN - dataProjection%dataProjectionResults(dataPointIdx)%exitTag=DATA_PROJECTION_CANCELLED - dataProjection%dataProjectionResults(dataPointIdx)%elementNumber=0 - dataProjection%dataProjectionResults(dataPointIdx)%elementLineFaceNumber=0 - dataProjection%dataProjectionResults(dataPointIdx)%xi=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%elementXi=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%distance=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%projectionVector=0.0_DP + CALL DataProjection_DataProjectionResultCancel(dataProjection%dataProjectionResults(dataPointIdx),err,error,*999) ENDIF ENDDO !dataPointIdx CASE(DATA_PROJECTION_DISTANCE_LESS) DO dataPointIdx=1,dataPoints%numberOfDataPoints IF(dataProjection%dataProjectionResults(dataPointIdx)%distanceSets the projection element number for a data point identified by a given global number. + SUBROUTINE DataProjection_ResultElementNumberSet(dataProjection,dataPointUserNumber,projectionElementUserNumber,err,error,*) + + !Argument variables + TYPE(DataProjectionType), POINTER :: dataProjection !Gets the projection element face number for a data point identified by a given global number. SUBROUTINE DataProjection_ResultElementFaceNumberGet(dataProjection,dataPointUserNumber,projectionElementFaceNumber & & ,err,error,*) @@ -4940,6 +5076,86 @@ END SUBROUTINE DataProjection_ResultElementFaceNumberGet !================================================================================================================================ ! + !>Sets the projection element face number for a data point identified by a given global number. + SUBROUTINE DataProjection_ResultElementFaceNumberSet(dataProjection,dataPointUserNumber,elementUserNumber,localFaceNormal, & + & err,error,*) + + !Argument variables + TYPE(DataProjectionType), POINTER :: dataProjection !Gets the projection element line number for a data point identified by a given global number. SUBROUTINE DataProjection_ResultElementLineNumberGet(dataProjection,dataPointUserNumber,projectionElementLineNumber & & ,err,error,*) @@ -4978,6 +5194,91 @@ END SUBROUTINE DataProjection_ResultElementLineNumberGet !================================================================================================================================ ! + !>Sets the projection element line number for a data point identified by a given global number. + SUBROUTINE DataProjection_ResultElementLineNumberSet(dataProjection,dataPointUserNumber,elementUserNumber,localLineNormals, & + & err,error,*) + + !Argument variables + TYPE(DataProjectionType), POINTER :: dataProjection !Gets the projection exit tag for a data point identified by a given global number. SUBROUTINE DataProjection_ResultExitTagGet(dataProjection,dataPointUserNumber,projectionExitTag,err,error,*) @@ -5066,7 +5367,7 @@ SUBROUTINE DataProjection_ResultXiSet(dataProjection,dataPointUserNumber,project IF(.NOT.ASSOCIATED(dataProjection)) CALL FlagError("Data projection is not associated.",err,error,*999) IF(.NOT.dataProjection%dataProjectionFinished) CALL FlagError("Data projection has not been finished.",err,error,*999) - IF(.NOT.dataProjection%dataProjectionProjected) CALL FlagError("Data projection has not been projected.",err,error,*999) + IF(dataProjection%dataProjectionProjected) CALL FlagError("Data projection has already been projected.",err,error,*999) IF(SIZE(projectionXi,1)/=dataProjection%numberOfXi) THEN localError="The specified projection xi has size of "//TRIM(NumberToVString(SIZE(projectionXi,1),"*",err,error))// & & "but it needs to have size of "//TRIM(NumberToVString(dataProjection%numberOfXi,"*",err,error))//"." diff --git a/src/interface_routines.F90 b/src/interface_routines.F90 index 426c36cc..94e5a7d1 100644 --- a/src/interface_routines.F90 +++ b/src/interface_routines.F90 @@ -2194,7 +2194,7 @@ SUBROUTINE InterfacePointsConnectivity_UpdateFromProjection(InterfacePointsConne InterfacePointsConnectivity%pointsConnectivity(dataPointIdx,coupledMeshIndex)%reducedXi(:)= & & dataProjectionResult%xi InterfacePointsConnectivity%pointsConnectivity(dataPointIdx,coupledMeshIndex)%coupledMeshElementNumber= & - & dataProjectionResult%elementNumber + & dataProjectionResult%elementLocalNumber InterfacePointsConnectivity%pointsConnectivity(dataPointIdx,coupledMeshIndex)%elementLineFaceNumber= & & dataProjectionResult%elementLineFaceNumber ENDDO !dataPointIdx diff --git a/src/mesh_routines.F90 b/src/mesh_routines.F90 index 15fbfe48..34232d65 100644 --- a/src/mesh_routines.F90 +++ b/src/mesh_routines.F90 @@ -8567,7 +8567,7 @@ SUBROUTINE MeshTopology_DataPointsCalculateProjection(mesh,dataProjection,err,er !Calculate number of projected data points on an element DO dataPointIdx=1,dataPoints%numberOfDataPoints dataProjectionResult=>dataProjection%dataProjectionResults(dataPointIdx) - elementNumber=dataProjectionResult%elementNumber + elementNumber=dataProjectionResult%elementGlobalNumber exitTag=dataProjectionResult%exitTag IF(exitTag/=DATA_PROJECTION_CANCELLED) THEN DO elementIdx=1,elements%NUMBER_OF_ELEMENTS @@ -8593,7 +8593,7 @@ SUBROUTINE MeshTopology_DataPointsCalculateProjection(mesh,dataProjection,err,er dataPointsTopology%totalNumberOfProjectedData=0 DO dataPointIdx=1,dataPoints%numberOfDataPoints dataProjectionResult=>dataProjection%dataProjectionResults(dataPointIdx) - elementNumber=dataProjectionResult%elementNumber + elementNumber=dataProjectionResult%elementGlobalNumber exitTag=dataProjectionResult%exitTag IF(exitTag/=DATA_PROJECTION_CANCELLED) THEN DO elementIdx=1,elements%NUMBER_OF_ELEMENTS diff --git a/src/opencmiss_iron.F90 b/src/opencmiss_iron.F90 index 300cf38f..62ec0edf 100644 --- a/src/opencmiss_iron.F90 +++ b/src/opencmiss_iron.F90 @@ -1840,6 +1840,7 @@ MODULE OpenCMISS_Iron !> \see OpenCMISS::Iron::DataProjection,OpenCMISS !>@{ INTEGER(INTG), PARAMETER :: CMFE_DATA_PROJECTION_CANCELLED = DATA_PROJECTION_CANCELLED !Returns the projection element number for a data point identified by a given user number. + !>Returns the projection element number for a data point. INTERFACE cmfe_DataProjection_ResultElementNumberGet MODULE PROCEDURE cmfe_DataProjection_ResultElementNumberGetNumber MODULE PROCEDURE cmfe_DataProjection_ResultElementNumberGetObj END INTERFACE cmfe_DataProjection_ResultElementNumberGet - !>Returns the projection element face number for a data point identified by a given user number. + !>Sets the projection element number for a data point. + INTERFACE cmfe_DataProjection_ResultElementNumberSet + MODULE PROCEDURE cmfe_DataProjection_ResultElementNumberSetNumber + MODULE PROCEDURE cmfe_DataProjection_ResultElementNumberSetObj + END INTERFACE cmfe_DataProjection_ResultElementNumberSet + + !>Returns the projection element face number for a data point. INTERFACE cmfe_DataProjection_ResultElementFaceNumberGet MODULE PROCEDURE cmfe_DataProjection_ResultElementFaceNumberGetNumber MODULE PROCEDURE cmfe_DataProjection_ResultElementFaceNumberGetObj END INTERFACE cmfe_DataProjection_ResultElementFaceNumberGet - !>Returns the projection element line number for a data point identified by a given user number. + !>Sets the projection element face number for a data point. + INTERFACE cmfe_DataProjection_ResultElementFaceNumberSet + MODULE PROCEDURE cmfe_DataProjection_ResultElementFaceNumberSetNumber + MODULE PROCEDURE cmfe_DataProjection_ResultElementFaceNumberSetObj + END INTERFACE cmfe_DataProjection_ResultElementFaceNumberSet + + !>Returns the projection element line number for a data point. INTERFACE cmfe_DataProjection_ResultElementLineNumberGet MODULE PROCEDURE cmfe_DataProjection_ResultElementLineNumberGetNumber MODULE PROCEDURE cmfe_DataProjection_ResultElementLineNumberGetObj END INTERFACE cmfe_DataProjection_ResultElementLineNumberGet - !>Returns the projection exit tag for a data point identified by a given user number. + !>Sets the projection element line number for a data point. + INTERFACE cmfe_DataProjection_ResultElementLineNumberSet + MODULE PROCEDURE cmfe_DataProjection_ResultElementLineNumberSetNumber + MODULE PROCEDURE cmfe_DataProjection_ResultElementLineNumberSetObj + END INTERFACE cmfe_DataProjection_ResultElementLineNumberSet + + !>Returns the projection exit tag for a data point. INTERFACE cmfe_DataProjection_ResultExitTagGet MODULE PROCEDURE cmfe_DataProjection_ResultExitTagGetNumber MODULE PROCEDURE cmfe_DataProjection_ResultExitTagGetObj @@ -2175,8 +2194,8 @@ MODULE OpenCMISS_Iron PUBLIC CMFE_DATA_PROJECTION_BOUNDARY_LINES_PROJECTION_TYPE,CMFE_DATA_PROJECTION_BOUNDARY_FACES_PROJECTION_TYPE, & & CMFE_DATA_PROJECTION_ALL_ELEMENTS_PROJECTION_TYPE - PUBLIC CMFE_DATA_PROJECTION_CANCELLED,CMFE_DATA_PROJECTION_EXIT_TAG_CONVERGED,CMFE_DATA_PROJECTION_EXIT_TAG_BOUNDS, & - & CMFE_DATA_PROJECTION_EXIT_TAG_MAX_ITERATION,CMFE_DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + PUBLIC CMFE_DATA_PROJECTION_CANCELLED,CMFE_DATA_PROJECTION_USER_SPECIFIED,CMFE_DATA_PROJECTION_EXIT_TAG_CONVERGED, & + & CMFE_DATA_PROJECTION_EXIT_TAG_BOUNDS,CMFE_DATA_PROJECTION_EXIT_TAG_MAX_ITERATION,CMFE_DATA_PROJECTION_EXIT_TAG_NO_ELEMENT PUBLIC CMFE_DATA_PROJECTION_DISTANCE_GREATER,CMFE_DATA_PROJECTION_DISTANCE_GREATER_EQUAL,CMFE_DATA_PROJECTION_DISTANCE_LESS, & & CMFE_DATA_PROJECTION_DISTANCE_LESS_EQUAL @@ -2227,9 +2246,13 @@ MODULE OpenCMISS_Iron PUBLIC cmfe_DataProjection_ResultAnalysisOutput - PUBLIC cmfe_DataProjection_ResultDistanceGet,cmfe_DataProjection_ResultElementNumberGet + PUBLIC cmfe_DataProjection_ResultDistanceGet + + PUBLIC cmfe_DataProjection_ResultElementNumberGet,cmfe_DataProjection_ResultElementNumberSet + + PUBLIC cmfe_DataProjection_ResultElementFaceNumberGet,cmfe_DataProjection_ResultElementFaceNumberSet - PUBLIC cmfe_DataProjection_ResultElementFaceNumberGet,cmfe_DataProjection_ResultElementLineNumberGet + PUBLIC cmfe_DataProjection_ResultElementLineNumberGet,cmfe_DataProjection_ResultElementLineNumberSet PUBLIC cmfe_DataProjection_ResultExitTagGet @@ -22794,6 +22817,75 @@ END SUBROUTINE cmfe_DataProjection_ResultElementNumberGetObj !================================================================================================================================ ! + !>Sets the projection element number for a data point in a set of data points identified by user number. + SUBROUTINE cmfe_DataProjection_ResultElementNumberSetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & + & dataPointUserNumber,projectionElementUserNumber,err) + !DLLEXPORT(cmfe_DataProjection_ResultElementNumberSetNumber) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: regionUserNumber !Sets the projection element number for a data point in a set of data points identified by an object. + SUBROUTINE cmfe_DataProjection_ResultElementNumberSetObj(dataProjection,dataPointUserNumber,projectionElementUserNumber,err) + !DLLEXPORT(cmfe_DataProjection_ResultElementNumberSetObj) + + !Argument variables + TYPE(cmfe_DataProjectionType), INTENT(IN) :: dataProjection !Returns the projection element face number for a data point in a set of data points identified by user number. SUBROUTINE cmfe_DataProjection_ResultElementFaceNumberGetNumber(regionUserNumber,dataPointsUserNumber, & & dataProjectionUserNumber,dataPointUserNumber,projectionElementFaceNumber,err) @@ -22865,6 +22957,79 @@ END SUBROUTINE cmfe_DataProjection_ResultElementFaceNumberGetObj !================================================================================================================================ ! + !>Sets the projection element face number for a data point in a set of data points identified by user number. + SUBROUTINE cmfe_DataProjection_ResultElementFaceNumberSetNumber(regionUserNumber,dataPointsUserNumber, & + & dataProjectionUserNumber,dataPointUserNumber,projectionElementUserNumber,localFaceNormal,err) + !DLLEXPORT(cmfe_DataProjection_ResultElementFaceNumberSetNumber) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: regionUserNumber !Returns the projection element face number for a data point in a set of data points identified by an object. + SUBROUTINE cmfe_DataProjection_ResultElementFaceNumberSetObj(dataProjection,dataPointUserNumber, & + & projectionElementUserNumber,localFaceNormal,err) + !DLLEXPORT(cmfe_DataProjection_ResultElementFaceNumberSetObj) + + !Argument variables + TYPE(cmfe_DataProjectionType), INTENT(IN) :: dataProjection !Returns the projection element line number for a data point in a set of data points identified by user number. SUBROUTINE cmfe_DataProjection_ResultElementLineNumberGetNumber(regionUserNumber,dataPointsUserNumber, & & dataProjectionUserNumber,dataPointUserNumber,projectionElementLineNumber,err) @@ -22936,6 +23101,79 @@ END SUBROUTINE cmfe_DataProjection_ResultElementLineNumberGetObj !================================================================================================================================ ! + !>Sets the projection element line number for a data point in a set of data points identified by user number. + SUBROUTINE cmfe_DataProjection_ResultElementLineNumberSetNumber(regionUserNumber,dataPointsUserNumber, & + & dataProjectionUserNumber,dataPointUserNumber,projectionElementUserNumber,localLineNormals,err) + !DLLEXPORT(cmfe_DataProjection_ResultElementLineNumberSetNumber) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: regionUserNumber !Sets the projection element line number for a data point in a set of data points identified by an object. + SUBROUTINE cmfe_DataProjection_ResultElementLineNumberSetObj(dataProjection,dataPointUserNumber, & + & projectionElementUserNumber,localLineNormals,err) + !DLLEXPORT(cmfe_DataProjection_ResultElementLineNumberSetObj) + + !Argument variables + TYPE(cmfe_DataProjectionType), INTENT(IN) :: dataProjection !Returns the projection exit tag for a data point in a set of data points identified by user number. SUBROUTINE cmfe_DataProjection_ResultExitTagGetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & & dataPointUserNumber,projectionExitTag,err) diff --git a/src/types.F90 b/src/types.F90 index f31c42b5..ec30d228 100644 --- a/src/types.F90 +++ b/src/types.F90 @@ -288,7 +288,9 @@ MODULE Types TYPE DataProjectionResultType INTEGER(INTG) :: userNumber ! Date: Sat, 18 May 2019 15:37:29 +1200 Subject: [PATCH 02/17] Fix for no error when using wrong projection result routine. --- src/data_projection_routines.F90 | 47 ++++++++++++++++++++------------ 1 file changed, 30 insertions(+), 17 deletions(-) diff --git a/src/data_projection_routines.F90 b/src/data_projection_routines.F90 index b3082bc4..8c08d12b 100644 --- a/src/data_projection_routines.F90 +++ b/src/data_projection_routines.F90 @@ -5009,24 +5009,37 @@ SUBROUTINE DataProjection_ResultElementNumberSet(dataProjection,dataPointUserNum IF(.NOT.dataProjection%dataProjectionFinished) CALL FlagError("Data projection has not been finished.",err,error,*999) IF(dataProjection%dataProjectionProjected) CALL FlagError("Data projection has already been projected.",err,error,*999) - CALL DataProjection_DataPointGlobalNumberGet(dataProjection,dataPointUserNumber,dataPointGlobalNumber,err,error,*999) - NULLIFY(decomposition) - CALL DataProjection_DecompositionGet(dataProjection,decomposition,err,error,*999) - NULLIFY(decompositionTopology) - CALL Decomposition_TopologyGet(decomposition,decompositionTopology,err,error,*999) - CALL DecompositionTopology_ElementCheckExists(decompositionTopology,projectionElementUserNumber,elementExists, & - & elementLocalNumber,ghostElement,err,error,*999) - IF(elementExists) THEN - IF(.NOT.ghostElement) THEN - dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementLocalNumber=elementLocalNumber - dataProjection%dataProjectionResults(dataPointGlobalNumber)%exitTag=DATA_PROJECTION_USER_SPECIFIED + SELECT CASE(dataProjection%projectionType) + CASE(DATA_PROJECTION_BOUNDARY_FACES_PROJECTION_TYPE) + CALL FlagError("Use the result element faces routine not the element routine for a boundary faces projection.", & + & err,error,*999) + CASE(DATA_PROJECTION_BOUNDARY_LINES_PROJECTION_TYPE) + CALL FlagError("Use the result element lines routine not the element routine for a boundary lines projection.", & + & err,error,*999) + CASE(DATA_PROJECTION_ALL_ELEMENTS_PROJECTION_TYPE) + CALL DataProjection_DataPointGlobalNumberGet(dataProjection,dataPointUserNumber,dataPointGlobalNumber,err,error,*999) + NULLIFY(decomposition) + CALL DataProjection_DecompositionGet(dataProjection,decomposition,err,error,*999) + NULLIFY(decompositionTopology) + CALL Decomposition_TopologyGet(decomposition,decompositionTopology,err,error,*999) + CALL DecompositionTopology_ElementCheckExists(decompositionTopology,projectionElementUserNumber,elementExists, & + & elementLocalNumber,ghostElement,err,error,*999) + IF(elementExists) THEN + IF(.NOT.ghostElement) THEN + dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementLocalNumber=elementLocalNumber + dataProjection%dataProjectionResults(dataPointGlobalNumber)%exitTag=DATA_PROJECTION_USER_SPECIFIED + ENDIF + ELSE + localError="Element user number "//TRIM(NumberToVString(projectionElementUserNumber,"*",err,error))// & + & " does not exist." + CALL FlagError(localError,err,error,*999) ENDIF - ELSE - localError="Element user number "//TRIM(NumberToVString(projectionElementUserNumber,"*",err,error))// & - & " does not exist." + CASE DEFAULT + localError="The data projection type of "//TRIM(NumberToVString(dataProjection%projectionType,"",err,error))// & + & " is invalid." CALL FlagError(localError,err,error,*999) - ENDIF - + END SELECT + EXITS("DataProjection_ResultElementNumberSet") RETURN 999 ERRORSEXITS("DataProjection_ResultElementNumberSet",err,error) @@ -5229,7 +5242,7 @@ SUBROUTINE DataProjection_ResultElementLineNumberSet(dataProjection,dataPointUse SELECT CASE(dataProjection%projectionType) CASE(DATA_PROJECTION_BOUNDARY_FACES_PROJECTION_TYPE) - CALL FlagError("Use the result element faces routine not the element lines routine for a boundary faces projection.", & + CALL FlagError("Use the result element faces routine not the element lines routine for a boundary lines projection.", & & err,error,*999) CASE(DATA_PROJECTION_ALL_ELEMENTS_PROJECTION_TYPE) CALL FlagError("Use the result elements routine not the element lines routine for an all elements projection.", & From ed936d4b29be1499a89885d961fe75288701fea2 Mon Sep 17 00:00:00 2001 From: Chris Bradley Date: Wed, 22 May 2019 11:43:13 +1200 Subject: [PATCH 03/17] Separating out Sobolev smoothing parameters for each direction. --- src/fitting_routines.F90 | 739 +++++++++++++++++++++++---------------- 1 file changed, 446 insertions(+), 293 deletions(-) diff --git a/src/fitting_routines.F90 b/src/fitting_routines.F90 index b9650c12..2d882415 100644 --- a/src/fitting_routines.F90 +++ b/src/fitting_routines.F90 @@ -147,8 +147,8 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* REAL(DP), POINTER :: independentVectorParameters(:),independentWeightParameters(:) REAL(DP) :: projectionXi(3) - REAL(DP) :: porosity0, porosity, permOverVisParam0, permOverVisParam,tauParam,kappaParam - REAL(DP) :: tension,curvature + REAL(DP) :: porosity0, porosity, permOverVisParam0, permOverVisParam + REAL(DP) :: tau1,tau2,tau3,kappa11,kappa22,kappa33,kappa12,kappa13,kappa23,tension,curvature REAL(DP) :: materialFact REAL(DP) :: dXdY(3,3), dXdXi(3,3), dYdXi(3,3), dXidY(3,3), dXidX(3,3) REAL(DP) :: Jxy, Jyxi @@ -339,8 +339,19 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* !Get Sobolev smoothing parameters from interpolated material field CALL Field_InterpolateGauss(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,gaussPointIdx, & & equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) - tauParam=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(1,NO_PART_DERIV) - kappaParam=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(2,NO_PART_DERIV) + tau1=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(1,NO_PART_DERIV) + kappa11=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(2,NO_PART_DERIV) + IF(numberOfXi>1) THEN + tau2=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(3,NO_PART_DERIV) + kappa22=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(4,NO_PART_DERIV) + kappa12=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(5,NO_PART_DERIV) + IF(numberOfXi>2) THEN + tau3=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(6,NO_PART_DERIV) + kappa33=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(7,NO_PART_DERIV) + kappa13=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(8,NO_PART_DERIV) + kappa23=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(9,NO_PART_DERIV) + ENDIF + ENDIF jacobianGaussWeight=equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr%jacobian* & & quadratureScheme%GAUSS_WEIGHTS(gaussPointIdx) @@ -366,50 +377,35 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* dependentParameterColumnIdx=dependentParameterColumnIdx+1 !Calculate Sobolev surface tension and curvature smoothing terms - tension = tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S1, & - & gaussPointIdx)) - curvature = kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S1, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S1_S1, & - & gaussPointIdx)) + tension = tau1*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1, & + & gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx, & + & PART_DERIV_S1,gaussPointIdx) + curvature = kappa11*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx, & + & PART_DERIV_S1_S1,gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS( & + & dependentElementParameterColumnIdx,PART_DERIV_S1_S1,gaussPointIdx) IF(numberOfXi > 1) THEN - tension = tension + tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S2, & - & gaussPointIdx)) - curvature = curvature + kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2_S2, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S2_S2, & - & gaussPointIdx) + & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S2, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S1_S2, & - & gaussPointIdx)) + tension = tension + tau2*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx, & + & PART_DERIV_S2,gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS( & + & dependentElementParameterColumnIdx,PART_DERIV_S2,gaussPointIdx) + curvature = curvature + kappa22*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx, & + & PART_DERIV_S2_S2,gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS( & + & dependentElementParameterColumnIdx,PART_DERIV_S2_S2,gaussPointIdx) + & + & kappa12*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S2, & + & gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx, & + & PART_DERIV_S1_S2,gaussPointIdx) IF(numberOfXi > 2) THEN - tension = tension + tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S3, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S3, & - & gaussPointIdx)) - curvature = curvature + kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S3_S3, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S3_S3, & - & gaussPointIdx)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S3, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S1_S3, & - & gaussPointIdx)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2_S3, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S2_S3, & - & gaussPointIdx)) + tension = tension + tau3*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx, & + & PART_DERIV_S3,gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS( & + & dependentElementParameterColumnIdx,PART_DERIV_S3,gaussPointIdx) + curvature = curvature + kappa33*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,& + & PART_DERIV_S3_S3,gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS( & + & dependentElementParameterColumnIdx,PART_DERIV_S3_S3,gaussPointIdx)+ & + & kappa13*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S3, & + & gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx, & + & PART_DERIV_S1_S3,gaussPointIdx)+ & + & kappa23*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2_S3, & + & gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx, & + & PART_DERIV_S2_S3,gaussPointIdx) ENDIF ! 3D ENDIF ! 2 or 3D sum = (tension + curvature) * jacobianGaussWeight @@ -558,89 +554,92 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* CASE(EQUATIONS_SET_FITTING_NO_SMOOTHING) !Do nothing CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING) - !=========================================== - ! S o b o l e v S m o o t h i n g - !=========================================== - !Loop over gauss points - DO ng=1,quadratureScheme%NUMBER_OF_GAUSS - CALL Field_InterpolateGauss(SECOND_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% & - & geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) - CALL Field_InterpolateGauss(SECOND_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% & - & dependentInterpPoint(dependentVariableType)%ptr,err,error,*999) - CALL Field_InterpolateGauss(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% & - & materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) - CALL Field_InterpolatedPointMetricsCalculate(geometricBasis%NUMBER_OF_XI,equations%interpolation% & - & geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) - tauParam=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(1,NO_PART_DERIV) - kappaParam=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(2,NO_PART_DERIV) - !Loop over field components - jacobianGaussWeight=equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr%jacobian* & - & quadratureScheme%GAUSS_WEIGHTS(ng) - - mhs=0 - DO mh=1,dependentVariable%NUMBER_OF_COMPONENTS - !Loop over element rows - meshComponent1=dependentVariable%components(mh)%MESH_COMPONENT_NUMBER - dependentBasisRow=>dependentField%decomposition%domain(meshComponent1)%ptr% & - & topology%elements%elements(elementNumber)%basis - quadratureSchemeRow=>dependentBasisRow%quadrature%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr - DO ms=1,dependentBasisRow%NUMBER_OF_ELEMENT_PARAMETERS - mhs=mhs+1 - nhs=0 - IF(equationsMatrix%updateMatrix) THEN - !Loop over element columns - DO nh=1,dependentVariable%NUMBER_OF_COMPONENTS - meshComponent2=dependentVariable%components(nh)%MESH_COMPONENT_NUMBER - dependentBasisColumn=>dependentField%decomposition%domain(meshComponent2)%ptr% & - & topology%elements%elements(elementNumber)%basis - quadratureSchemeColumn=>dependentBasisColumn%quadrature%QUADRATURE_SCHEME_MAP( & - & BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr - DO ns=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS - nhs=nhs+1 - sum = 0.0_DP - - !Calculate sobolev surface tension and curvature smoothing terms - tension = tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1,ng)) - curvature = kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S1,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S1,ng)) - - IF(dependentVariable%NUMBER_OF_COMPONENTS > 1) THEN - tension = tension + tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2,ng)) - curvature = curvature + kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S2,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S2,ng) + & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S2,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S2,ng)) - - IF(dependentVariable%NUMBER_OF_COMPONENTS > 2) THEN - tension = tension + tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3,ng)) - curvature = curvature + kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3_S3,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S3,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S3,ng)) - ENDIF ! 3D - ENDIF ! 2 or 3D - - ! Add in smoothing terms to the element matrix - equationsMatrix%elementMatrix%matrix(mhs,nhs) = & - & equationsMatrix%elementMatrix%matrix(mhs,nhs) + (tension + curvature) * jacobianGaussWeight + !=========================================== + ! S o b o l e v S m o o t h i n g + !=========================================== + !Loop over gauss points + DO ng=1,quadratureScheme%NUMBER_OF_GAUSS + CALL Field_InterpolateGauss(SECOND_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% & + & geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + CALL Field_InterpolateGauss(SECOND_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% & + & dependentInterpPoint(dependentVariableType)%ptr,err,error,*999) + CALL Field_InterpolateGauss(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,ng,equations%interpolation% & + & materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + CALL Field_InterpolatedPointMetricsCalculate(geometricBasis%NUMBER_OF_XI,equations%interpolation% & + & geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + tau1=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(1,NO_PART_DERIV) + kappa11=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(2,NO_PART_DERIV) + IF(numberOfXi>1) THEN + tau2=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(3,NO_PART_DERIV) + kappa22=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(4,NO_PART_DERIV) + kappa12=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(5,NO_PART_DERIV) + IF(numberOfXi>2) THEN + tau3=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(6,NO_PART_DERIV) + kappa33=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(7,NO_PART_DERIV) + kappa13=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(8,NO_PART_DERIV) + kappa23=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(9,NO_PART_DERIV) + ENDIF + ENDIF + !Loop over field components + jacobianGaussWeight=equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr%jacobian* & + & quadratureScheme%GAUSS_WEIGHTS(ng) + + mhs=0 + DO mh=1,dependentVariable%NUMBER_OF_COMPONENTS + !Loop over element rows + meshComponent1=dependentVariable%components(mh)%MESH_COMPONENT_NUMBER + dependentBasisRow=>dependentField%decomposition%domain(meshComponent1)%ptr% & + & topology%elements%elements(elementNumber)%basis + quadratureSchemeRow=>dependentBasisRow%quadrature%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr + DO ms=1,dependentBasisRow%NUMBER_OF_ELEMENT_PARAMETERS + mhs=mhs+1 + nhs=0 + IF(equationsMatrix%updateMatrix) THEN + !Loop over element columns + DO nh=1,dependentVariable%NUMBER_OF_COMPONENTS + meshComponent2=dependentVariable%components(nh)%MESH_COMPONENT_NUMBER + dependentBasisColumn=>dependentField%decomposition%domain(meshComponent2)%ptr% & + & topology%elements%elements(elementNumber)%basis + quadratureSchemeColumn=>dependentBasisColumn%quadrature%QUADRATURE_SCHEME_MAP( & + & BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr + DO ns=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS + nhs=nhs+1 + sum = 0.0_DP + + !Calculate Sobolev surface tension and curvature smoothing terms + tension = tau1*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1,ng) + curvature = kappa11*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S1,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S1,ng) + IF(numberOfXi > 1) THEN + tension = tension + tau2*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2,ng) + curvature = curvature + kappa22*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S2,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S2,ng) + & + & kappa12*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S2,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S2,ng) + IF(numberOfXi > 2) THEN + tension = tension + tau3*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3,ng) + curvature = curvature + kappa33*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3_S3,ng)+ & + & kappa13*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S3,ng)+ & + & kappa23*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S3,ng) + ENDIF ! 3D + ENDIF ! 2 or 3D - ENDDO !ns - ENDDO !nh - ENDIF ! update matrix - ENDDO !ms - ENDDO !mh - ENDDO !ng + ! Add in smoothing terms to the element matrix + equationsMatrix%elementMatrix%matrix(mhs,nhs) = & + & equationsMatrix%elementMatrix%matrix(mhs,nhs) + (tension + curvature) * jacobianGaussWeight + + ENDDO !ns + ENDDO !nh + ENDIF ! update matrix + ENDDO !ms + ENDDO !mh + ENDDO !ng CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_FITTING_STRAIN_ENERGY_SMOOTHING) @@ -843,6 +842,7 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* dependentVariableType=fieldVariable%VARIABLE_TYPE dependentBasis=>dependentField%decomposition%domain(dependentField%decomposition%MESH_COMPONENT_NUMBER)%ptr% & & topology%elements%elements(elementNumber)%basis + numberOfXi=dependentBasis%NUMBER_OF_XI geometricBasis=>geometricField%decomposition%domain(geometricField%decomposition%MESH_COMPONENT_NUMBER)%ptr% & & topology%elements%elements(elementNumber)%basis sourceBasis=>sourceField%decomposition%domain(sourceField%decomposition%MESH_COMPONENT_NUMBER)%ptr% & @@ -866,9 +866,20 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* & materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) CALL Field_InterpolatedPointMetricsCalculate(geometricBasis%NUMBER_OF_XI,equations%interpolation% & & geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) - tauParam=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(1,NO_PART_DERIV) - kappaParam=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(2,NO_PART_DERIV) - ! WRITE(*,*)'tauParam ',tauParam + tau1=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(1,NO_PART_DERIV) + kappa11=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(2,NO_PART_DERIV) + IF(numberOfXi>1) THEN + tau2=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(3,NO_PART_DERIV) + kappa22=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(4,NO_PART_DERIV) + kappa12=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(5,NO_PART_DERIV) + IF(numberOfXi>2) THEN + tau3=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(6,NO_PART_DERIV) + kappa33=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(7,NO_PART_DERIV) + kappa13=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(8,NO_PART_DERIV) + kappa23=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(9,NO_PART_DERIV) + ENDIF + ENDIF + ! WRITE(*,*)'tauParam ',tauParam uValue=0.0_DP IF(sourceVector%updateVector) THEN CALL Field_InterpolationParametersElementGet(FIELD_VALUES_SET_TYPE,elementNumber, & @@ -934,29 +945,29 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* CASE(EQUATIONS_SET_FITTING_NO_SMOOTHING) !Do nothing CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING) - sum = sum + ( & - & tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3,ng)) +& - & kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S1,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S1,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S2,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S2,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3_S3,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S2,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S2,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S3,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S3,ng))) !& - ! no weighting either? - ! & * rwg + tension = tau1*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1,ng) + curvature = kappa11*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S1,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S1,ng) + IF(numberOfXi > 1) THEN + tension = tension + tau2*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2,ng) + curvature = curvature + kappa22*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S2,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S2,ng) + & + & kappa12*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S2,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S2,ng) + IF(numberOfXi > 2) THEN + tension = tension + tau3*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3,ng) + curvature = curvature + kappa33*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3_S3,ng)+ & + & kappa13*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S3,ng)+ & + & kappa23*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S3,ng) + ENDIF ! 3D + ENDIF ! 2 or 3D + sum = sum + (tension + curvature)*rwg CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) @@ -979,32 +990,29 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* CASE(EQUATIONS_SET_FITTING_NO_SMOOTHING) !Do nothing CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING) - !REDUCED SOBOLEV SMOOTHING - !This stiffness matrix contribution is with "integration" means ng=ng in fact! - sum = sum + ( & - & tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3,ng)) +& - & kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S1,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S1,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S2,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S2,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3_S3,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S2,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S2,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S3,ng)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S3,ng)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S3,ng))) !& - ! no weighting either? - ! & * rwg - + tension = tau1*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1,ng) + curvature = kappa11*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S1,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S1,ng) + IF(numberOfXi > 1) THEN + tension = tension + tau2*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2,ng) + curvature = curvature + kappa22*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S2,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S2,ng) + & + & kappa12*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S2,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S2,ng) + IF(numberOfXi > 2) THEN + tension = tension + tau3*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3,ng) + curvature = curvature + kappa33*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S3_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S3_S3,ng)+ & + & kappa13*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S1_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S1_S3,ng)+ & + & kappa23*quadratureSchemeRow%GAUSS_BASIS_FNS(ms,PART_DERIV_S2_S3,ng)* & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(ns,PART_DERIV_S2_S3,ng) + ENDIF ! 3D + ENDIF ! 2 or 3D + sum = sum + (tension + curvature)*rwg CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) @@ -1148,8 +1156,19 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* !Get Sobolev smoothing data from interpolated fields CALL Field_InterpolateGauss(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,gaussPointIdx, & & equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) - tauParam=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(1,NO_PART_DERIV) - kappaParam=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(2,NO_PART_DERIV) + tau1=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(1,NO_PART_DERIV) + kappa11=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(2,NO_PART_DERIV) + IF(numberOfXi>1) THEN + tau2=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(3,NO_PART_DERIV) + kappa22=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(4,NO_PART_DERIV) + kappa12=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(5,NO_PART_DERIV) + IF(numberOfXi>2) THEN + tau3=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(6,NO_PART_DERIV) + kappa33=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(7,NO_PART_DERIV) + kappa13=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(8,NO_PART_DERIV) + kappa23=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(9,NO_PART_DERIV) + ENDIF + ENDIF CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_FITTING_STRAIN_ENERGY_SMOOTHING) @@ -1194,53 +1213,40 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* !Do nothing CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING) !Calculate Sobolev surface tension and curvature smoothing terms - tension = tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S1, & - & gaussPointIdx)) - curvature = kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S1, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S1_S1, & - & gaussPointIdx)) + tension = tau1*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1, & + & gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx, & + & PART_DERIV_S1,gaussPointIdx) + curvature = kappa11*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx, & + & PART_DERIV_S1_S1,gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS( & + & dependentElementParameterColumnIdx,PART_DERIV_S1_S1,gaussPointIdx) IF(numberOfXi > 1) THEN - tension = tension + tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S2, & - & gaussPointIdx)) - curvature = curvature + kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2_S2, & - & gaussPointIdx)* & + tension = tension + tau2*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx, & + & PART_DERIV_S2,gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS( & + & dependentElementParameterColumnIdx,PART_DERIV_S2,gaussPointIdx) + curvature = curvature + kappa22*quadratureSchemeRow%GAUSS_BASIS_FNS( & + & dependentElementParameterRowIdx,PART_DERIV_S2_S2,gaussPointIdx)* & & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S2_S2, & & gaussPointIdx) + & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S2, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S1_S2, & - & gaussPointIdx)) + & kappa12*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S2, & + & gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx, & + & PART_DERIV_S1_S2,gaussPointIdx) IF(numberOfXi > 2) THEN - tension = tension + tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S3, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S3, & - & gaussPointIdx)) - curvature = curvature + kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S3_S3, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S3_S3, & - & gaussPointIdx)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S3, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S1_S3, & - & gaussPointIdx)+ & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2_S3, & - & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S2_S3, & - & gaussPointIdx)) + tension = tension + tau3*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx, & + & PART_DERIV_S3,gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS( & + & dependentElementParameterColumnIdx,PART_DERIV_S3,gaussPointIdx) + curvature = curvature + kappa33*quadratureSchemeRow%GAUSS_BASIS_FNS( & + & dependentElementParameterRowIdx, & + & PART_DERIV_S3_S3,gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS( & + & dependentElementParameterColumnIdx,PART_DERIV_S3_S3,gaussPointIdx)+ & + & kappa13*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S3, & + & gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx, & + & PART_DERIV_S1_S3,gaussPointIdx)+ & + & kappa23*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2_S3, & + & gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx, & + & PART_DERIV_S2_S3,gaussPointIdx) ENDIF ! 3D ENDIF ! 2 or 3D - sum = sum + (tension + curvature) * jacobianGaussWeight + sum = (tension + curvature) * jacobianGaussWeight CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_FITTING_STRAIN_ENERGY_SMOOTHING) @@ -2192,7 +2198,7 @@ SUBROUTINE Fitting_EquationsSetGaussSetup(equationsSet,equationsSetSetup,err,err TYPE(VARYING_STRING), INTENT(OUT) :: error ! Date: Sat, 1 Jun 2019 16:01:20 +1200 Subject: [PATCH 04/17] Fixing fit plus adding specifying fields for export. --- src/field_routines.F90 | 75 ++++++++++++++++++++++++++++++++++++++++ src/fitting_routines.F90 | 12 +++---- src/opencmiss_iron.F90 | 28 ++++++++++++++- 3 files changed, 108 insertions(+), 7 deletions(-) diff --git a/src/field_routines.F90 b/src/field_routines.F90 index d275080a..814629b0 100644 --- a/src/field_routines.F90 +++ b/src/field_routines.F90 @@ -1290,6 +1290,8 @@ MODULE FIELD_ROUTINES PUBLIC FieldVariable_ParameterSetGet + PUBLIC Fields_AddField + PUBLIC Fields_Finalise,Fields_Initialise PUBLIC MESH_EMBEDDING_PUSH_DATA, MESH_EMBEDDING_PULL_GAUSS_POINT_DATA, FIELD_PARAMETER_SET_GET_GAUSS_POINT_COORD @@ -32067,6 +32069,79 @@ END SUBROUTINE FieldVariable_ParameterSetsCopyIfExists !================================================================================================================================ ! + !>Adds a field to a fields list + SUBROUTINE Fields_AddField(fields,field,err,error,*) + + !Argument variables + TYPE(FIELDS_TYPE), POINTER :: fields !0) THEN + IF(ASSOCIATED(fields%region)) THEN + IF(ASSOCIATED(field%region)) THEN + IF(.NOT.ASSOCIATED(fields%region,field%region)) & + & CALL FlagError("The specified field does not have the same region as the specified fields.",err,error,*999) + ELSE IF(ASSOCIATED(field%INTERFACE)) THEN + CALL FlagError("Can not add a field from an interface to fields that are from a region.",err,error,*999) + ELSE + CALL FlagError("Field does not have an associated region or interface.",err,error,*999) + ENDIF + ELSE IF(ASSOCIATED(fields%INTERFACE)) THEN + IF(ASSOCIATED(field%interface)) THEN + IF(.NOT.ASSOCIATED(fields%interface,field%interface)) & + & CALL FlagError("The specified field does not have the same interface as the specified fields.",err,error,*999) + ELSE IF(ASSOCIATED(field%region)) THEN + CALL FlagError("Can not add a field from a region to fields that are from an interface.",err,error,*999) + ELSE + CALL FlagError("Field does not have an associated region or interface.",err,error,*999) + ENDIF + ELSE + CALL FlagError("Fields does not have a region or interface.",err,error,*999) + ENDIF + ELSE + IF(ASSOCIATED(field%region)) THEN + fields%region=>field%region + ELSE IF(ASSOCIATED(field%INTERFACE)) THEN + fields%interface=>field%interface + ELSE + CALL FlagError("Field does not have an associated region or interface.",err,error,*999) + ENDIF + ENDIF + DO fieldIdx=1,fields%NUMBER_OF_FIELDS + newFields(fieldIdx)%ptr=>fields%fields(fieldIdx)%ptr + ENDDO !fieldIdx + newFields(fields%NUMBER_OF_FIELDS+1)%ptr=>field + IF(ASSOCIATED(fields%fields)) DEALLOCATE(fields%fields) + fields%fields=>newFields + fields%NUMBER_OF_FIELDS=fields%NUMBER_OF_FIELDS+1 + + EXITS("Fields_AddField") + RETURN +999 ERRORSEXITS("Fields_AddField",err,error) + RETURN 1 + + END SUBROUTINE Fields_AddField + + ! + !================================================================================================================================ + ! + !>Finalises the fields and deallocates all memory. SUBROUTINE FIELDS_FINALISE(FIELDS,ERR,ERROR,*) diff --git a/src/fitting_routines.F90 b/src/fitting_routines.F90 index 2d882415..b12f4388 100644 --- a/src/fitting_routines.F90 +++ b/src/fitting_routines.F90 @@ -257,12 +257,12 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* dataPointGlobalNumber = dataPoints%elementDataPoint(elementNumber)%dataIndices(dataPointIdx)%globalNumber ! Need to use global number to get the correct projection results projectionXi(1:numberOfXi) = dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementXi(1:numberOfXi) - CALL Field_InterpolateXi(FIRST_PART_DERIV,projectionXi,equations%interpolation% & - & geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) - CALL Field_InterpolateXi(FIRST_PART_DERIV,projectionXi,equations%interpolation% & - & dependentInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) - CALL Field_InterpolatedPointMetricsCalculate(geometricBasis%NUMBER_OF_XI,equations%interpolation% & - & geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + !CALL Field_InterpolateXi(FIRST_PART_DERIV,projectionXi,equations%interpolation% & + ! & geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + !CALL Field_InterpolateXi(FIRST_PART_DERIV,projectionXi,equations%interpolation% & + ! & dependentInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + !CALL Field_InterpolatedPointMetricsCalculate(geometricBasis%NUMBER_OF_XI,equations%interpolation% & + ! & geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) !Get data point vector value and weight DO componentIdx=1,numberOfDataComponents localDof=dataVariable%components(componentIdx)%PARAM_TO_DOF_MAP% & diff --git a/src/opencmiss_iron.F90 b/src/opencmiss_iron.F90 index 62ec0edf..9419b604 100644 --- a/src/opencmiss_iron.F90 +++ b/src/opencmiss_iron.F90 @@ -380,7 +380,7 @@ MODULE OpenCMISS_Iron PUBLIC cmfe_FieldType,cmfe_Field_Finalise,cmfe_Field_Initialise - PUBLIC cmfe_FieldsType,cmfe_Fields_Create,cmfe_Fields_Finalise,cmfe_Fields_Initialise + PUBLIC cmfe_FieldsType,cmfe_Fields_AddField,cmfe_Fields_Create,cmfe_Fields_Finalise,cmfe_Fields_Initialise PUBLIC cmfe_GeneratedMeshType,cmfe_GeneratedMesh_Finalise,cmfe_GeneratedMesh_Initialise @@ -8577,6 +8577,32 @@ END SUBROUTINE cmfe_Field_Initialise !================================================================================================================================ ! + !>Creates a cmfe_FieldsType object for an inteface by an object reference. + SUBROUTINE cmfe_Fields_AddField(fields,field,err) + !DLLEXPORT(cmfe_Fields_CreateInterface) + + !Argument variables + TYPE(cmfe_FieldsType), INTENT(INOUT) :: fields !Creates a cmfe_FieldsType object for an inteface by an object reference. SUBROUTINE cmfe_Fields_CreateInterface(interface,fields,err) !DLLEXPORT(cmfe_Fields_CreateInterface) From d08ba05635293a3a434620e008182c6112b14b4c Mon Sep 17 00:00:00 2001 From: Chris Bradley Date: Sun, 2 Jun 2019 18:55:35 +1200 Subject: [PATCH 05/17] Correct symmetric Hessian entries to avoid NaNs --- src/data_projection_routines.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/data_projection_routines.F90 b/src/data_projection_routines.F90 index 8c08d12b..1c35fbeb 100644 --- a/src/data_projection_routines.F90 +++ b/src/data_projection_routines.F90 @@ -2541,6 +2541,7 @@ SUBROUTINE DataProjection_NewtonElementsEvaluate_2(dataProjection,interpolatedPo & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1_S2))- & & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1), & & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + functionHessian(2,1)=functionHessian(1,2) functionHessian(2,2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2_S2))- & & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2), & @@ -2790,6 +2791,7 @@ SUBROUTINE DataProjection_NewtonElementsEvaluate_3(dataProjection,interpolatedPo & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1_S3))- & & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1), & & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S3))) + functionHessian(2,1)=functionHessian(1,2) functionHessian(2,2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2_S2))- & & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2), & @@ -2798,6 +2800,8 @@ SUBROUTINE DataProjection_NewtonElementsEvaluate_3(dataProjection,interpolatedPo & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2_S3))- & & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2), & & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S3))) + functionHessian(3,1)=functionHessian(1,3) + functionHessian(3,2)=functionHessian(2,3) functionHessian(3,3)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S3_S3))- & & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S3), & From b372f22f6222bb6e24a092b5954bf64b26cf195b Mon Sep 17 00:00:00 2001 From: Chris Bradley Date: Wed, 5 Jun 2019 10:38:18 +1200 Subject: [PATCH 06/17] Adding Sobolev difference smoothing --- src/fitting_routines.F90 | 94 ++++++++++++++++++++++++++++++---------- 1 file changed, 72 insertions(+), 22 deletions(-) diff --git a/src/fitting_routines.F90 b/src/fitting_routines.F90 index b12f4388..1ac6819d 100644 --- a/src/fitting_routines.F90 +++ b/src/fitting_routines.F90 @@ -26,7 +26,7 @@ !> Auckland, the University of Oxford and King's College, London. !> All Rights Reserved. !> -!> Contributor(s): Chris Bradley +!> Contributor(s): Chris Bradley,Sebastian Krittian !> !> Alternatively, the contents of this file may be used under the terms of !> either the GNU General Public License Version 2 or later (the "GPL"), or @@ -119,7 +119,7 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* INTEGER(INTG) :: ng,mh,mhs,ms,nh,nhs,ns,mi,ni INTEGER(INTG) :: dependentComponentColumnIdx,dependentComponentRowIdx,dependentElementParameterColumnIdx, & & dependentElementParameterRowIdx,dependentParameterColumnIdx,dependentParameterRowIdx,gaussPointIdx, & - & meshComponentRow,meshComponentColumn,numberOfDataComponents + & geometricDerivative,meshComponentRow,meshComponentColumn,numberOfDataComponents REAL(DP) :: rwg,sum,jacobianGaussWeight REAL(DP) :: basisFunctionRow,basisFunctionColumn,PGM,PGN,PGMSI(3),PGNSI(3) REAL(DP) :: uValue(3) @@ -145,22 +145,21 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* TYPE(QUADRATURE_SCHEME_TYPE), POINTER :: quadratureScheme,quadratureSchemeColumn,quadratureSchemeRow TYPE(VARYING_STRING) :: localError - REAL(DP), POINTER :: independentVectorParameters(:),independentWeightParameters(:) - REAL(DP) :: projectionXi(3) - REAL(DP) :: porosity0, porosity, permOverVisParam0, permOverVisParam - REAL(DP) :: tau1,tau2,tau3,kappa11,kappa22,kappa33,kappa12,kappa13,kappa23,tension,curvature - REAL(DP) :: materialFact - REAL(DP) :: dXdY(3,3), dXdXi(3,3), dYdXi(3,3), dXidY(3,3), dXidX(3,3) - REAL(DP) :: Jxy, Jyxi - REAL(DP) :: dataPointWeight(99),dataPointVector(99) INTEGER(INTG) :: derivative_idx, component_idx, xi_idx, numberOfDimensions,smoothingType INTEGER(INTG) :: dataPointIdx,dataPointUserNumber,dataPointLocalNumber,dataPointGlobalNumber INTEGER(INTG) :: numberOfXi INTEGER(INTG) :: componentIdx INTEGER(INTG) :: dependentVariableType,variableType,localDof - INTEGER(INTG) numberDofs INTEGER(INTG) meshComponent1,meshComponent2 + REAL(DP) :: projectionXi(3) + REAL(DP) :: porosity0, porosity, permOverVisParam0, permOverVisParam + REAL(DP) :: tau1,tau2,tau3,kappa11,kappa22,kappa33,kappa12,kappa13,kappa23,tension,curvature,rhsTension,rhsCurvature + REAL(DP) :: materialFact + REAL(DP) :: dXdY(3,3), dXdXi(3,3), dYdXi(3,3), dXidY(3,3), dXidX(3,3) + REAL(DP) :: Jxy, Jyxi + REAL(DP) :: dataPointWeight(99),dataPointVector(99) + REAL(DP), POINTER :: independentVectorParameters(:),independentWeightParameters(:) ENTERS("Fitting_FiniteElementCalculate",err,error,*999) @@ -322,15 +321,21 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* SELECT CASE(smoothingType) CASE(EQUATIONS_SET_FITTING_NO_SMOOTHING) !Do nothing - CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING) + CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING, & + & EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) IF(equationsMatrix%updateMatrix) THEN materialsField=>equations%interpolation%materialsField CALL Field_InterpolationParametersElementGet(FIELD_VALUES_SET_TYPE,elementNumber,equations%interpolation% & & materialsInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + IF(smoothingType==EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) THEN + geometricDerivative=SECOND_PART_DERIV + ELSE + geometricDerivative=FIRST_PART_DERIV + ENDIF !Loop over Gauss points DO gaussPointIdx=1,quadratureScheme%NUMBER_OF_GAUSS !Interpolate fields - CALL Field_InterpolateGauss(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,gaussPointIdx, & + CALL Field_InterpolateGauss(geometricDerivative,BASIS_DEFAULT_QUADRATURE_SCHEME,gaussPointIdx, & & equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) CALL Field_InterpolateGauss(SECOND_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,gaussPointIdx, & & equations%interpolation%dependentInterpPoint(dependentVariableType)%ptr,err,error,*999) @@ -408,20 +413,67 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* & PART_DERIV_S2_S3,gaussPointIdx) ENDIF ! 3D ENDIF ! 2 or 3D - sum = (tension + curvature) * jacobianGaussWeight + sum = 2.0_DP*(tension + curvature) * jacobianGaussWeight equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)= & equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)+sum ENDDO !dependentElementParameterColumnIdx ENDDO !dependentComponentColumnIdx + IF(smoothingType==EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) THEN + IF(rhsVector%updateVector) THEN + !Calculate Sobolev surface tension and curvature smoothing terms + rhsTension=0.0_DP + rhsCurvature=0.0_DP + DO dependentComponentColumnIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + rhsTension = rhsTension+ & + & tau1*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1, & + & gaussPointIdx)*equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr% & + & values(dependentComponentColumnIdx,PART_DERIV_S1) + rhsCurvature = rhsCurvature + & + & kappa11*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S1, & + & gaussPointIdx)*equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr% & + & values(dependentComponentColumnIdx,PART_DERIV_S1_S1) + IF(numberOfXi > 1) THEN + rhsTension = rhsTension + & + & tau2*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2, & + & gaussPointIdx)*equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr% & + & values(dependentComponentColumnIdx,PART_DERIV_S2) + rhsCurvature = rhsCurvature + & + & kappa22*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2_S2, & + & gaussPointIdx)*equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr% & + & values(dependentComponentColumnIdx,PART_DERIV_S2_S2) + & + & kappa12*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S2, & + & gaussPointIdx)*equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr% & + & values(dependentComponentColumnIdx,PART_DERIV_S1_S2) + IF(numberOfXi > 2) THEN + rhsTension = rhsTension + & + & tau3*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S3, & + & gaussPointIdx)*equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr% & + & values(dependentComponentColumnIdx,PART_DERIV_S3) + rhsCurvature = rhsCurvature + & + & kappa33*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S3_S3, & + & gaussPointIdx)*equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr% & + & values(dependentComponentColumnIdx,PART_DERIV_S3_S3) + & + & kappa13*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S3, & + & gaussPointIdx)*equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr% & + & values(dependentComponentColumnIdx,PART_DERIV_S1_S3) + & + & kappa23*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2_S3, & + & gaussPointIdx)*equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr% & + & values(dependentComponentColumnIdx,PART_DERIV_S2_S3) + ENDIF ! 3D + ENDIF ! 2 or 3D + ENDDO !dependentComponentColumnIdx + sum = 2.0_DP*(rhsTension + rhsCurvature) * jacobianGaussWeight + rhsVector%elementVector%vector(dependentParameterRowIdx)= & + & rhsVector%elementVector%vector(dependentParameterRowIdx)+sum + ENDIF + ENDIF ENDDO !dependentElementParameterRowIdx ENDDO !dependentComponentRowIdx ENDDO !gaussPointIdx ENDIF - CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) - CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_FITTING_STRAIN_ENERGY_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) CASE DEFAULT @@ -1152,7 +1204,8 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* SELECT CASE(smoothingType) CASE(EQUATIONS_SET_FITTING_NO_SMOOTHING) !Do nothing - CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING) + CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING, & + & EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) !Get Sobolev smoothing data from interpolated fields CALL Field_InterpolateGauss(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,gaussPointIdx, & & equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) @@ -1169,8 +1222,6 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* kappa23=equations%interpolation%materialsInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr%values(9,NO_PART_DERIV) ENDIF ENDIF - CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) - CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_FITTING_STRAIN_ENERGY_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) CASE DEFAULT @@ -1211,7 +1262,8 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* SELECT CASE(smoothingType) CASE(EQUATIONS_SET_FITTING_NO_SMOOTHING) !Do nothing - CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING) + CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING, & + & EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) !Calculate Sobolev surface tension and curvature smoothing terms tension = tau1*quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1, & & gaussPointIdx)*quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx, & @@ -1247,8 +1299,6 @@ SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,* ENDIF ! 3D ENDIF ! 2 or 3D sum = (tension + curvature) * jacobianGaussWeight - CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) - CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_FITTING_STRAIN_ENERGY_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) CASE DEFAULT From 249060b99fbb478f68f8b8d00ceb121be89a700c Mon Sep 17 00:00:00 2001 From: Chris Bradley Date: Fri, 7 Jun 2019 13:30:33 +1200 Subject: [PATCH 07/17] Output number of canceled data points --- src/data_projection_routines.F90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/data_projection_routines.F90 b/src/data_projection_routines.F90 index 1c35fbeb..f2a607c2 100644 --- a/src/data_projection_routines.F90 +++ b/src/data_projection_routines.F90 @@ -4780,7 +4780,8 @@ SUBROUTINE DataProjection_ResultAnalysisOutput(dataProjection,filename,err,error TYPE(VARYING_STRING), INTENT(OUT) :: error ! Date: Sat, 8 Jun 2019 18:27:02 +1200 Subject: [PATCH 08/17] Fixes to data projection analysis output. --- src/data_projection_routines.F90 | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/data_projection_routines.F90 b/src/data_projection_routines.F90 index f2a607c2..b839ae5e 100644 --- a/src/data_projection_routines.F90 +++ b/src/data_projection_routines.F90 @@ -26,7 +26,7 @@ !> Auckland, the University of Oxford and King's College, London. !> All Rights Reserved. !> -!> Contributor(s): Chris Bradley, Kumar Mithraratne, Xiani (Nancy) Yan, Prasad Babarenda Gamage +!> Contributor(s): Chris Bradley, Tim Wu, Kumar Mithraratne, Xiani (Nancy) Yan, Prasad Babarenda Gamage !> !> Alternatively, the contents of this file may be used under the terms of !> either the GNU General Public License Version 2 or later (the "GPL"), or @@ -4927,19 +4927,23 @@ SUBROUTINE DataProjection_ResultAnalysisOutput(dataProjection,filename,err,error ENDIF ENDDO !dataPointIdx CALL WriteString(outputID,"",err,error,*999) - CALL WriteStringValue(outputID,"Number of canceled data points = ",numberOfCanceledDataPoints,err,error,*999) + CALL WriteStringValue(outputID," Number of canceled data points = ",numberOfCanceledDataPoints,err,error,*999) CALL WriteString(outputID,"",err,error,*999) CALL WriteString(outputID," Errors:",err,error,*999) CALL WriteString(outputID,"",err,error,*999) CALL WriteString(outputID," Error type Value Data user#",err,error,*999) WRITE(localString,'(" RMS error ",2X,E12.5)') dataProjection%rmsError CALL WriteString(outputID,localString,err,error,*999) - WRITE(localString,'(" Maximum error",2X,E12.5,4X,I8)') dataProjection%maximumError, & - & dataProjection%dataProjectionResults(dataProjection%maximumErrorDataPoint)%userNumber - CALL WriteString(outputID,localString,err,error,*999) - WRITE(localString,'(" Minimum error",2X,E12.5,4X,I8)') dataProjection%minimumError, & - & dataProjection%dataProjectionResults(dataProjection%minimumErrorDataPoint)%userNumber - CALL WriteString(outputID,localString,err,error,*999) + IF(dataProjection%maximumErrorDataPoint/=0) THEN + WRITE(localString,'(" Maximum error",2X,E12.5,4X,I8)') dataProjection%maximumError, & + & dataProjection%dataProjectionResults(dataProjection%maximumErrorDataPoint)%userNumber + CALL WriteString(outputID,localString,err,error,*999) + ENDIF + IF(dataProjection%minimumErrorDataPoint/=0) THEN + WRITE(localString,'(" Minimum error",2X,E12.5,4X,I8)') dataProjection%minimumError, & + & dataProjection%dataProjectionResults(dataProjection%minimumErrorDataPoint)%userNumber + CALL WriteString(outputID,localString,err,error,*999) + ENDIF CALL WriteString(outputID,"",err,error,*999) ENDIF IF(filenameLength>=1) THEN From 5817d3628f1a581d76c5e454c456f0def817a595 Mon Sep 17 00:00:00 2001 From: Chris Bradley Date: Sun, 9 Jun 2019 15:13:07 +1200 Subject: [PATCH 09/17] Allow for interpolating a field at a projected data point. --- src/data_projection_routines.F90 | 75 ++++++++++++++++ src/opencmiss_iron.F90 | 141 +++++++++++++++++++++++++++++++ 2 files changed, 216 insertions(+) diff --git a/src/data_projection_routines.F90 b/src/data_projection_routines.F90 index b839ae5e..262ce134 100644 --- a/src/data_projection_routines.F90 +++ b/src/data_projection_routines.F90 @@ -142,6 +142,8 @@ MODULE DataProjectionRoutines PUBLIC DataProjection_Destroy + PUBLIC DataProjection_DataPointFieldEvaluate + PUBLIC DataProjection_DataPointsProjectionEvaluate PUBLIC DataProjection_DataPointsPositionEvaluate @@ -1325,6 +1327,79 @@ SUBROUTINE DataProjection_Initialise(dataProjection,err,error,*) END SUBROUTINE DataProjection_Initialise + ! + !================================================================================================================================ + ! + + !>Evaluate the data point in a field based on data projection + SUBROUTINE DataProjection_DataPointFieldEvaluate(dataProjection,field,fieldVariableType,fieldParameterSetType, & + & dataPointUserNumber,fieldResult,err,error,*) + + !Argument variables + TYPE(DataProjectionType), POINTER :: dataProjection !interpolatedPoints(fieldVariableType)%ptr + + elementNumber=dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementLocalNumber + CALL Field_InterpolationParametersElementGet(fieldParameterSetType,elementNumber, & + & interpolationParameters(fieldVariableType)%ptr,err,error,*999) + CALL Field_InterpolateXi(NO_PART_DERIV,dataProjection%dataProjectionResults(dataPointGlobalNumber)%xi, & + & interpolatedPoint,err,error,*999) + DO coordinateIdx=1,fieldVariable%NUMBER_OF_COMPONENTS + fieldResult(coordinateIdx)=interpolatedPoint%values(coordinateIdx,NO_PART_DERIV) + ENDDO !coordinateIdx + + CALL Field_InterpolatedPointsFinalise(interpolatedPoints,err,error,*999) + CALL Field_InterpolationParametersFinalise(interpolationParameters,err,error,*999) + + EXITS("DataProjection_DataPointFieldEvaluate") + RETURN +999 ERRORSEXITS("DataProjection_DataPointFieldEvaluate",err,error) + RETURN 1 + + END SUBROUTINE DataProjection_DataPointFieldEvaluate + ! !================================================================================================================================ ! diff --git a/src/opencmiss_iron.F90 b/src/opencmiss_iron.F90 index 9419b604..66aeeefc 100644 --- a/src/opencmiss_iron.F90 +++ b/src/opencmiss_iron.F90 @@ -1893,6 +1893,13 @@ MODULE OpenCMISS_Iron MODULE PROCEDURE cmfe_DataProjection_DestroyObj END INTERFACE cmfe_DataProjection_Destroy + !>Evaluate A data point in a field based on data projection + INTERFACE cmfe_DataProjection_DataPointFieldEvaluate + MODULE PROCEDURE cmfe_DataProjection_DataPointFieldEvaluateRegionNumber + MODULE PROCEDURE cmfe_DataProjection_DataPointFieldEvaluateInterfaceNumber + MODULE PROCEDURE cmfe_DataProjection_DataPointFieldEvaluateObj + END INTERFACE cmfe_DataProjection_DataPointFieldEvaluate + !>Evaluate the data points position in a field based on data projection INTERFACE cmfe_DataProjection_DataPointsPositionEvaluate MODULE PROCEDURE cmfe_DataProjection_DataPointsPositionEvaluateRegionNumber @@ -2206,6 +2213,8 @@ MODULE OpenCMISS_Iron PUBLIC cmfe_DataProjection_Destroy + PUBLIC cmfe_DataProjection_DataPointFieldEvaluate + PUBLIC cmfe_DataProjection_DataPointsPositionEvaluate PUBLIC cmfe_DataProjection_ProjectionCancelByDataPoints @@ -20179,6 +20188,138 @@ END SUBROUTINE cmfe_DataProjection_DestroyObj !================================================================================================================================ ! + !>Evaluate the data point in a field based on data projection in a region, identified by user number + SUBROUTINE cmfe_DataProjection_DataPointFieldEvaluateRegionNumber(regionUserNumber,dataPointsUserNumber, & + & dataProjectionUserNumber,fieldUserNumber,fieldVariableType,fieldParameterSetType,dataPointUserNumber,fieldResult,err) + !DLLEXPORT(cmfe_DataProjection_DataPointFieldEvaluateRegionNumber) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: regionUserNumber !Evaluate the data point in a field based on data projection in an interface, identified by user number + SUBROUTINE cmfe_DataProjection_DataPointFieldEvaluateInterfaceNumber(parentRegionUserNumber,interfaceUserNumber, & + & dataPointsUserNumber,dataProjectionUserNumber,fieldUserNumber,fieldVariableType,fieldParameterSetType, & + & dataPointUserNumber,fieldResult,err) + !DLLEXPORT(cmfe_DataProjection_DataPointFieldEvaluateInterfaceNumber) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: parentRegionUserNumber !Evaluate the data point in a field based on data projection, identified by object + SUBROUTINE cmfe_DataProjection_DataPointFieldEvaluateObj(dataProjection,field,fieldVariableType,fieldParameterSetType, & + & dataPointUserNumber,fieldResult,err) + !DLLEXPORT(cmfe_DataProjection_DataPointFieldEvaluateObj) + + !Argument variables + TYPE(cmfe_DataProjectionType), INTENT(INOUT) :: dataProjection !Evaluate the data points position in a field based on data projection in a region, identified by user number SUBROUTINE cmfe_DataProjection_DataPointsPositionEvaluateRegionNumber(regionUserNumber,dataPointsUserNumber, & & dataProjectionUserNumber,fieldUserNumber,fieldVariableType,fieldParameterSetType,err) From 45f709587eaaaa57052759050fa8b4df64007021 Mon Sep 17 00:00:00 2001 From: Chris Bradley Date: Mon, 10 Jun 2019 09:49:34 +1200 Subject: [PATCH 10/17] Data field evaluate bug fix. --- src/data_projection_routines.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/data_projection_routines.F90 b/src/data_projection_routines.F90 index 262ce134..c6e30468 100644 --- a/src/data_projection_routines.F90 +++ b/src/data_projection_routines.F90 @@ -1384,7 +1384,7 @@ SUBROUTINE DataProjection_DataPointFieldEvaluate(dataProjection,field,fieldVaria elementNumber=dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementLocalNumber CALL Field_InterpolationParametersElementGet(fieldParameterSetType,elementNumber, & & interpolationParameters(fieldVariableType)%ptr,err,error,*999) - CALL Field_InterpolateXi(NO_PART_DERIV,dataProjection%dataProjectionResults(dataPointGlobalNumber)%xi, & + CALL Field_InterpolateXi(NO_PART_DERIV,dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementXi, & & interpolatedPoint,err,error,*999) DO coordinateIdx=1,fieldVariable%NUMBER_OF_COMPONENTS fieldResult(coordinateIdx)=interpolatedPoint%values(coordinateIdx,NO_PART_DERIV) From e6eed064191305385cdf5b9d2ed3360fe6303e86 Mon Sep 17 00:00:00 2001 From: Chris Bradley Date: Tue, 11 Jun 2019 13:10:51 +1200 Subject: [PATCH 11/17] Renaming DataProjection_ResultXiGet and DataProjection_ResultXiSet to DataProjection_ResultProjectionXiGet and DataProjection_ResultProjectionXiSet and adding in DataProjection_ResultElementXiGet. --- src/data_projection_routines.F90 | 63 ++++++++++--- src/opencmiss_iron.F90 | 154 +++++++++++++++++++++++-------- 2 files changed, 169 insertions(+), 48 deletions(-) diff --git a/src/data_projection_routines.F90 b/src/data_projection_routines.F90 index c6e30468..303b5928 100644 --- a/src/data_projection_routines.F90 +++ b/src/data_projection_routines.F90 @@ -190,11 +190,13 @@ MODULE DataProjectionRoutines PUBLIC DataProjection_ResultElementLineNumberGet,DataProjection_ResultElementLineNumberSet + PUBLIC DataProjection_ResultElementXiGet + PUBLIC DataProjection_ResultExitTagGet PUBLIC DataProjection_ResultProjectionVectorGet - PUBLIC DataProjection_ResultXiGet,DataProjection_ResultXiSet + PUBLIC DataProjection_ResultProjectionXiGet,DataProjection_ResultProjectionXiSet PUBLIC DataProjection_StartingXiGet,DataProjection_StartingXiSet @@ -5380,6 +5382,45 @@ END SUBROUTINE DataProjection_ResultElementLineNumberSet !================================================================================================================================ ! + !>Gets the element xi for a data point identified by a given global number. + SUBROUTINE DataProjection_ResultElementXiGet(dataProjection,dataPointUserNumber,elementXi,err,error,*) + + !Argument variables + TYPE(DataProjectionType), POINTER :: dataProjection != "//TRIM(NumberToVString(numberOfElementXi,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + elementXi(1:numberOfElementXi)=dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementXi(1:numberOfElementXi) + + EXITS("DataProjection_ResultElementXiGet") + RETURN +999 ERRORSEXITS("DataProjection_ResultElementXiGet",err,error) + RETURN 1 + + END SUBROUTINE DataProjection_ResultElementXiGet + + ! + !================================================================================================================================ + ! + !>Gets the projection exit tag for a data point identified by a given global number. SUBROUTINE DataProjection_ResultExitTagGet(dataProjection,dataPointUserNumber,projectionExitTag,err,error,*) @@ -5413,7 +5454,7 @@ END SUBROUTINE DataProjection_ResultExitTagGet ! !>Gets the projection xi for a data point identified by a given global number. - SUBROUTINE DataProjection_ResultXiGet(dataProjection,dataPointUserNumber,projectionXi,err,error,*) + SUBROUTINE DataProjection_ResultProjectionXiGet(dataProjection,dataPointUserNumber,projectionXi,err,error,*) !Argument variables TYPE(DataProjectionType), POINTER :: dataProjection !Sets the projection xi for a data point identified by a given global number. - SUBROUTINE DataProjection_ResultXiSet(dataProjection,dataPointUserNumber,projectionXi,err,error,*) + SUBROUTINE DataProjection_ResultProjectionXiSet(dataProjection,dataPointUserNumber,projectionXi,err,error,*) !Argument variables TYPE(DataProjectionType), POINTER :: dataProjection !Returns the element xi for a data point identified by a given user number. + INTERFACE cmfe_DataProjection_ResultElementXiGet + MODULE PROCEDURE cmfe_DataProjection_ResultElementXiGetNumber + MODULE PROCEDURE cmfe_DataProjection_ResultElementXiGetObj + END INTERFACE cmfe_DataProjection_ResultElementXiGet !>Returns the projection exit tag for a data point. INTERFACE cmfe_DataProjection_ResultExitTagGet @@ -2181,16 +2187,16 @@ MODULE OpenCMISS_Iron END INTERFACE cmfe_DataProjection_ResultRMSErrorGet !>Returns the projection xi for a data point identified by a given user number. - INTERFACE cmfe_DataProjection_ResultXiGet - MODULE PROCEDURE cmfe_DataProjection_ResultXiGetNumber - MODULE PROCEDURE cmfe_DataProjection_ResultXiGetObj - END INTERFACE cmfe_DataProjection_ResultXiGet + INTERFACE cmfe_DataProjection_ResultProjectionXiGet + MODULE PROCEDURE cmfe_DataProjection_ResultProjectionXiGetNumber + MODULE PROCEDURE cmfe_DataProjection_ResultProjectionXiGetObj + END INTERFACE cmfe_DataProjection_ResultProjectionXiGet !>Sets the projection xi for a data point identified by a given user number. - INTERFACE cmfe_DataProjection_ResultXiSet - MODULE PROCEDURE cmfe_DataProjection_ResultXiSetNumber - MODULE PROCEDURE cmfe_DataProjection_ResultXiSetObj - END INTERFACE cmfe_DataProjection_ResultXiSet + INTERFACE cmfe_DataProjection_ResultProjectionXiSet + MODULE PROCEDURE cmfe_DataProjection_ResultProjectionXiSetNumber + MODULE PROCEDURE cmfe_DataProjection_ResultProjectionXiSetObj + END INTERFACE cmfe_DataProjection_ResultProjectionXiSet !>Returns the projection vector for a data point identified by a given user number. INTERFACE cmfe_DataProjection_ResultProjectionVectorGet @@ -2263,13 +2269,15 @@ MODULE OpenCMISS_Iron PUBLIC cmfe_DataProjection_ResultElementLineNumberGet,cmfe_DataProjection_ResultElementLineNumberSet + PUBLIC cmfe_DataProjection_ResultElementXiGet + PUBLIC cmfe_DataProjection_ResultExitTagGet PUBLIC cmfe_DataProjection_ResultMaximumErrorGet,cmfe_DataProjection_ResultMinimumErrorGet PUBLIC cmfe_DataProjection_ResultRMSErrorGet - PUBLIC cmfe_DataProjection_ResultXiGet,cmfe_DataProjection_ResultXiSet + PUBLIC cmfe_DataProjection_ResultProjectionXiGet,cmfe_DataProjection_ResultProjectionXiSet PUBLIC cmfe_DataProjection_ResultProjectionVectorGet @@ -23408,6 +23416,74 @@ END SUBROUTINE cmfe_DataProjection_ResultExitTagGetObj !================================================================================================================================ ! + !>Returns the element xi for a data point in a set of data points identified by user number. + SUBROUTINE cmfe_DataProjection_ResultElementXiGetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & + & dataPointUserNumber,elementXi,err) + !DLLEXPORT(cmfe_DataProjection_ResultElementXiGetNumber) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: regionUserNumber !Returns the element xi for a data point in a set of data points identified by an object. + SUBROUTINE cmfe_DataProjection_ResultElementXiGetObj(dataProjection,dataPointUserNumber,elementXi,err) + !DLLEXPORT(cmfe_DataProjection_ResultElementXiGetObj) + + !Argument variables + TYPE(cmfe_DataProjectionType), INTENT(IN) :: dataProjection !Returns the maximum error for a data projection given by numbers. SUBROUTINE cmfe_DataProjection_ResultMaximumErrorGetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & & maximumDataPointUserNumber,maximumError,err) @@ -23611,9 +23687,9 @@ END SUBROUTINE cmfe_DataProjection_ResultRMSErrorGetObj ! !>Returns the projection xi for a data point in a set of data points identified by user number. - SUBROUTINE cmfe_DataProjection_ResultXiGetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & + SUBROUTINE cmfe_DataProjection_ResultProjectionXiGetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & & dataPointUserNumber,projectionXi,err) - !DLLEXPORT(cmfe_DataProjection_ResultXiGetNumber) + !DLLEXPORT(cmfe_DataProjection_ResultProjectionXiGetNumber) !Argument variables INTEGER(INTG), INTENT(IN) :: regionUserNumber !Returns the projection xi for a data point in a set of data points identified by an object. - SUBROUTINE cmfe_DataProjection_ResultXiGetObj(dataProjection,dataPointUserNumber,projectionXi,err) - !DLLEXPORT(cmfe_DataProjection_ResultXiGetObj) + SUBROUTINE cmfe_DataProjection_ResultProjectionXiGetObj(dataProjection,dataPointUserNumber,projectionXi,err) + !DLLEXPORT(cmfe_DataProjection_ResultProjectionXiGetObj) !Argument variables TYPE(cmfe_DataProjectionType), INTENT(IN) :: dataProjection !Sets the projection xi for a data point in a set of data points identified by user number. - SUBROUTINE cmfe_DataProjection_ResultXiSetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & + SUBROUTINE cmfe_DataProjection_ResultProjectionXiSetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & & dataPointUserNumber,projectionXi,err) - !DLLEXPORT(cmfe_DataProjection_ResultXiSetNumber) + !DLLEXPORT(cmfe_DataProjection_ResultProjectionXiSetNumber) !Argument variables INTEGER(INTG), INTENT(IN) :: regionUserNumber !Sets the projection xi for a data point in a set of data points identified by an object. - SUBROUTINE cmfe_DataProjection_ResultXiSetObj(dataProjection,dataPointUserNumber,projectionXi,err) - !DLLEXPORT(cmfe_DataProjection_ResultXiSetObj) + SUBROUTINE cmfe_DataProjection_ResultProjectionXiSetObj(dataProjection,dataPointUserNumber,projectionXi,err) + !DLLEXPORT(cmfe_DataProjection_ResultProjectionXiSetObj) !Argument variables TYPE(cmfe_DataProjectionType), INTENT(IN) :: dataProjection ! Date: Thu, 13 Jun 2019 10:47:47 +1200 Subject: [PATCH 12/17] Bug fix for data point fields in finite elasticity derived variable calculate. --- src/finite_elasticity_routines.F90 | 48 ++++++++++++++++-------------- 1 file changed, 26 insertions(+), 22 deletions(-) diff --git a/src/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90 index c5b8c450..24e104f1 100644 --- a/src/finite_elasticity_routines.F90 +++ b/src/finite_elasticity_routines.F90 @@ -3702,7 +3702,7 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field INTEGER(INTG) :: componentIdx,dependentNumberOfComponents,elementIdx,elementNumber,fieldVariableType,gaussIdx, & & meshComponentNumber,numberOfComponents,numberOfDimensions,numberOfGauss,numberOfTimes,numberOfXi,partIdx, & & startIdx,finishIdx,fieldInterpolation,dataPointNumber,numberOfDataPoints,dataPointIdx,residualVariableType, & - & fieldVarType + & fieldVarType,elementUserNumber REAL(DP) :: dZdNu(3,3),Fg(3,3),Fe(3,3),J,Jg,Je,C(3,3),f(3,3),E(3,3),growthValues(3),xi(3),values(3,3) REAL(SP) :: elementUserElapsed,elementSystemElapsed,systemElapsed,systemTime1(1),systemTime2(1),systemTime3(1), & & systemTime4(1),userElapsed,userTime1(1),userTime2(1),userTime3(1),userTime4(1) @@ -3711,6 +3711,7 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field TYPE(DataProjectionType), POINTER :: dataProjection TYPE(DecompositionDataPointsType), POINTER :: dataPoints TYPE(DECOMPOSITION_TYPE), POINTER :: decomposition + TYPE(DECOMPOSITION_ELEMENTS_TYPE), POINTER :: decompositionElements TYPE(DECOMPOSITION_TOPOLOGY_TYPE), POINTER :: decompositionTopology TYPE(DOMAIN_TYPE), POINTER :: domain TYPE(DOMAIN_ELEMENTS_TYPE), POINTER :: domainElements @@ -3811,6 +3812,8 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field CALL Field_DecompositionGet(dependentField,decomposition,err,error,*999) NULLIFY(decompositionTopology) CALL Decomposition_TopologyGet(decomposition,decompositionTopology,err,error,*999) + NULLIFY(decompositionElements) + CALL DecompositionTopology_ElementsGet(decompositionTopology,decompositionElements,err,error,*999) NULLIFY(domain) CALL Decomposition_DomainGet(decomposition,0,domain,err,error,*999) NULLIFY(domainMappings) @@ -3890,6 +3893,7 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field numberOfTimes=numberOfTimes+1 elementNumber=elementsMappings%DOMAIN_LIST(elementIdx) + elementUserNumber=decompositionElements%elements(elementNumber)%USER_NUMBER IF(diagnostics1) THEN CALL WriteStringValue(DIAGNOSTIC_OUTPUT_TYPE," Element number = ",elementNumber,err,error,*999) @@ -4099,31 +4103,31 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field CASE(3) ! 3 dimensional problem ! ORDER OF THE COMPONENTS: U_11, U_12, U_13, U_22, U_23, U_33 (upper triangular matrix) - CALL Field_ParameterSetUpdateLocalElement(field,fieldVarType,FIELD_VALUES_SET_TYPE, & - & elementNumber,1,values(1,1),err,error,*999) - CALL Field_ParameterSetUpdateLocalElement(field,fieldVarType,FIELD_VALUES_SET_TYPE, & - & elementNumber,2,values(1,2),err,error,*999) - CALL Field_ParameterSetUpdateLocalElement(field,fieldVarType,FIELD_VALUES_SET_TYPE, & - & elementNumber,3,values(1,3),err,error,*999) - CALL Field_ParameterSetUpdateLocalElement(field,fieldVarType,FIELD_VALUES_SET_TYPE, & - & elementNumber,4,values(2,2),err,error,*999) - CALL Field_ParameterSetUpdateLocalElement(field,fieldVarType,FIELD_VALUES_SET_TYPE, & - & elementNumber,5,values(2,3),err,error,*999) - CALL Field_ParameterSetUpdateLocalElement(field,fieldVarType,FIELD_VALUES_SET_TYPE, & - & elementNumber,6,values(3,3),err,error,*999) + CALL Field_ParameterSetUpdateElementDataPoint(field,fieldVarType,FIELD_VALUES_SET_TYPE, & + & elementUserNumber,dataPointIdx,1,values(1,1),err,error,*999) + CALL Field_ParameterSetUpdateElementDataPoint(field,fieldVarType,FIELD_VALUES_SET_TYPE, & + & elementUserNumber,dataPointIdx,2,values(1,2),err,error,*999) + CALL Field_ParameterSetUpdateElementDataPoint(field,fieldVarType,FIELD_VALUES_SET_TYPE, & + & elementUserNumber,dataPointIdx,3,values(1,3),err,error,*999) + CALL Field_ParameterSetUpdateElementDataPoint(field,fieldVarType,FIELD_VALUES_SET_TYPE, & + & elementUserNumber,dataPointIdx,4,values(2,2),err,error,*999) + CALL Field_ParameterSetUpdateElementDataPoint(field,fieldVarType,FIELD_VALUES_SET_TYPE, & + & elementUserNumber,dataPointIdx,5,values(2,3),err,error,*999) + CALL Field_ParameterSetUpdateElementDataPoint(field,fieldVarType,FIELD_VALUES_SET_TYPE, & + & elementUserNumber,dataPointIdx,6,values(3,3),err,error,*999) CASE(2) ! 2 dimensional problem ! ORDER OF THE COMPONENTS: U_11, U_12, U_22 (upper triangular matrix) - CALL Field_ParameterSetUpdateLocalElement(field,fieldVarType,FIELD_VALUES_SET_TYPE, & - & elementNumber,1,values(1,1),err,error,*999) - CALL Field_ParameterSetUpdateLocalElement(field,fieldVarType,FIELD_VALUES_SET_TYPE, & - & elementNumber,2,values(1,2),err,error,*999) - CALL Field_ParameterSetUpdateLocalElement(field,fieldVarType,FIELD_VALUES_SET_TYPE, & - & elementNumber,3,values(2,2),err,error,*999) + CALL Field_ParameterSetUpdateElementDataPoint(field,fieldVarType,FIELD_VALUES_SET_TYPE, & + & elementUserNumber,dataPointIdx,1,values(1,1),err,error,*999) + CALL Field_ParameterSetUpdateElementDataPoint(field,fieldVarType,FIELD_VALUES_SET_TYPE, & + & elementUserNumber,dataPointIdx,2,values(1,2),err,error,*999) + CALL Field_ParameterSetUpdateElementDataPoint(field,fieldVarType,FIELD_VALUES_SET_TYPE, & + & elementUserNumber,dataPointIdx,3,values(2,2),err,error,*999) CASE(1) ! 1 dimensional problem - CALL Field_ParameterSetUpdateLocalElement(field,fieldVarType,FIELD_VALUES_SET_TYPE, & - & elementNumber,1,values(1,1),err,error,*999) + CALL Field_ParameterSetUpdateElementDataPoint(field,fieldVarType,FIELD_VALUES_SET_TYPE, & + & elementUserNumber,dataPointIdx,1,values(1,1),err,error,*999) CASE DEFAULT localError="The number of dimensions of "//TRIM(NumberToVString(numberofDimensions,"*",err,error))// & & " is invalid." @@ -4769,7 +4773,7 @@ SUBROUTINE FiniteElasticity_TensorInterpolateXi(equationsSet,tensorEvaluateType, EXITS("FiniteElasticity_TensorInterpolateXi") RETURN -999 ERRORSEXITS("FiniteElasticity_TensorInterpolateXi",err,error) +999 ERRORSEXITS("Fin7iteElasticity_TensorInterpolateXi",err,error) RETURN 1 END SUBROUTINE FiniteElasticity_TensorInterpolateXi From b2007c16385cd1bacd177f72062f89710bb2e540 Mon Sep 17 00:00:00 2001 From: Chris Bradley Date: Sat, 15 Jun 2019 11:23:10 +1200 Subject: [PATCH 13/17] Return full deformation gradient tensor --- src/finite_elasticity_routines.F90 | 403 ++++++++++++++++++++--------- 1 file changed, 280 insertions(+), 123 deletions(-) diff --git a/src/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90 index 24e104f1..46525954 100644 --- a/src/finite_elasticity_routines.F90 +++ b/src/finite_elasticity_routines.F90 @@ -3700,9 +3700,9 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field TYPE(VARYING_STRING), INTENT(OUT) :: error ! Date: Sun, 11 Aug 2019 00:26:16 +1200 Subject: [PATCH 14/17] Allow different starting xi for datapoints in a dataprojection #176 #177 --- src/data_projection_routines.F90 | 138 +++++++++++++++++++------------ src/opencmiss_iron.F90 | 8 +- src/types.F90 | 2 +- 3 files changed, 89 insertions(+), 59 deletions(-) diff --git a/src/data_projection_routines.F90 b/src/data_projection_routines.F90 index 303b5928..1c0e3d2e 100644 --- a/src/data_projection_routines.F90 +++ b/src/data_projection_routines.F90 @@ -270,7 +270,7 @@ END SUBROUTINE DataProjection_AbsoluteToleranceSet !>Find the closest elements to a data point based on starting xi guess. SUBROUTINE DataProjection_ClosestElementsFind(dataProjection,interpolatedPoint,dataPointLocation,numberOfCandidates, & - & candidateElements,closestElements,closestDistances,err,error,*) + & candidateElements,closestElements,startingXi,closestDistances,err,error,*) !Argument variables TYPE(DataProjectionType), POINTER :: dataProjection !Find the closest faces to a data point base on starting xi guess. SUBROUTINE DataProjection_ClosestFacesFind(dataProjection,interpolatedPoint,dataPointLocation,numberOfCandidates, & - & candidateElements,candidateElementFaces,closestElements,closestElementFaces,closestDistances,err,error,*) + & candidateElements,candidateElementFaces,startingXi,closestElements,closestElementFaces,closestDistances,err,error,*) !Argument variables TYPE(DataProjectionType), POINTER :: dataProjection !Find the closest lines to a data point base on starting xi guess. SUBROUTINE DataProjection_ClosestLinesFind(dataProjection,interpolatedPoint,dataPointLocation,numberOfCandidates, & - & candidateElements,candidateElementLines,closestElements,closestElementLines,closestDistances,err,error,*) + & candidateElements,candidateElementLines,startingXi,closestElements,closestElementLines,closestDistances,err,error,*) !Argument variables TYPE(DataProjectionType), POINTER :: dataProjection !=1 and <= 3." CALL FlagError(localError,err,error,*999) END SELECT - ALLOCATE(dataProjection%startingXi(dataProjection%numberOfXi),STAT=err) + ALLOCATE(dataProjection%startingXi(dataProjection%datapoints%numberOfDatapoints,dataProjection%numberOfXi),STAT=err) IF(err/=0) CALL FlagError("Could not allocate data points data projection starting xi.",err,error,*999) dataProjection%startingXi=0.5_DP !0) THEN + IF(numberOfClosestCandidates>0) THEN + projectedXi(:,dataPointIdx) = dataProjection%startingXi(dataPointIdx,:) CALL DataProjection_NewtonLinesEvaluate(dataProjection,interpolatedPoint, & & dataPoints%dataPoints(dataPointIdx)%position,closestElements( & & dataPointIdx,1:numberOfClosestCandidates),closestLinesFaces(dataPointIdx,1: & @@ -1874,7 +1884,8 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection !Newton project to closest faces, and find miminum projection DO dataPointIdx=1,numberOfDataPoints numberOfClosestCandidates=globalToLocalNumberOfClosestCandidates(dataPointIdx) - IF(numberOfClosestCandidates>0) THEN + IF(numberOfClosestCandidates>0) THEN + projectedXi(:,dataPointIdx) = dataProjection%startingXi(dataPointIdx,:) CALL DataProjection_NewtonFacesEvaluate(dataProjection,interpolatedPoint, & & dataPoints%dataPoints(dataPointIdx)%position,closestElements( & & dataPointIdx,1:numberOfClosestCandidates),closestLinesFaces(dataPointIdx, & @@ -1891,7 +1902,8 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CASE(1) !1D element DO dataPointIdx=1,numberOfDataPoints numberOfClosestCandidates=globalToLocalNumberOfClosestCandidates(dataPointIdx) - IF(numberOfClosestCandidates>0) THEN + IF(numberOfClosestCandidates>0) THEN + projectedXi(:,dataPointIdx) = dataProjection%startingXi(dataPointIdx,:) CALL DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx, & & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx), & @@ -1904,7 +1916,8 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CASE(2) !2D element DO dataPointIdx=1,numberOfDataPoints numberOfClosestCandidates=globalToLocalNumberOfClosestCandidates(dataPointIdx) - IF(numberOfClosestCandidates>0) THEN + IF(numberOfClosestCandidates>0) THEN + projectedXi(:,dataPointIdx) = dataProjection%startingXi(dataPointIdx,:) CALL DataProjection_NewtonElementsEvaluate_2(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx, & & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx), & @@ -1918,7 +1931,8 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CASE(3) !3D element DO dataPointIdx=1,numberOfDataPoints numberOfClosestCandidates=globalToLocalNumberOfClosestCandidates(dataPointIdx) - IF(numberOfClosestCandidates>0) THEN + IF(numberOfClosestCandidates>0) THEN + projectedXi(:,dataPointIdx) = dataProjection%startingXi(dataPointIdx,:) CALL DataProjection_NewtonElementsEvaluate_3(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx, & & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx), & @@ -2033,6 +2047,9 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CASE(DATA_PROJECTION_BOUNDARY_FACES_PROJECTION_TYPE) !Newton project to closest faces, and find miminum projection DO dataPointIdx=1,numberOfDataPoints + !IF (dataPointIdx==132) THEN + ! WRITE(*,*) dataPointIdx + !ENDIF CALL DataProjection_NewtonFacesEvaluate(dataProjection,interpolatedPoint,dataPoints%dataPoints( & & dataPointIdx)%position,closestElements(dataPointIdx,:),closestLinesFaces(dataPointIdx,:), & & dataProjection%dataProjectionResults(dataPointIdx)%exitTag,dataProjection% & @@ -2366,7 +2383,7 @@ SUBROUTINE DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPo INTEGER(INTG), INTENT(OUT) :: projectionExitTag !SIZE(dataProjection%startingXi,1)) THEN - startingXi(1:SIZE(dataProjection%startingXi,1))=dataProjection%startingXi(1:SIZE(dataProjection%startingXi,1)) - startingXi(SIZE(dataProjection%startingXi,1):dataProjection%numberOfXi)=0.5_DP + startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:SIZE(dataProjection%startingXi,1))= & + & dataProjection%startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:SIZE(dataProjection%startingXi,1)) + startingXi(1:dataProjection%datapoints%numberOfDatapoints,SIZE(dataProjection%startingXi,1):dataProjection%numberOfXi)= & + & 0.5_DP ELSE - startingXi(1:SIZE(dataProjection%startingXi,1))=dataProjection%startingXi(1:dataProjection%numberOfXi) + startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:SIZE(dataProjection%startingXi,1))= & + & dataProjection%startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:dataProjection%numberOfXi) ENDIF IF(ALLOCATED(dataProjection%startingXi)) DEALLOCATE(dataProjection%startingXi) - ALLOCATE(dataProjection%startingXi(dataProjection%numberOfXi),STAT=err) + ALLOCATE(dataProjection%startingXi(1:dataProjection%datapoints%numberOfDatapoints,dataProjection%numberOfXi),STAT=err) IF(err/=0) CALL FlagError("Could not allocate data projection starting xi.",err,error,*999) - dataProjection%startingXi(1:dataProjection%numberOfXi)=startingXi(1:dataProjection%numberOfXi) + dataProjection%startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:dataProjection%numberOfXi)= & + & startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:dataProjection%numberOfXi) DEALLOCATE(startingXi) ENDIF @@ -4578,7 +4599,7 @@ SUBROUTINE DataProjection_StartingXiGet(dataProjection,startingXi,err,error,*) !Argument variables TYPE(DataProjectionType), POINTER :: dataProjection != "//TRIM(NumberToVString(SIZE(dataProjection%startingXi,1),"*",err,error))//"." CALL FlagError(localError,err,error,*999) ENDIF - startingXi(1:SIZE(dataProjection%startingXi,1))=dataProjection%startingXi(1:SIZE(dataProjection%startingXi,1)) + startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:SIZE(dataProjection%startingXi,1))=& + & dataProjection%startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:SIZE(dataProjection%startingXi,1)) EXITS("DataProjection_StartingXiGet") RETURN @@ -4612,11 +4634,11 @@ SUBROUTINE DataProjection_StartingXiSet(dataProjection,startingXi,err,error,*) !Argument variables TYPE(DataProjectionType), POINTER :: dataProjection !1.0_DP)) THEN - validInput=.FALSE. - EXIT !this do - ENDIF - ENDDO !xiIdx + DO datapoint=1,dataProjection%datapoints%numberOfDatapoints + DO xiIdx=1,dataProjection%numberOfXi + IF((startingXi(datapoint,xiIdx)<0.0_DP).OR.(startingXi(datapoint,xiIdx)>1.0_DP)) THEN + validInput=.FALSE. + EXIT !this do + ENDIF + ENDDO !xiIdx + ENDDO !datapoint IF(.NOT.validInput) CALL FlagError("Data projection starting xi must be between 0.0 and 1.0.",err,error,*999) - dataProjection%startingXi(1:SIZE(startingXi))=startingXi(1:SIZE(startingXi)) + dataProjection%startingXi(1:SIZE(startingXi,1), 1:SIZE(startingXi,2))=startingXi EXITS("DataProjection_StartingXiSet") RETURN -999 ERRORSEXITS("DataProjection_StartingXiSet",err,error) +999 ERRORSEXITS("DataProjection_StartingXiSet",err,error) RETURN 1 END SUBROUTINE DataProjection_StartingXiSet diff --git a/src/opencmiss_iron.F90 b/src/opencmiss_iron.F90 index 510774e2..deae8cb4 100644 --- a/src/opencmiss_iron.F90 +++ b/src/opencmiss_iron.F90 @@ -24356,7 +24356,7 @@ SUBROUTINE cmfe_DataProjection_StartingXiGetNumber(regionUserNumber,dataPointsUs INTEGER(INTG), INTENT(IN) :: regionUserNumber ! Date: Sun, 11 Aug 2019 01:19:55 +1200 Subject: [PATCH 15/17] Adding EquationsSetSubtypes.REFERENCE_STATE_MOONEY_RIVLIN back to finite elasticity routines #183 --- src/finite_elasticity_routines.F90 | 70 ++++++++++++++++++++++-------- 1 file changed, 53 insertions(+), 17 deletions(-) diff --git a/src/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90 index c5b8c450..3b8ff4a0 100644 --- a/src/finite_elasticity_routines.F90 +++ b/src/finite_elasticity_routines.F90 @@ -732,7 +732,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT SELECT CASE(EQUATIONS_SET%specification(3)) CASE(EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, & - & EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, & + & EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE) LOCAL_ERROR="Analytic Jacobian has not been validated for the Mooney-Rivlin equations, please use finite differences instead." CALL FlagWarning(LOCAL_ERROR,ERR,ERROR,*999) @@ -929,7 +929,8 @@ SUBROUTINE FiniteElasticity_FiniteElementJacobianEvaluate(EQUATIONS_SET,ELEMENT_ nonlinearMatrices=>vectorMatrices%nonlinearMatrices jacobianMatrix=>nonlinearMatrices%jacobians(1)%ptr IF(jacobianMatrix%updateJacobian) THEN - IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN DEPENDENT_FIELD=>equations%interpolation%geometricField GEOMETRIC_FIELD=>equations%interpolation%dependentField ELSE @@ -980,7 +981,8 @@ SUBROUTINE FiniteElasticity_FiniteElementJacobianEvaluate(EQUATIONS_SET,ELEMENT_ !Point interpolation pointer geometricInterpPoint=>equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr geometricInterpPointMetrics=>equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr - IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN dependentInterpPoint=>equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr dependentInterpPointMetrics=>equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr geometricInterpPoint=>equations%interpolation%dependentInterpPoint(FIELD_VAR_TYPE)%ptr @@ -1513,7 +1515,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN rhsVector=>vectorMatrices%rhsVector vectorMapping =>vectorEquations%vectorMapping - IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN DEPENDENT_FIELD =>equations%interpolation%geometricField GEOMETRIC_FIELD =>equations%interpolation%dependentField ELSE @@ -1570,7 +1573,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN !Grab interpolation parameters FIELD_VARIABLE=>EQUATIONS_SET%equations%vectorEquations%vectorMapping%nonlinearMapping%residualVariables(1)%ptr FIELD_VAR_TYPE=FIELD_VARIABLE%VARIABLE_TYPE - IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN GEOMETRIC_INTERPOLATION_PARAMETERS=>equations%interpolation%dependentInterpParameters(FIELD_VAR_TYPE)%ptr DEPENDENT_INTERPOLATION_PARAMETERS=>equations%interpolation%geometricInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr ELSE @@ -1631,7 +1635,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN ! END IF !Point interpolation pointer - IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN DEPENDENT_INTERPOLATED_POINT=>equations%interpolation%geometricInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr DEPENDENT_INTERPOLATED_POINT_METRICS=>equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr GEOMETRIC_INTERPOLATED_POINT=>equations%interpolation%dependentInterpPoint(FIELD_VAR_TYPE)%ptr @@ -1667,7 +1672,7 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN SELECT CASE(EQUATIONS_SET_SUBTYPE) ! --------------------------------------------------------------- CASE(EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, & - & EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE) + & EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) !Loop over gauss points and add residuals DO gauss_idx=1,DEPENDENT_NUMBER_OF_GAUSS_POINTS !Interpolate dependent, geometric, fibre and materials fields @@ -2138,7 +2143,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN !Hydrostatic pressure component (skip for membrane problems) IF (EQUATIONS_SET_SUBTYPE /= EQUATIONS_SET_MEMBRANE_SUBTYPE) THEN - IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN HYDROSTATIC_PRESSURE_COMPONENT=GEOMETRIC_FIELD%VARIABLES(var1)%NUMBER_OF_COMPONENTS DEPENDENT_COMPONENT_INTERPOLATION_TYPE=GEOMETRIC_FIELD%VARIABLES(var1)%COMPONENTS( & & HYDROSTATIC_PRESSURE_COMPONENT)%INTERPOLATION_TYPE @@ -2153,7 +2159,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN TEMPTERM1=GAUSS_WEIGHT*Jxxi*(Jznu-Jg) ENDIF IF(DEPENDENT_COMPONENT_INTERPOLATION_TYPE==FIELD_NODE_BASED_INTERPOLATION) THEN !node based - IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN COMPONENT_BASIS=>GEOMETRIC_FIELD%VARIABLES(var1)%COMPONENTS(HYDROSTATIC_PRESSURE_COMPONENT)%DOMAIN% & & TOPOLOGY%ELEMENTS%ELEMENTS(elementNumber)%BASIS ELSE @@ -2524,7 +2531,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN !Hydrostatic pressure component (skip for membrane problems) IF (EQUATIONS_SET_SUBTYPE /= EQUATIONS_SET_MEMBRANE_SUBTYPE) THEN - IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN HYDROSTATIC_PRESSURE_COMPONENT=GEOMETRIC_FIELD%VARIABLES(var1)%NUMBER_OF_COMPONENTS DEPENDENT_COMPONENT_INTERPOLATION_TYPE=GEOMETRIC_FIELD%VARIABLES(var1)%COMPONENTS( & & HYDROSTATIC_PRESSURE_COMPONENT)%INTERPOLATION_TYPE @@ -2539,7 +2547,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN TEMPTERM1=GAUSS_WEIGHT*Jxxi*(Je-1.0_DP) ENDIF IF(DEPENDENT_COMPONENT_INTERPOLATION_TYPE==FIELD_NODE_BASED_INTERPOLATION) THEN !node based - IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN COMPONENT_BASIS=>GEOMETRIC_FIELD%VARIABLES(var1)%COMPONENTS(HYDROSTATIC_PRESSURE_COMPONENT)%DOMAIN% & & TOPOLOGY%ELEMENTS%ELEMENTS(elementNumber)%BASIS ELSE @@ -3555,6 +3564,7 @@ SUBROUTINE FiniteElasticity_FiniteElementPreResidualEvaluate(equationsSet,err,er & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_STATIC_INRIA_SUBTYPE, & & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_HOLMES_MOW_SUBTYPE, & & EQUATIONS_SET_ELASTICITY_FLUID_PRES_HOLMES_MOW_ACTIVE_SUBTYPE, & @@ -3623,6 +3633,7 @@ SUBROUTINE FiniteElasticity_FiniteElementPostResidualEvaluate(EQUATIONS_SET,err, & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE,& + & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_STATIC_INRIA_SUBTYPE, & & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_HUMPHREY_YIN_SUBTYPE,& & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_HOLMES_MOW_SUBTYPE, & @@ -5032,7 +5043,8 @@ SUBROUTINE FiniteElasticity_SurfacePressureResidualEvaluate(EQUATIONS_SET,ELEMEN nonlinearMatrices=>vectorEquations%vectorMatrices%nonlinearMatrices EQUATIONS_SET_SUBTYPE = EQUATIONS_SET%SPECIFICATION(3) - IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN DEPENDENT_FIELD=>equations%interpolation%geometricField ELSE DEPENDENT_FIELD=>equations%interpolation%dependentField @@ -5071,7 +5083,8 @@ SUBROUTINE FiniteElasticity_SurfacePressureResidualEvaluate(EQUATIONS_SET,ELEMEN DEPENDENT_FACE_BASIS=>DECOMPOSITION%DOMAIN(MESH_COMPONENT_NUMBER)%ptr%TOPOLOGY%FACES%FACES(face_number)%BASIS FACE_QUADRATURE_SCHEME=>DEPENDENT_FACE_BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr - IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN + IF (EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE == EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN FACE_DEPENDENT_INTERPOLATION_PARAMETERS=>equations%interpolation%geometricInterpParameters(FIELD_VAR_U_TYPE)%ptr FACE_DEPENDENT_INTERPOLATED_POINT=>equations%interpolation%geometricInterpPoint(FIELD_VAR_U_TYPE)%ptr FACE_DEPENDENT_INTERPOLATED_POINT_METRICS=>equations%interpolation% & @@ -6038,7 +6051,17 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO & EQUATIONS_SET_INCOMPRESSIBLE_FINITE_ELASTICITY_DARCY_SUBTYPE,EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE, & & EQUATIONS_SET_1D3D_MONODOMAIN_ELASTICITY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_POLYNOMIAL_SUBTYPE, & & EQUATIONS_SET_MONODOMAIN_ELASTICITY_W_TITIN_SUBTYPE,EQUATIONS_SET_MONODOMAIN_ELASTICITY_VELOCITY_SUBTYPE, & - & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE,EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE) + & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_ACTIVE_SUBTYPE,EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) + + IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN + PRESSURE_COMPONENT=GEOMETRIC_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS + P=GEOMETRIC_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV) + ELSE + PRESSURE_COMPONENT=DEPENDENT_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS + P=DEPENDENT_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV) + ENDIF + !Form of constitutive model is: ! W=c1*(I1-3)+c2*(I2-3)+p*(I3-1) !Also assumed I3 = det(AZL) = 1.0 @@ -7022,10 +7045,17 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO SELECT CASE(EQUATIONS_SET%specification(3)) CASE(EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, & - & EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, & + & EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_MR_AND_GROWTH_LAW_IN_CELLML_SUBTYPE) - PRESSURE_COMPONENT=DEPENDENT_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS - P=DEPENDENT_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV) + + IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE) THEN + PRESSURE_COMPONENT=GEOMETRIC_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS + P=GEOMETRIC_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV) + ELSE + PRESSURE_COMPONENT=DEPENDENT_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS + P=DEPENDENT_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV) + ENDIF + !Form of constitutive model is: !W=c1*(I1-3)+c2*(I2-3)+p/2*(I3-1) @@ -7397,6 +7427,7 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_LAW_IN_CELLML_EVALUATE_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, & & EQUATIONS_SET_GROWTH_LAW_IN_CELLML_SUBTYPE, & @@ -7481,6 +7512,7 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET SELECT CASE(EQUATIONS_SET_SUBTYPE) CASE(EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE, & & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, & & EQUATIONS_SET_COMPRESSIBLE_FINITE_ELASTICITY_SUBTYPE,& @@ -7591,6 +7623,7 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET & EQUATIONS_SET_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_NO_SUBTYPE,EQUATIONS_SET_MEMBRANE_SUBTYPE, & & EQUATIONS_SET_ORTHOTROPIC_MATERIAL_HOLZAPFEL_OGDEN_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_HUMPHREY_YIN_SUBTYPE, & & EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE,EQUATIONS_SET_1D3D_MONODOMAIN_ELASTICITY_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & @@ -10167,6 +10200,7 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET NUMBER_OF_FLUID_COMPONENTS=0 SELECT CASE(EQUATIONS_SET_SUBTYPE) CASE(EQUATIONS_SET_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NO_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_FINITE_ELASTICITY_DARCY_SUBTYPE, & & EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE, & @@ -10820,6 +10854,7 @@ SUBROUTINE FiniteElasticity_EquationsSetSolutionMethodSet(EQUATIONS_SET,SOLUTION & EQUATIONS_SET_INCOMPRESSIBLE_ELASTICITY_DRIVEN_DARCY_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_LAW_IN_CELLML_EVALUATE_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, & & EQUATIONS_SET_GROWTH_LAW_IN_CELLML_SUBTYPE, & @@ -10914,6 +10949,7 @@ SUBROUTINE FiniteElasticity_EquationsSetSpecificationSet(equationsSet,specificat & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_LAW_IN_CELLML_EVALUATE_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, & & EQUATIONS_SET_GROWTH_LAW_IN_CELLML_SUBTYPE, & From d14e9149353db3d5a20d0b83c455536bcb4d2aa5 Mon Sep 17 00:00:00 2001 From: Prasad Date: Tue, 26 Nov 2019 10:13:22 +1300 Subject: [PATCH 16/17] Active contraction (latest updates from Mario) #163 --- src/equations_set_constants.F90 | 19 +- src/equations_set_routines.F90 | 8 + src/finite_elasticity_routines.F90 | 486 +++++++++++++++++++++++++---- src/opencmiss_iron.F90 | 18 +- 4 files changed, 462 insertions(+), 69 deletions(-) diff --git a/src/equations_set_constants.F90 b/src/equations_set_constants.F90 index 53a9bbac..e009e78c 100644 --- a/src/equations_set_constants.F90 +++ b/src/equations_set_constants.F90 @@ -154,6 +154,8 @@ MODULE EquationsSetConstants INTEGER(INTG), PARAMETER :: EQUATIONS_SET_ANISOTROPIC_POLYNOMIAL_ACTIVE_SUBTYPE=30 INTEGER(INTG), PARAMETER :: EQUATIONS_SET_HOLZAPFEL_OGDEN_ACTIVECONTRACTION_SUBTYPE=31 INTEGER(INTG), PARAMETER :: EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE=32 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE=44 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE=45 INTEGER(INTG), PARAMETER :: EQUATIONS_SET_CONSTITUTIVE_AND_GROWTH_LAW_IN_CELLML_SUBTYPE=33 INTEGER(INTG), PARAMETER :: EQUATIONS_SET_GROWTH_LAW_IN_CELLML_SUBTYPE=34 INTEGER(INTG), PARAMETER :: EQUATIONS_SET_ACTIVE_STRAIN_SUBTYPE=35 @@ -406,14 +408,17 @@ MODULE EquationsSetConstants !> \see EquationsSetConstants !>@{ INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DEFORMATION_GRADIENT_TENSOR=1 !@} - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_NUMBER_OF_TENSOR_TYPES=7 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_NUMBER_OF_TENSOR_TYPES=10 !> \addtogroup EquationsSetConstants_DynamicMatrixTypes EquationsSetConstants::DynamicMatrixTypes !> \brief Type of matrix in a dynamic equations set diff --git a/src/equations_set_routines.F90 b/src/equations_set_routines.F90 index 10437fa8..9fd95202 100644 --- a/src/equations_set_routines.F90 +++ b/src/equations_set_routines.F90 @@ -2772,6 +2772,8 @@ SUBROUTINE EquationsSet_FiniteElementJacobianEvaluate(equationsSet,elementNumber CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) CALL ClassicalField_FiniteElementJacobianEvaluate(equationsSet,elementNumber,err,error,*999) + CASE(EQUATIONS_SET_FITTING_CLASS) + CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_MODAL_CLASS) @@ -2907,6 +2909,10 @@ SUBROUTINE EquationsSet_FiniteElementJacobianEvaluateFD(equationsSet,elementNumb ! so let's just loop over component, node/el, derivative column=0 ! element jacobian matrix column number DO componentIdx=1,columnVariable%NUMBER_OF_COMPONENTS + ! adjusts the pertubation for the hydrostatic pressure to not use the L2Norm + IF(componentIdx==4) THEN + delta=(1.0_DP)*equationsSet%equations%jacobianFiniteDifferenceStepSize + ENDIF elementsTopology=>columnVariable%COMPONENTS(componentIdx)%DOMAIN%TOPOLOGY%ELEMENTS componentInterpolationType=columnVariable%COMPONENTS(componentIdx)%INTERPOLATION_TYPE SELECT CASE(componentInterpolationType) @@ -3012,6 +3018,8 @@ SUBROUTINE EquationsSet_FiniteElementResidualEvaluate(equationsSet,elementNumber CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) CALL ClassicalField_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*999) + CASE(EQUATIONS_SET_FITTING_CLASS) + CALL Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*999) CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_MODAL_CLASS) diff --git a/src/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90 index 10e3e3d3..d3e48a97 100644 --- a/src/finite_elasticity_routines.F90 +++ b/src/finite_elasticity_routines.F90 @@ -794,6 +794,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)+P CASE(EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE,EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) PRESSURE_COMPONENT=DEPENDENT_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS P=DEPENDENT_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV) @@ -803,30 +804,7 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR(EQUATIONS_SET,DEPENDENT_INT TEMPTERM1=0.5_DP*C(1)*EXP(0.5_DP*DOT_PRODUCT(E,DQ_DE)) !Calculate 2nd Piola tensor (in Voigt form) STRESS_TENSOR=TEMPTERM1*DQ_DE + P*AZUv - !lambda = 1.0_DP - !lambda(1) = SQRT(AZL(1,1)) - IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE) THEN - !add active contraction stress values - CALL Field_VariableGet(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE,FIELD_VARIABLE,ERR,ERROR,*999) - DO component_idx=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS - SELECT CASE(FIELD_VARIABLE%COMPONENTS(component_idx)%INTERPOLATION_TYPE) - CASE(FIELD_CONSTANT_INTERPOLATION) - VALUE=0.1_DP - CASE(FIELD_GAUSS_POINT_BASED_INTERPOLATION) - dof_idx=FIELD_VARIABLE%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP% & - & GAUSS_POINTS(GAUSS_POINT_NUMBER,ELEMENT_NUMBER) - CALL FIELD_PARAMETER_SET_GET_LOCAL_DOF(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & - & FIELD_VALUES_SET_TYPE,dof_idx,VALUE,err,error,*999) - CASE(FIELD_ELEMENT_BASED_INTERPOLATION) - VALUE=INDEPENDENT_INTERPOLATED_POINT%VALUES(component_idx,NO_PART_DERIV) - CASE DEFAULT - LOCAL_ERROR="This independent field variable interpolation type is not supported." - CALL FlagError(LOCAL_ERROR,err,error,*999) - END SELECT - STRESS_TENSOR(component_idx)=STRESS_TENSOR(component_idx)+VALUE - !STRESS_TENSOR(component_idx)=STRESS_TENSOR(component_idx)+VALUE*(1.0_DP+1.45_DP*(lambda(component_idx)-1.0_DP)) - ENDDO - ENDIF + !\todo blas has routines specifically for symmetric matrices, so it would be worth to check if these could give some speedup. @@ -1433,8 +1411,9 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN INTEGER(INTG) :: var2 ! Variable number corresponding to 'DELUDLEN' in single physics case INTEGER(INTG), POINTER :: EQUATIONS_SET_FIELD_DATA(:) REAL(DP) :: DZDNU(3,3),DZDNUT(3,3),dzdx(3,3),AZL(3,3),AZU(3,3),Fe(3,3),FeT(3,3),Fg(3,3),C(3,3),f(3,3),E(3,3),I3,P, & - & piolaTensor(3,3),TEMP(3,3),prevdzdx(3,3),prevdZdNu(3,3),invPrevdZdNu(3,3) - REAL(DP) :: cauchyTensor(3,3),JGW_CAUCHY_TENSOR(3,3),kirchoffTensor(3,3),STRESS_TENSOR(6),growthValues(3) + & piolaTensor(3,3),TEMP(3,3),prevdzdx(3,3),prevdZdNu(3,3),invPrevdZdNu(3,3), DNUODNU(3,3) + REAL(DP) :: cauchyTensor(3,3),JGW_CAUCHY_TENSOR(3,3),kirchoffTensor(3,3),STRESS_TENSOR(6),growthValues(3), & + & cauchyTensorFbr(3,3) REAL(DP) :: deformationGradientTensor(3,3),growthTensor(3,3),growthTensorInverse(3,3),growthTensorInverseTranspose(3,3), & & fibreGrowth,sheetGrowth,normalGrowth,fibreVector(3),sheetVector(3),normalVector(3) REAL(DP) :: dNudXi(3,3),dXidNu(3,3) @@ -1594,7 +1573,9 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN IF(DARCY_DEPENDENT) THEN DARCY_DEPENDENT_INTERPOLATION_PARAMETERS=>equations%interpolation%dependentInterpParameters(FIELD_V_VARIABLE_TYPE)%ptr ELSE IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE .OR. & - & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE) THEN + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE) THEN INDEPENDENT_INTERPOLATION_PARAMETERS=>equations%interpolation%independentInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr ENDIF ! IF(ASSOCIATED(SOURCE_FIELD)) THEN @@ -1625,7 +1606,9 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,elementNumber, & & DARCY_DEPENDENT_INTERPOLATION_PARAMETERS,err,error,*999) ELSE IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE .OR. & - & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE) THEN + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE) THEN CALL FIELD_INTERPOLATION_PARAMETERS_ELEMENT_GET(FIELD_VALUES_SET_TYPE,elementNumber, & & INDEPENDENT_INTERPOLATION_PARAMETERS,err,error,*999) ENDIF @@ -1661,7 +1644,9 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN IF(DARCY_DEPENDENT) THEN DARCY_DEPENDENT_INTERPOLATED_POINT=>equations%interpolation%dependentInterpPoint(FIELD_V_VARIABLE_TYPE)%ptr ELSE IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE .OR. & - & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE) THEN + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE) THEN INDEPENDENT_INTERPOLATED_POINT=>equations%interpolation%independentInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr ENDIF IF(ASSOCIATED(SOURCE_FIELD)) THEN @@ -1831,6 +1816,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN & EQUATIONS_SET_ISOTROPIC_EXPONENTIAL_SUBTYPE, & & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_EXPONENTIAL_SUBTYPE, & & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE,EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_ORTHOTROPIC_MATERIAL_COSTA_SUBTYPE, EQUATIONS_SET_ACTIVECONTRACTION_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_FINITE_ELASTICITY_DARCY_SUBTYPE,EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE, & @@ -1868,7 +1855,9 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,gauss_idx, & & DARCY_DEPENDENT_INTERPOLATED_POINT,err,error,*999) ELSE IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE .OR. & - & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE) THEN + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE .OR. & + & EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE) THEN CALL FIELD_INTERPOLATE_GAUSS(NO_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,gauss_idx, & & INDEPENDENT_INTERPOLATED_POINT,err,error,*999) ENDIF @@ -1902,7 +1891,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN !Calculate Sigma=1/Jznu.FTF', the Cauchy stress tensor at the gauss point CALL FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPOLATED_POINT, & & MATERIALS_INTERPOLATED_POINT,GEOMETRIC_INTERPOLATED_POINT,DARCY_DEPENDENT_INTERPOLATED_POINT, & - & INDEPENDENT_INTERPOLATED_POINT,cauchyTensor,Jznu,DZDNU,elementNumber,gauss_idx,err,error,*999) + & INDEPENDENT_INTERPOLATED_POINT,cauchyTensor,cauchyTensorFbr,Jznu,DZDNU,DZDX,DNUODNU, & + & elementNumber,gauss_idx,err,error,*999) IF(DIAGNOSTICS1) THEN CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"",err,error,*999) @@ -3230,7 +3220,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN !Calculate Sigma=1/Jznu.FTF', the Cauchy stress tensor at the gauss point CALL FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPOLATED_POINT, & & MATERIALS_INTERPOLATED_POINT,GEOMETRIC_INTERPOLATED_POINT,DARCY_DEPENDENT_INTERPOLATED_POINT, & - & INDEPENDENT_INTERPOLATED_POINT,cauchyTensor,Jznu,DZDNU,elementNumber,gauss_idx,err,error,*999) + & INDEPENDENT_INTERPOLATED_POINT,cauchyTensor,cauchyTensorFbr,Jznu,DZDNU,DZDX,DNUODNU, & + & elementNumber,gauss_idx,err,error,*999) !Calculate dPhi/dZ at the gauss point, Phi is the basis function CALL FINITE_ELASTICITY_GAUSS_DFDZ(DEPENDENT_INTERPOLATED_POINT,elementNumber,gauss_idx,numberOfDimensions, & @@ -3341,7 +3332,8 @@ SUBROUTINE FiniteElasticity_FiniteElementResidualEvaluate(EQUATIONS_SET,elementN !Calculate Cauchy stress tensor at the gauss point CALL FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPOLATED_POINT, & & MATERIALS_INTERPOLATED_POINT,GEOMETRIC_INTERPOLATED_POINT,DARCY_DEPENDENT_INTERPOLATED_POINT, & - & INDEPENDENT_INTERPOLATED_POINT,cauchyTensor,Jznu,DZDNU,elementNumber,gauss_idx,err,error,*999) + & INDEPENDENT_INTERPOLATED_POINT,cauchyTensor,cauchyTensorFbr,Jznu,DZDNU,DZDX,DNUODNU, & + & elementNumber,gauss_idx,err,error,*999) !Calculate dF/DZ at the gauss point CALL FINITE_ELASTICITY_GAUSS_DFDZ(DEPENDENT_INTERPOLATED_POINT,elementNumber,gauss_idx,numberOfDimensions, & @@ -3563,8 +3555,10 @@ SUBROUTINE FiniteElasticity_FiniteElementPreResidualEvaluate(equationsSet,err,er & EQUATIONS_SET_ORTHOTROPIC_MATERIAL_HOLZAPFEL_OGDEN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & - & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_STATIC_INRIA_SUBTYPE, & & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_HOLMES_MOW_SUBTYPE, & & EQUATIONS_SET_ELASTICITY_FLUID_PRES_HOLMES_MOW_ACTIVE_SUBTYPE, & @@ -3632,8 +3626,10 @@ SUBROUTINE FiniteElasticity_FiniteElementPostResidualEvaluate(EQUATIONS_SET,err, & EQUATIONS_SET_ORTHOTROPIC_MATERIAL_HOLZAPFEL_OGDEN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & - & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE,& + & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE,& & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_STATIC_INRIA_SUBTYPE, & & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_HUMPHREY_YIN_SUBTYPE,& & EQUATIONS_SET_ELASTICITY_FLUID_PRESSURE_HOLMES_MOW_SUBTYPE, & @@ -4386,8 +4382,8 @@ SUBROUTINE FiniteElasticity_StressStrainPoint(equationsSet,evaluateType,numberOf TYPE(VARYING_STRING), INTENT(OUT) :: error !Evaluates the Cauchy stress tensor at a given Gauss point SUBROUTINE FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPOLATED_POINT, & & MATERIALS_INTERPOLATED_POINT,GEOMETRIC_INTERPOLATED_POINT,DARCY_DEPENDENT_INTERPOLATED_POINT, & - & INDEPENDENT_INTERPOLATED_POINT,CAUCHY_TENSOR,Jznu,DZDNU,ELEMENT_NUMBER,GAUSS_POINT_NUMBER,err,error,*) + & INDEPENDENT_INTERPOLATED_POINT,CAUCHY_TENSOR,CAUCHY_TENSOR_FBR,Jznu,DZDNU,DZDX,DNUODNU, & + & ELEMENT_NUMBER,GAUSS_POINT_NUMBER,err,error,*) !Argument variables TYPE(EQUATIONS_SET_TYPE), POINTER, INTENT(IN) :: EQUATIONS_SET !EQUATIONS_SET% & + & EQUATIONS%INTERPOLATION%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr + ! deformed metric + NULLIFY(DEPENDENT_INTERPOLATED_POINT_METRICS) + DEPENDENT_INTERPOLATED_POINT_METRICS=>EQUATIONS_SET% & + & EQUATIONS%INTERPOLATION%dependentInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr + XR=GEOMETRIC_INTERPOLATED_POINT_METRICS%DX_DXI !reference spatial directions + XD=DEPENDENT_INTERPOLATED_POINT_METRICS%DX_DXI !deformed spatial directions + ! fibre interpolation + NULLIFY(FIBRE_INTERPOLATED_POINT) + FIBRE_INTERPOLATED_POINT=>EQUATIONS_SET%EQUATIONS% & + & INTERPOLATION%fibreInterpPoint(FIELD_U_VARIABLE_TYPE)%ptr + numberOfXDimensions=GEOMETRIC_INTERPOLATED_POINT_METRICS%NUMBER_OF_X_DIMENSIONS + numberOfXiDimensions=GEOMETRIC_INTERPOLATED_POINT_METRICS%NUMBER_OF_XI_DIMENSIONS + numberOfNuDimensions=SIZE(FIBRE_INTERPOLATED_POINT%VALUES,1) + ! call this function to get the transformation tensors between material and spatial coordinates + CALL IdentityMatrix(DXDNU,err,error,*999) ! initialise rotation matrix + CALL Coordinates_MaterialSystemCalculate(GEOMETRIC_INTERPOLATED_POINT_METRICS,FIBRE_INTERPOLATED_POINT, & + & DNUDX,DXDNU,DNUDXI(1:numberOfXDimensions,1:numberOfXiDimensions), & + & DXIDNU(1:numberOfXiDimensions,1:numberOfXDimensions),err,error,*999) + DO component_idx=1,numberOfXDimensions + ! normalise reference and deformed spatial vectors + CALL Normalise(XR(1:numberOfXDimensions,component_idx), & + & XR(1:numberOfXDimensions,component_idx),err,error,*999) + CALL Normalise(XD(1:numberOfXDimensions,component_idx), & + & XD(1:numberOfXDimensions,component_idx),err,error,*999) + ENDDO + ! calculate deformed material coordinates + ! dnu/dxi = dz/dNu * dX/dxi + CALL MatrixVectorProduct(DZDNU,XR(:,1),fibre_def,err,error,*999) + CALL MatrixVectorProduct(DZDNU,XR(:,2),sheet_def,err,error,*999) + CALL Normalise(fibre_def,fibre_def,err,error,*999) + CALL Normalise(sheet_def,sheet_def,err,error,*999) + ! calculate normal vector on the plane spanned by the fibre and sheet vector + CALL CrossProduct(fibre_def,sheet_def,normal_def,err,error,*999) ! orthogonal material coorindate 3 + ! calculate orthogonal sheet vector + CALL CrossProduct(fibre_def,normal_def,sheet_def,err,error,*999) ! orthogonal material coorindate 2 + CALL IdentityMatrix(NUDO,err,error,*999) ! initialise deformed fibre vector matrix + NUDO(1:numberOfXDimensions,1)=fibre_def + NUDO(1:numberOfXDimensions,2)=sheet_def + NUDO(1:numberOfXDimensions,3)=normal_def + ! normalise the deformed orthogonal material vectors + CALL Normalise(NUDO,NUDO,err,error,*999) + DO component_idx=1,numberOfXDimensions + CALL Normalise(NUDO(1:numberOfXDimensions,component_idx), & + & NUDO(1:numberOfXDimensions,component_idx),err,error,*999) + ENDDO + + ! construct the tensor to transform from geometric to orthogonal material coordinates + ! we have dnuo/dxi and dxi/dz + ! q^-1 = dnuo/dz = dnuo/dxi * dxi/dz + CALL IdentityMatrix(XID,err,error,*999) ! initialise dxi/dz + XID=DEPENDENT_INTERPOLATED_POINT_METRICS%DXI_DX + CALL Normalise(XID,XID,err,error,*999) + DO component_idx=1,numberOfXDimensions + CALL Normalise(XID(1:numberOfXDimensions,component_idx), & + & XID(1:numberOfXDimensions,component_idx),err,error,*999) + ENDDO + + CALL IdentityMatrix(DNUODZ,err,error,*999) ! initialise q^-1 + CALL MatrixProduct(NUDO(1:numberOfXDimensions,1:numberOfXiDimensions), & + & XID(1:numberOfXDimensions,1:numberOfXiDimensions), & + & DNUODZ(1:numberOfXDimensions,1:numberOfXDimensions),err,error,*999) + ! Calculate q + CALL IdentityMatrix(DZDNUO,err,error,*999) ! initialise q + CALL Invert(DNUODZ,DZDNUO,DET_DNUODZ,err,error,*999) + ! Calculate q^-T + CALL IdentityMatrix(DNUODZT,err,error,*999) ! initialise q^-T + CALL MatrixTranspose(DNUODZ,DNUODZT,err,error,*999) + ! Calculate q^T + CALL IdentityMatrix(DZDNUOT,err,error,*999) ! initialise q^T + CALL MatrixTranspose(DZDNUO,DZDNUOT,err,error,*999) + + ! Add active stress + ! Hunter-McCulloch-ter Keurs constitutive model + ! T_a = T_Ca * [1 + beta * (lambda - 1)] + ! our CAUCHY_TENSOR_DEFGEO equals to T_Ca + ! lamda(stretch) wrt reference coordinates: + ! lamda = sqrt(C) ... C: right Cauchy-Green strain + ! lamda(stretch) wrt deformed coordinates: + ! lambda = sqrt(1/(1-2e)) + ! with: e(x) = 1/2 * (g - c(x)) + ! c(x) = F^-T * F^-1 + ! hence we can express lambda wrt to deformed coordinates + ! calculate Cauchy deformation tensor, using the expression + ! F^-T * F^-1 = (F * F^T)^-1 + CALL IdentityMatrix(AZL_INV,err,error,*999) ! initialise left Cauchy-Green + CALL IdentityMatrix(CDT,err,error,*999) ! initialise right Cauchy-Green + CALL MatrixProduct(DZDNU,DZDNUT,AZL_INV,err,error,*999) + CALL Invert(AZL_INV,CDT,DET_CDT,err,error,*999) + CALL IdentityMatrix(CDT_TEMP,err,error,*999) + ! rotate Cauchy deformation tensor to be with respect to deformed material coordinates + ! c(nuo) = q^-1 * c(x) * q^-T + CALL MatrixProduct(CDT(1:numberOfXDimensions,1:numberOfXDimensions), & + & DNUODZT(1:numberOfXDimensions,1:numberOfXDimensions), & + & CDT_TEMP(1:numberOfXDimensions,1:numberOfXDimensions),err,error,*999) + CALL MatrixProduct(DNUODZ(1:numberOfXDimensions,1:numberOfXDimensions), & + & CDT_TEMP(1:numberOfXDimensions,1:numberOfXDimensions), & + & CDT_FBR(1:numberOfXDimensions,1:numberOfXDimensions),err,error,*999) + + IF(FIELD_VARIABLE%COMPONENTS(1)%INTERPOLATION_TYPE==FIELD_CONSTANT_INTERPOLATION) THEN + DO component_idx=1,numberOfXDimensions + DO j=1,numberOfXDimensions + lambda_g(component_idx,j)=1.0_DP + ENDDO + ENDDO + ELSE + DO component_idx=1,numberOfXDimensions ! initialise the matrix + DO j=1,numberOfXDimensions + lambda_g(component_idx,j)=1.0_DP + ENDDO + ENDDO + ! lambda = sqrt(1/c(nuo)) + DO component_idx=1,numberOfXDimensions + TEMP_LAMBDA=CDT(component_idx,component_idx) + IF(TEMP_LAMBDA .LT. 0.0_DP) THEN ! avoid negative square root + WRITE(*,*) '>>> Warning: Negative square root!' + IF(ABS(TEMP_LAMBDA) .LT. 1.e-6) THEN ! avoid division by zero + TEMP_LAMBDA=(10.0_DP)**(-6.0_DP) + ENDIF + TEMP_LAMBDA=ABS(TEMP_LAMBDA) + ENDIF + IF(TEMP_LAMBDA==0.0_DP) THEN ! avoid division by zero + WRITE(*,*) '>>> Warning: Division by zero!' + TEMP_LAMBDA=(10.0_DP)**(-6.0_DP) + ENDIF + lambda_g(component_idx,component_idx)=SQRT(1.0_DP/TEMP_LAMBDA) + ENDDO + ENDIF + + ! transform the Cauchy stress wrt orthogonal fibres to be wrt spatial coordinates + ! sigma(x) = q * sigma(nuo) * q^T + CALL IdentityMatrix(TEMP_ROT,err,error,*999) ! initialise temporary variable + CALL IdentityMatrix(CAUCHY_TENSOR_DEFGEO,err,error,*999) ! initialise Cauchy tensor wrt deformed spatial + CALL MatrixProduct(CAUCHY_TENSOR_DEFIBRE(1:numberOfXDimensions,1:numberOfXDimensions), & + & DZDNUOT(1:numberOfXDimensions,1:numberOfXDimensions), & + & TEMP_ROT(1:numberOfXDimensions,1:numberOfXDimensions),err,error,*999) + CALL MatrixProduct(DZDNUO(1:numberOfXDimensions,1:numberOfXDimensions), & + & TEMP_ROT(1:numberOfXDimensions,1:numberOfXDimensions), & + & CAUCHY_TENSOR_DEFGEO(1:numberOfXDimensions,1:numberOfXDimensions),err,error,*999) + ! calculate active stress components + DO component_idx=1,numberOfXDimensions + DO j=1,numberOfXDimensions + CAUCHY_TENSOR_DEFGEO(component_idx,j)=CAUCHY_TENSOR_DEFGEO(component_idx,j)* & + & (1.0_DP+1.45_DP*(lambda_g(component_idx,j)-1.0_DP)) + ENDDO + ENDDO + + ! add active stress components + DO component_idx=1,FIELD_VARIABLE%NUMBER_OF_COMPONENTS + DO j=1, FIELD_VARIABLE%NUMBER_OF_COMPONENTS + CAUCHY_TENSOR(component_idx,j)=CAUCHY_TENSOR(component_idx,j)+ & + & CAUCHY_TENSOR_DEFGEO(component_idx,j) + ENDDO + ENDDO + ENDIF + CAUCHY_TENSOR=CAUCHY_TENSOR/Jznu + + IF(EQUATIONS_SET_SUBTYPE==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE) THEN + ! rotate the cauchy stress wrt spatial coordinates to be wrt deformed orthogonal fibres + ! for exporting and evaluating the stress in fibre direction + ! sigma(x) = q * sigma(nuo) * q^T + ! sigma(nuo) = q^-1 * sigma(x) * q^-T + CALL MatrixProduct(CAUCHY_TENSOR(1:numberOfXDimensions,1:numberOfXDimensions), & + & DZDNUOT(1:numberOfXDimensions,1:numberOfXDimensions), & + & TEMP_ROT2(1:numberOfXDimensions,1:numberOfXDimensions),err,error,*999) + CALL MatrixProduct(DZDNUO(1:numberOfXDimensions,1:numberOfXDimensions), & + & TEMP_ROT2(1:numberOfXDimensions,1:numberOfXDimensions), & + & CAUCHY_TENSOR_FBR(1:numberOfXDimensions,1:numberOfXDimensions),err,error,*999) + ! get the deformation gradient which maps reference and deformed orthogonal fibre coordinates + ! dnuo/dnu = dz/dnu * dnuo/dz + CALL MatrixProduct(DZDNU(1:numberOfXDimensions,1:numberOfXDimensions), & + & DNUODZ(1:numberOfXDimensions,1:numberOfXDimensions), & + & DNUODNU(1:numberOfXDimensions,1:numberOfXDimensions),err,error,*999) + ENDIF IF(DIAGNOSTICS1) THEN CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," ELEMENT_NUMBER = ",ELEMENT_NUMBER,err,error,*999) CALL WRITE_STRING_VALUE(DIAGNOSTIC_OUTPUT_TYPE," gauss_idx = ",GAUSS_POINT_NUMBER,err,error,*999) @@ -7250,7 +7582,8 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO STRESS_TENSOR(1:3)=STRESS_TENSOR(1:3)-ONETHIRD_TRACE+P CASE(EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE,EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & - & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE,EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE) IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE) THEN PRESSURE_COMPONENT=GEOMETRIC_INTERPOLATED_POINT%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS P=GEOMETRIC_INTERPOLATED_POINT%VALUES(PRESSURE_COMPONENT,NO_PART_DERIV) @@ -7264,7 +7597,9 @@ SUBROUTINE FINITE_ELASTICITY_GAUSS_STRESS_TENSOR(EQUATIONS_SET,DEPENDENT_INTERPO TEMPTERM1=0.5_DP*C(1)*EXP(0.5_DP*DOT_PRODUCT(E,DQ_DE)) ! Calculate isochoric fictitious 2nd Piola tensor (in Voigt form) STRESS_TENSOR=TEMPTERM1*DQ_DE - IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE) THEN + IF(EQUATIONS_SET%specification(3)==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE .OR. & + & EQUATIONS_SET%specification(3)==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE .OR. & + & EQUATIONS_SET%specification(3)==EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE) THEN !add active contraction stress values !Be aware for modified DZDNU, should active contraction be added here? Normally should be okay as modified DZDNU and DZDNU !converge during the Newton iteration. @@ -7586,8 +7921,10 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET & EQUATIONS_SET_ORTHOTROPIC_MATERIAL_HOLZAPFEL_OGDEN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & - & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_LAW_IN_CELLML_EVALUATE_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, & & EQUATIONS_SET_GROWTH_LAW_IN_CELLML_SUBTYPE, & @@ -7697,6 +8034,8 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_LAW_IN_CELLML_EVALUATE_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, & & EQUATIONS_SET_GROWTH_LAW_IN_CELLML_SUBTYPE, & @@ -7782,8 +8121,11 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET & EQUATIONS_SET_COMPRESSIBLE_ACTIVECONTRACTION_SUBTYPE, & & EQUATIONS_SET_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_NO_SUBTYPE,EQUATIONS_SET_MEMBRANE_SUBTYPE, & & EQUATIONS_SET_ORTHOTROPIC_MATERIAL_HOLZAPFEL_OGDEN_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & - & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_HUMPHREY_YIN_SUBTYPE, & & EQUATIONS_SET_STANDARD_MONODOMAIN_ELASTICITY_SUBTYPE,EQUATIONS_SET_1D3D_MONODOMAIN_ELASTICITY_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & @@ -9703,7 +10045,8 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET !Mooney Rivlin, St Venant Kirchoff and Compressible active contraction subtype CASE(EQUATIONS_SET_MOONEY_RIVLIN_ACTIVECONTRACTION_SUBTYPE,EQUATIONS_SET_STVENANT_KIRCHOFF_ACTIVECONTRACTION_SUBTYPE, & & EQUATIONS_SET_COMPRESSIBLE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_HOLZAPFEL_OGDEN_ACTIVECONTRACTION_SUBTYPE, & - & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE) + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE) NUMBER_OF_COMPONENTS = 3 !one contractile stress value for each of the three directions IF(EQUATIONS_SET%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN CALL FIELD_CREATE_START(EQUATIONS_SET_SETUP%FIELD_USER_NUMBER,EQUATIONS_SET%REGION,EQUATIONS_SET%INDEPENDENT% & @@ -10419,6 +10762,8 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET NUMBER_OF_COMPONENTS=3; CASE(EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_TRANSVERSE_ISOTROPIC_HUMPHREY_YIN_SUBTYPE) NUMBER_OF_COMPONENTS=4; @@ -10863,6 +11208,18 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET CALL FIELD_VARIABLE_LABEL_SET(EQUATIONS_SET%derived%derivedField,variableType,"Strain",err,error,*999) CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%derived%derivedField,variableType, & & numberOfFullTensorComponents,err,error,*999) + CASE(EQUATIONS_SET_DEFORMATION_GRADIENT_TENSOR_SPATIAL) + CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%derived%derivedField,variableType, & + & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) + CALL FIELD_VARIABLE_LABEL_SET(EQUATIONS_SET%derived%derivedField,variableType,"Strain",err,error,*999) + CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%derived%derivedField,variableType, & + & numberOfFullTensorComponents,err,error,*999) + CASE(EQUATIONS_SET_DEFORMATION_GRADIENT_TENSOR_FIBRE) + CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%derived%derivedField,variableType, & + & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) + CALL FIELD_VARIABLE_LABEL_SET(EQUATIONS_SET%derived%derivedField,variableType,"Strain",err,error,*999) + CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%derived%derivedField,variableType, & + & numberOfFullTensorComponents,err,error,*999) CASE(EQUATIONS_SET_R_CAUCHY_GREEN_DEFORMATION_TENSOR) CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%derived%derivedField,variableType, & & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) @@ -10887,6 +11244,12 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET CALL FIELD_VARIABLE_LABEL_SET(EQUATIONS_SET%derived%derivedField,variableType,"Stress",err,error,*999) CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%derived%derivedField,variableType, & & numberOfSymmetricTensorComponents,err,error,*999) + CASE(EQUATIONS_SET_CAUCHY_STRESS_TENSOR_FIBRE) + CALL FIELD_DIMENSION_SET_AND_LOCK(EQUATIONS_SET%derived%derivedField,variableType, & + & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) + CALL FIELD_VARIABLE_LABEL_SET(EQUATIONS_SET%derived%derivedField,variableType,"Stress",err,error,*999) + CALL FIELD_NUMBER_OF_COMPONENTS_SET_AND_LOCK(EQUATIONS_SET%derived%derivedField,variableType, & + & numberOfSymmetricTensorComponents,err,error,*999) CASE(EQUATIONS_SET_FIRST_PK_STRESS_TENSOR) CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_SECOND_PK_STRESS_TENSOR) @@ -11033,6 +11396,8 @@ SUBROUTINE FiniteElasticity_EquationsSetSolutionMethodSet(EQUATIONS_SET,SOLUTION & EQUATIONS_SET_MULTISCALE_ACTIVE_STRAIN_SUBTYPE, & & EQUATIONS_SET_1D3D_MONODOMAIN_ACTIVE_STRAIN_SUBTYPE, & & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & & EQUATIONS_SET_RATE_BASED_SMOOTH_MODEL_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_SMOOTH_MODEL_SUBTYPE, & & EQUATIONS_SET_RATE_BASED_GROWTH_MODEL_SUBTYPE,EQUATIONS_SET_COMPRESSIBLE_RATE_BASED_GROWTH_MODEL_SUBTYPE) SELECT CASE(SOLUTION_METHOD) @@ -11109,8 +11474,11 @@ SUBROUTINE FiniteElasticity_EquationsSetSpecificationSet(equationsSet,specificat & EQUATIONS_SET_ORTHOTROPIC_MATERIAL_HOLZAPFEL_OGDEN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE,EQUATIONS_SET_NEARLY_INCOMPRESSIBLE_MOONEY_RIVLIN_SUBTYPE, & & EQUATIONS_SET_INCOMPRESSIBLE_ELAST_MULTI_COMP_DARCY_SUBTYPE,EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & - & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & & EQUATIONS_SET_REFERENCE_STATE_MOONEY_RIVLIN_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & + & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_LAW_IN_CELLML_EVALUATE_SUBTYPE, & & EQUATIONS_SET_CONSTITUTIVE_AND_GROWTH_LAW_IN_CELLML_SUBTYPE, & & EQUATIONS_SET_GROWTH_LAW_IN_CELLML_SUBTYPE, & diff --git a/src/opencmiss_iron.F90 b/src/opencmiss_iron.F90 index deae8cb4..b5c44b93 100644 --- a/src/opencmiss_iron.F90 +++ b/src/opencmiss_iron.F90 @@ -2591,6 +2591,10 @@ MODULE OpenCMISS_Iron & EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE !< Transverse isotropic Guccione constitutive law for finite elasticity equstions set subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE = & & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE !< Transverse isotropic Guccione constitutive law with active contraction subtype \see OpenCMISS_EquationsSetSubtypes,OpenCMISS + INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE = & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE !< Transverse isotropic Guccione constitutive law with active contraction subtype without length dependence \see OpenCMISS_EquationsSetSubtypes,OpenCMISS + INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE = & + & EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE !< Transverse isotropic Guccione constitutive law with active contraction subtype for benchmark problem \see OpenCMISS_EquationsSetSubtypes,OpenCMISS INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_INCOMPRESS_FINITE_ELASTICITY_DARCY_SUBTYPE= & & EQUATIONS_SET_INCOMPRESSIBLE_FINITE_ELASTICITY_DARCY_SUBTYPE ! \see OpenCMISS::Iron::EquationsSet,OpenCMISS !>@{ INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_DERIVED_DEFORMATION_GRADIENT = EQUATIONS_SET_DEFORMATION_GRADIENT_TENSOR !@} @@ -3136,6 +3145,8 @@ MODULE OpenCMISS_Iron & CMFE_EQUATIONS_SET_INCOMPRESSIBLE_ELASTICITY_DRIVEN_MR_SUBTYPE, & & CMFE_EQUATIONS_SET_INCOMPRESS_ELAST_MULTI_COMP_DARCY_SUBTYPE,CMFE_EQUATIONS_SET_TRANSVERSE_ISOTROPIC_GUCCIONE_SUBTYPE, & & CMFE_EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_SUBTYPE, & + & CMFE_EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_NOLENDEP_SUBTYPE, & + & CMFE_EQUATIONS_SET_GUCCIONE_ACTIVECONTRACTION_2NDPIOLA_SUBTYPE, & & CMFE_EQUATIONS_SET_ACTIVE_STRAIN_SUBTYPE, & & CMFE_EQUATIONS_SET_MULTISCALE_ACTIVE_STRAIN_SUBTYPE, & & CMFE_EQUATIONS_SET_REFERENCE_STATE_TRANSVERSE_GUCCIONE_SUBTYPE, & @@ -3236,10 +3247,11 @@ MODULE OpenCMISS_Iron & CMFE_EQUATIONS_SET_FV_SOLUTION_METHOD,CMFE_EQUATIONS_SET_GFEM_SOLUTION_METHOD,CMFE_EQUATIONS_SET_GFD_SOLUTION_METHOD, & & CMFE_EQUATIONS_SET_GFV_SOLUTION_METHOD - PUBLIC CMFE_EQUATIONS_SET_DERIVED_DEFORMATION_GRADIENT,CMFE_EQUATIONS_SET_DERIVED_R_CAUCHY_GREEN_DEFORMATION, & + PUBLIC CMFE_EQUATIONS_SET_DERIVED_DEFORMATION_GRADIENT,CMFE_EQUATIONS_SET_DERIVED_DEFORMATION_GRADIENT_SPATIAL, & + & CMFE_EQUATIONS_SET_DERIVED_DEFORMATION_GRADIENT_FIBRE,CMFE_EQUATIONS_SET_DERIVED_R_CAUCHY_GREEN_DEFORMATION, & & CMFE_EQUATIONS_SET_DERIVED_L_CAUCHY_GREEN_DEFORMATION,CMFE_EQUATIONS_SET_DERIVED_GREEN_LAGRANGE_STRAIN, & - & CMFE_EQUATIONS_SET_DERIVED_CAUCHY_STRESS,CMFE_EQUATIONS_SET_DERIVED_FIRST_PK_STRESS, & - & CMFE_EQUATIONS_SET_DERIVED_SECOND_PK_STRESS + & CMFE_EQUATIONS_SET_DERIVED_CAUCHY_STRESS,CMFE_EQUATIONS_SET_DERIVED_CAUCHY_STRESS_FIBRE, & + & CMFE_EQUATIONS_SET_DERIVED_FIRST_PK_STRESS, CMFE_EQUATIONS_SET_DERIVED_SECOND_PK_STRESS PUBLIC CMFE_EQUATIONS_MATRIX_STIFFNESS,CMFE_EQUATIONS_MATRIX_DAMPING,CMFE_EQUATIONS_MATRIX_MASS From 4f01130b6650eceb21fd441ff88c7b498949a532 Mon Sep 17 00:00:00 2001 From: Prasad Date: Tue, 26 Nov 2019 16:40:15 +1300 Subject: [PATCH 17/17] Updates for projections, fitting, and finite elasticity finite difference jacobian --- src/data_projection_routines.F90 | 49 +++++++++++++++++++++++--------- src/equations_set_routines.F90 | 4 +-- src/fitting_routines.F90 | 8 +++--- src/opencmiss_iron.F90 | 15 +++++----- 4 files changed, 49 insertions(+), 27 deletions(-) diff --git a/src/data_projection_routines.F90 b/src/data_projection_routines.F90 index 1c0e3d2e..74e72d63 100644 --- a/src/data_projection_routines.F90 +++ b/src/data_projection_routines.F90 @@ -4508,20 +4508,11 @@ SUBROUTINE DataProjection_ProjectionTypeSet(dataProjection,projectionType,err,er IF(dataProjection%numberOfXi/=SIZE(dataProjection%startingXi,1)) THEN ALLOCATE(startingXi(dataProjection%datapoints%numberOfDatapoints,dataProjection%numberOfXi),STAT=err) IF(err/=0) CALL FlagError("Could not allocate starting xi.",err,error,*999) - IF(dataProjection%numberOfXi>SIZE(dataProjection%startingXi,1)) THEN - startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:SIZE(dataProjection%startingXi,1))= & - & dataProjection%startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:SIZE(dataProjection%startingXi,1)) - startingXi(1:dataProjection%datapoints%numberOfDatapoints,SIZE(dataProjection%startingXi,1):dataProjection%numberOfXi)= & - & 0.5_DP - ELSE - startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:SIZE(dataProjection%startingXi,1))= & - & dataProjection%startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:dataProjection%numberOfXi) - ENDIF + startingXi=0.5_DP IF(ALLOCATED(dataProjection%startingXi)) DEALLOCATE(dataProjection%startingXi) ALLOCATE(dataProjection%startingXi(1:dataProjection%datapoints%numberOfDatapoints,dataProjection%numberOfXi),STAT=err) IF(err/=0) CALL FlagError("Could not allocate data projection starting xi.",err,error,*999) - dataProjection%startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:dataProjection%numberOfXi)= & - & startingXi(1:dataProjection%datapoints%numberOfDatapoints,1:dataProjection%numberOfXi) + dataProjection%startingXi=startingXi DEALLOCATE(startingXi) ENDIF @@ -5127,7 +5118,7 @@ SUBROUTINE DataProjection_ResultElementNumberSet(dataProjection,dataPointUserNum IF(.NOT.ASSOCIATED(dataProjection)) CALL FlagError("Data projection is not associated.",err,error,*999) IF(.NOT.dataProjection%dataProjectionFinished) CALL FlagError("Data projection has not been finished.",err,error,*999) - IF(dataProjection%dataProjectionProjected) CALL FlagError("Data projection has already been projected.",err,error,*999) + !IF(dataProjection%dataProjectionProjected) CALL FlagError("Data projection has already been projected.",err,error,*999) SELECT CASE(dataProjection%projectionType) CASE(DATA_PROJECTION_BOUNDARY_FACES_PROJECTION_TYPE) @@ -5147,6 +5138,7 @@ SUBROUTINE DataProjection_ResultElementNumberSet(dataProjection,dataPointUserNum IF(elementExists) THEN IF(.NOT.ghostElement) THEN dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementLocalNumber=elementLocalNumber + dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementGlobalNumber=elementLocalNumber dataProjection%dataProjectionResults(dataPointGlobalNumber)%exitTag=DATA_PROJECTION_USER_SPECIFIED ENDIF ELSE @@ -5532,14 +5524,20 @@ SUBROUTINE DataProjection_ResultProjectionXiSet(dataProjection,dataPointUserNumb INTEGER(INTG), INTENT(OUT) :: err !dataProjection%decomposition + NULLIFY(decompositionTopology) + CALL Decomposition_TopologyGet(decomposition,decompositionTopology,err,error,*999) + NULLIFY(domain) + CALL Decomposition_DomainGet(decomposition,0,domain,err,error,*999) + NULLIFY(domainTopology) + CALL Domain_TopologyGet(domain,domainTopology,err,error,*999) + NULLIFY(domainElements) + CALL DomainTopology_ElementsGet(domainTopology,domainElements,err,error,*999) + + IF(dataProjection%numberOfXi==dataProjection%numberOfElementXi) THEN + dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementXi= & + & dataProjection%dataProjectionResults(dataPointGlobalNumber)%xi + ELSE + elementNumber=dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementLocalNumber + localLineFaceNumber=dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementLineFaceNumber + basis=>domainElements%elements(elementNumber)%basis + CALL Basis_BoundaryXiToXi(basis,localLineFaceNumber,dataProjection% & + & dataProjectionResults(dataPointGlobalNumber)%xi(1:dataProjection%numberOfXi),dataProjection% & + & dataProjectionResults(dataPointGlobalNumber)%elementXi,err,error,*999) + ENDIF + EXITS("DataProjection_ResultProjectionXiSet") RETURN 999 ERRORSEXITS("DataProjection_ResultProjectionXiSet",err,error) diff --git a/src/equations_set_routines.F90 b/src/equations_set_routines.F90 index 9fd95202..442b7882 100644 --- a/src/equations_set_routines.F90 +++ b/src/equations_set_routines.F90 @@ -2911,7 +2911,7 @@ SUBROUTINE EquationsSet_FiniteElementJacobianEvaluateFD(equationsSet,elementNumb DO componentIdx=1,columnVariable%NUMBER_OF_COMPONENTS ! adjusts the pertubation for the hydrostatic pressure to not use the L2Norm IF(componentIdx==4) THEN - delta=(1.0_DP)*equationsSet%equations%jacobianFiniteDifferenceStepSize + delta=(1.0_DP)*nonlinearMatrices%jacobians(jacobianNumber)%ptr%jacobianFiniteDifferenceStepSize ENDIF elementsTopology=>columnVariable%COMPONENTS(componentIdx)%DOMAIN%TOPOLOGY%ELEMENTS componentInterpolationType=columnVariable%COMPONENTS(componentIdx)%INTERPOLATION_TYPE @@ -3019,7 +3019,7 @@ SUBROUTINE EquationsSet_FiniteElementResidualEvaluate(equationsSet,elementNumber CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) CALL ClassicalField_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*999) CASE(EQUATIONS_SET_FITTING_CLASS) - CALL Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*999) + CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_MODAL_CLASS) diff --git a/src/fitting_routines.F90 b/src/fitting_routines.F90 index 1ac6819d..b73e13f6 100644 --- a/src/fitting_routines.F90 +++ b/src/fitting_routines.F90 @@ -3016,8 +3016,8 @@ SUBROUTINE Fitting_EquationsSetDataSetup(equationsSet,equationsSetSetup,err,erro CALL Field_NumberOfVariablesCheck(equationsSetSetup%field,2,err,error,*999) CALL Field_VariableTypesCheck(equationsSetSetup%field,[FIELD_U_VARIABLE_TYPE,FIELD_DELUDELN_VARIABLE_TYPE], & & err,error,*999) - CALL Field_DimensionCheck(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,FIELD_SCALAR_DIMENSION_TYPE,err,error,*999) - CALL Field_DimensionCheck(equationsSetSetup%field,FIELD_DELUDELN_VARIABLE_TYPE,FIELD_SCALAR_DIMENSION_TYPE, & + CALL Field_DimensionCheck(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) + CALL Field_DimensionCheck(equationsSetSetup%field,FIELD_DELUDELN_VARIABLE_TYPE,FIELD_VECTOR_DIMENSION_TYPE, & & err,error,*999) CALL Field_DataTypeCheck(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,FIELD_DP_TYPE,err,error,*999) CALL Field_DataTypeCheck(equationsSetSetup%field,FIELD_DELUDELN_VARIABLE_TYPE,FIELD_DP_TYPE,err,error,*999) @@ -3357,9 +3357,9 @@ SUBROUTINE Fitting_EquationsSetDataSetup(equationsSet,equationsSetSetup,err,erro CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD) DO componentIdx=1,numberOfComponents CALL Field_ComponentInterpolationCheck(equationsSetSetup%FIELD,FIELD_U_VARIABLE_TYPE,componentIdx, & - & FIELD_GAUSS_POINT_BASED_INTERPOLATION,err,error,*999) + & FIELD_DATA_POINT_BASED_INTERPOLATION,err,error,*999) CALL Field_ComponentInterpolationCheck(equationsSetSetup%FIELD,FIELD_V_VARIABLE_TYPE,componentIdx, & - & FIELD_GAUSS_POINT_BASED_INTERPOLATION,err,error,*999) + & FIELD_DATA_POINT_BASED_INTERPOLATION,err,error,*999) ENDDO !componentIdx CASE DEFAULT localError="The solution method of "//TRIM(NumberToVString(equationsSet%SOLUTION_METHOD, & diff --git a/src/opencmiss_iron.F90 b/src/opencmiss_iron.F90 index b5c44b93..48cf3fa0 100644 --- a/src/opencmiss_iron.F90 +++ b/src/opencmiss_iron.F90 @@ -21377,8 +21377,8 @@ SUBROUTINE cmfe_DataProjection_ProjectionDataCandidateElementsSetObj00(dataProje ENTERS("cmfe_DataProjection_ProjectionDataCandidateElementsSetObj00",err,error,*999) - CALL cmfe_DataProjection_ProjectionDataCandidateElementsSetObj11(dataProjection,[dataPointUserNumber], & - & [candidateElementUserNumber],err) + CALL DataProjection_ProjectionDataCandidateElementsSet(dataProjection%dataProjection,[dataPointUserNumber], & + & [candidateElementUserNumber],err,error,*999) EXITS("cmfe_DataProjection_ProjectionDataCandidateElementsSetObj00") RETURN @@ -21407,8 +21407,8 @@ SUBROUTINE cmfe_DataProjection_ProjectionDataCandidateElementsSetObj01(dataProje ENTERS("cmfe_DataProjection_ProjectionDataCandidateElementsSetObj01",err,error,*999) - CALL cmfe_DataProjection_ProjectionDataCandidateElementsSetObj11(dataProjection,[dataPointUserNumber], & - & candidateElementUserNumbers,err) + CALL DataProjection_ProjectionDataCandidateElementsSet(dataProjection%dataProjection,[dataPointUserNumber], & + & candidateElementUserNumbers,err,error,*999) EXITS("cmfe_DataProjection_ProjectionDataCandidateElementsSetObj01") RETURN @@ -21437,8 +21437,8 @@ SUBROUTINE cmfe_DataProjection_ProjectionDataCandidateElementsSetObj10(dataProje ENTERS("cmfe_DataProjection_ProjectionDataCandidateElementsSetObj10",err,error,*999) - CALL cmfe_DataProjection_ProjectionDataCandidateElementsSetObj11(dataProjection,dataPointUserNumbers, & - & [candidateElementUserNumber],err) + CALL DataProjection_ProjectionDataCandidateElementsSet(dataProjection%dataProjection,dataPointUserNumbers, & + & [candidateElementUserNumber],err,error,*999) EXITS("cmfe_DataProjection_ProjectionDataCandidateElementsSetObj10") RETURN @@ -21467,7 +21467,8 @@ SUBROUTINE cmfe_DataProjection_ProjectionDataCandidateElementsSetObj11(dataProje ENTERS("cmfe_DataProjection_ProjectionDataCandidateElementsSetObj11",err,error,*999) - CALL DataProjection_ProjectionCandidateElementsSet(dataProjection%dataProjection,candidateElementUserNumbers,err,error,*999) + CALL DataProjection_ProjectionDataCandidateElementsSet(dataProjection%dataProjection,dataPointUserNumbers, & + & candidateElementUserNumbers,err,error,*999) EXITS("cmfe_DataProjection_ProjectionDataCandidateElementsSetObj11") RETURN