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..74e72d63 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 @@ -142,6 +142,8 @@ MODULE DataProjectionRoutines PUBLIC DataProjection_Destroy + PUBLIC DataProjection_DataPointFieldEvaluate + PUBLIC DataProjection_DataPointsProjectionEvaluate PUBLIC DataProjection_DataPointsPositionEvaluate @@ -182,13 +184,19 @@ 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_ResultElementXiGet + PUBLIC DataProjection_ResultExitTagGet PUBLIC DataProjection_ResultProjectionVectorGet - PUBLIC DataProjection_ResultXiGet,DataProjection_ResultXiSet + PUBLIC DataProjection_ResultProjectionXiGet,DataProjection_ResultProjectionXiSet PUBLIC DataProjection_StartingXiGet,DataProjection_StartingXiSet @@ -262,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 !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 +946,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 +979,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 @@ -1034,7 +1076,7 @@ SUBROUTINE DataProjection_DataProjectionResultsInitialise(dataProjection,err,err CALL FlagError(localError,err,error,*999) ENDIF dataProjection%dataProjectionResults(dataPointIdx)%xi(1:dataProjection%numberOfXi)= & - & dataProjection%startingXi(1:dataProjection%numberOfXi) + & dataProjection%startingXi(dataPointIdx,1:dataProjection%numberOfXi) dataProjection%dataProjectionResults(dataPointIdx)%elementXi=0.0_DP dataProjection%dataProjectionResults(dataPointIdx)%projectionVector(1:dataProjection%numberOfCoordinates)=0.0_DP ENDDO !dataPointIdx @@ -1290,6 +1332,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)%elementXi, & + & 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 + ! !================================================================================================================================ ! @@ -1372,7 +1487,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, & @@ -1582,6 +1697,7 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection !find the clostest elements/faces/lines for each point in the current computational node base on starting xi !the clostest elements/faces/lines are required to shrink down on the list of possible projection candiates numberOfClosestCandidates=MIN(dataProjection%numberOfClosestElements,dataProjection%maxNumberOfCandidates) + numberOfClosestCandidates=dataProjection%maxNumberOfCandidates !Allocated and store he information for each data point. The information has to be stored in the corresponding !rows for them to be contiguous in memory for easy MPI access ALLOCATE(closestElements(numberOfDataPoints,numberOfClosestCandidates),STAT=err) @@ -1608,12 +1724,14 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection & dataPoints%dataPoints(dataPointIdx)%position,dataNumberOfCandidates,dataProjection% & & dataProjectionCandidates(dataPointIdx)%candidateElementNumbers, & & dataProjection%dataProjectionCandidates(dataPointIdx)%localFaceLineNumbers, & + & dataProjection%startingXi(dataPointIdx,:), & & closestElements(dataPointIdx,:),closestLinesFaces(dataPointIdx,:), & & closestDistances(dataPointIdx,:),err,error,*999) ELSE CALL DataProjection_ClosestLinesFind(dataProjection,interpolatedPoint, & & dataPoints%dataPoints(dataPointIdx)%position,numberOfCandidates,candidateElements, & - & candidateLinesFaces,closestElements(dataPointIdx,:),closestLinesFaces(dataPointIdx,:), & + & candidateLinesFaces,dataProjection%startingXi(dataPointIdx,:), & + & closestElements(dataPointIdx,:),closestLinesFaces(dataPointIdx,:), & & closestDistances(dataPointIdx,:),err,error,*999) ENDIF ENDDO !dataPointIdx @@ -1627,13 +1745,14 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection & dataPoints%dataPoints(dataPointIdx)%position,dataNumberOfCandidates,dataProjection% & & dataProjectionCandidates(dataPointIdx)%candidateElementNumbers, & & dataProjection%dataProjectionCandidates(dataPointIdx)%localFaceLineNumbers, & + & dataProjection%startingXi(dataPointIdx,:), & & closestElements(dataPointIdx,:),closestLinesFaces(dataPointIdx,:), & & closestDistances(dataPointIdx,:),err,error,*999) ELSE CALL DataProjection_ClosestFacesFind(dataProjection,interpolatedPoint, & & dataPoints%dataPoints(dataPointIdx)%position,numberOfCandidates,candidateElements, & - & candidateLinesFaces,closestElements(dataPointIdx,:),closestLinesFaces(dataPointIdx,:), & - & closestDistances(dataPointIdx,:),err,error,*999) + & candidateLinesFaces,dataProjection%startingXi(dataPointIdx,:),closestElements(dataPointIdx,:), & + & closestLinesFaces(dataPointIdx,:),closestDistances(dataPointIdx,:),err,error,*999) ENDIF ENDDO !dataPointIdx CASE(DATA_PROJECTION_ALL_ELEMENTS_PROJECTION_TYPE) @@ -1644,11 +1763,13 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_ClosestElementsFind(dataProjection,interpolatedPoint, & & dataPoints%dataPoints(dataPointIdx)%position,dataNumberOfCandidates, & & dataProjection%dataProjectionCandidates(dataPointIdx)%candidateElementNumbers, & - & closestElements(dataPointIdx,:),closestDistances(dataPointIdx,:),err,error,*999) + & closestElements(dataPointIdx,:),dataProjection%startingXi(dataPointIdx,:),closestDistances(dataPointIdx,:),& + & err,error,*999) ELSE CALL DataProjection_ClosestElementsFind(dataProjection,interpolatedPoint, & & dataPoints%dataPoints(dataPointIdx)%position,numberOfCandidates,candidateElements, & - & closestElements(dataPointIdx,:),closestDistances(dataPointIdx,:),err,error,*999) + & closestElements(dataPointIdx,:),dataProjection%startingXi(dataPointIdx,:),closestDistances(dataPointIdx,:),& + & err,error,*999) ENDIF ENDDO !dataPointIdx CASE DEFAULT @@ -1673,8 +1794,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 +1847,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 @@ -1743,30 +1868,32 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection !Newton project to closest lines, 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_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) !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, & - & 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) @@ -1775,41 +1902,44 @@ 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), & - & 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 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), & - & 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 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), & - & 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 +1971,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 +2005,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 +2017,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,21 +2037,28 @@ 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 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% & - & 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 +2068,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 +2118,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 +2142,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 +2164,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 +2225,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, & @@ -2230,7 +2383,7 @@ SUBROUTINE DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPo INTEGER(INTG), INTENT(OUT) :: projectionExitTag !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=projectionXi + 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=projectionXi + 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,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), & + & 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=projectionXi + 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 - 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) + 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=projectionXi + 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))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 + 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 - 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 - 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 + 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 - !OK model: keep the current region size + !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.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 + 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=projectionXi + 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)%distanceSIZE(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 - ELSE - startingXi(1:SIZE(dataProjection%startingXi,1))=dataProjection%startingXi(1:dataProjection%numberOfXi) - ENDIF + startingXi=0.5_DP 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=startingXi DEALLOCATE(startingXi) ENDIF @@ -4418,7 +4590,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 @@ -4452,11 +4625,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 @@ -4521,7 +4702,7 @@ SUBROUTINE DataProjection_ElementSet(dataProjection,dataPointUserNumber,elementU & " cannot be set as a projection element as it is a ghost element." CALL FlagError(localError,err,error,*999) ELSE - dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementNumber=elementUserNumber + dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementLocalNumber=localElementNumber ENDIF ELSE localError="The specificed user element number of "//TRIM(NumberToVString(elementUserNumber,"*",err,error))// & @@ -4697,7 +4878,8 @@ SUBROUTINE DataProjection_ResultAnalysisOutput(dataProjection,filename,err,error TYPE(VARYING_STRING), INTENT(OUT) :: error !=1) THEN @@ -4889,7 +5079,8 @@ SUBROUTINE DataProjection_ResultElementNumberGet(dataProjection,dataPointUserNum IF(.NOT.dataProjection%dataProjectionProjected) CALL FlagError("Data projection has not been projected.",err,error,*999) CALL DataProjection_DataPointGlobalNumberGet(dataProjection,dataPointUserNumber,dataPointGlobalNumber,err,error,*999) - projectionElementNumber=dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementNumber +!!TODO: this should return the element user number + projectionElementNumber=dataProjection%dataProjectionResults(dataPointGlobalNumber)%elementGlobalNumber EXITS("DataProjection_ResultElementNumberGet") RETURN @@ -4902,6 +5093,76 @@ END SUBROUTINE DataProjection_ResultElementNumberGet !================================================================================================================================ ! + !>Sets 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 +5201,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 +5319,130 @@ 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 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,*) @@ -5011,7 +5476,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 !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_ResultXiSet",err,error) +999 ERRORSEXITS("DataProjection_ResultProjectionXiSet",err,error) RETURN 1 - END SUBROUTINE DataProjection_ResultXiSet + END SUBROUTINE DataProjection_ResultProjectionXiSet ! !================================================================================================================================ 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..442b7882 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)*nonlinearMatrices%jacobians(jacobianNumber)%ptr%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 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/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/finite_elasticity_routines.F90 b/src/finite_elasticity_routines.F90 index c5b8c450..d3e48a97 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) @@ -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. @@ -929,7 +907,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 +959,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 @@ -1431,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) @@ -1513,7 +1494,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 +1552,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 @@ -1590,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 @@ -1621,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 @@ -1631,7 +1618,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 @@ -1656,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 @@ -1667,7 +1657,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 @@ -1826,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, & @@ -1863,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 @@ -1897,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) @@ -2138,7 +2133,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 +2149,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 +2521,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 +2537,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 @@ -3221,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, & @@ -3332,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, & @@ -3554,7 +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, & @@ -3622,7 +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, & @@ -3700,9 +3707,9 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field TYPE(VARYING_STRING), INTENT(OUT) :: error !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 +5270,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% & @@ -5218,25 +5418,52 @@ END SUBROUTINE FiniteElasticity_GaussDeformationGradientTensor !>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) @@ -7022,10 +7536,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) @@ -7061,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) @@ -7075,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. @@ -7322,7 +7846,8 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET !Local Variables INTEGER(INTG) :: GEOMETRIC_MESH_COMPONENT,GEOMETRIC_SCALING_TYPE,NUMBER_OF_COMPONENTS, & & numberOfDimensions,NUMBER_OF_DARCY_COMPONENTS,GEOMETRIC_COMPONENT_NUMBER,NUMBER_OF_COMPONENTS_2,component_idx, & - & componentIdx,derivedIdx,varIdx,variableType,NUMBER_OF_FLUID_COMPONENTS,numberOfTensorComponents + & componentIdx,derivedIdx,varIdx,variableType,NUMBER_OF_FLUID_COMPONENTS,numberOfSymmetricTensorComponents, & + & numberOfFullTensorComponents TYPE(COORDINATE_SYSTEM_TYPE), POINTER :: coordinateSystem TYPE(DECOMPOSITION_TYPE), POINTER :: GEOMETRIC_DECOMPOSITION TYPE(FIELD_TYPE), POINTER :: ANALYTIC_FIELD,DEPENDENT_FIELD,GEOMETRIC_FIELD @@ -7396,7 +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, & @@ -7481,6 +8009,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,& @@ -7505,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, & @@ -7590,7 +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, & @@ -9510,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% & @@ -10167,6 +10703,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, & @@ -10225,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; @@ -10651,7 +11190,8 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET VARIABLE_TYPES(varIdx)=EQUATIONS_SET%derived%variableTypes(derivedIdx) END IF END DO - numberOfTensorComponents=NUMBER_OF_VOIGT(numberOfDimensions) + numberOfSymmetricTensorComponents=NUMBER_OF_VOIGT(numberOfDimensions) + numberOfFullTensorComponents=numberOfDimensions*numberOfDimensions IF(EQUATIONS_SET%derived%derivedFieldAutoCreated) THEN CALL FIELD_NUMBER_OF_VARIABLES_SET_AND_LOCK(EQUATIONS_SET%derived%derivedField, & & EQUATIONS_SET%derived%numberOfVariables,err,error,*999) @@ -10667,31 +11207,49 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET & 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, & - & numberOfTensorComponents,err,error,*999) + & 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) 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, & - & numberOfTensorComponents,err,error,*999) + & numberOfSymmetricTensorComponents,err,error,*999) CASE(EQUATIONS_SET_L_CAUCHY_GREEN_DEFORMATION_TENSOR) 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, & - & numberOfTensorComponents,err,error,*999) + & numberOfSymmetricTensorComponents,err,error,*999) CASE(EQUATIONS_SET_GREEN_LAGRANGE_STRAIN_TENSOR) 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, & - & numberOfTensorComponents,err,error,*999) + & numberOfSymmetricTensorComponents,err,error,*999) CASE(EQUATIONS_SET_CAUCHY_STRESS_TENSOR) 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, & - & numberOfTensorComponents,err,error,*999) + & 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) @@ -10718,27 +11276,27 @@ SUBROUTINE FINITE_ELASTICITY_EQUATIONS_SET_SETUP(EQUATIONS_SET,EQUATIONS_SET_SET CALL FIELD_DIMENSION_CHECK(EQUATIONS_SET%derived%derivedField,variableType, & & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET%derived%derivedField,variableType, & - & numberOfTensorComponents,err,error,*999) + & numberOfFullTensorComponents,err,error,*999) CASE(EQUATIONS_SET_R_CAUCHY_GREEN_DEFORMATION_TENSOR) CALL FIELD_DIMENSION_CHECK(EQUATIONS_SET%derived%derivedField,variableType, & & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET%derived%derivedField,variableType, & - & numberOfTensorComponents,err,error,*999) + & numberOfSymmetricTensorComponents,err,error,*999) CASE(EQUATIONS_SET_L_CAUCHY_GREEN_DEFORMATION_TENSOR) CALL FIELD_DIMENSION_CHECK(EQUATIONS_SET%derived%derivedField,variableType, & & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET%derived%derivedField,variableType, & - & numberOfTensorComponents,err,error,*999) + & numberOfSymmetricTensorComponents,err,error,*999) CASE(EQUATIONS_SET_GREEN_LAGRANGE_STRAIN_TENSOR) CALL FIELD_DIMENSION_CHECK(EQUATIONS_SET%derived%derivedField,variableType, & & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET%derived%derivedField,variableType, & - & numberOfTensorComponents,err,error,*999) + & numberOfSymmetricTensorComponents,err,error,*999) CASE(EQUATIONS_SET_CAUCHY_STRESS_TENSOR) CALL FIELD_DIMENSION_CHECK(EQUATIONS_SET%derived%derivedField,variableType, & & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) CALL FIELD_NUMBER_OF_COMPONENTS_CHECK(EQUATIONS_SET%derived%derivedField,variableType, & - & numberOfTensorComponents,err,error,*999) + & 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) @@ -10820,6 +11378,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, & @@ -10837,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) @@ -10913,7 +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/fitting_routines.F90 b/src/fitting_routines.F90 index b9650c12..b73e13f6 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,tauParam,kappaParam - REAL(DP) :: 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) @@ -257,12 +256,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% & @@ -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) @@ -339,8 +344,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,66 +382,98 @@ 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 + 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 @@ -558,89 +606,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 +894,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 +918,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 +997,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 +1042,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) @@ -1144,14 +1204,24 @@ 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) - 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) - CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) - CALL FlagError("Not implemented.",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 CASE(EQUATIONS_SET_FITTING_STRAIN_ENERGY_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) CASE DEFAULT @@ -1192,57 +1262,43 @@ 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 = 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 - CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) - CALL FlagError("Not implemented.",err,error,*999) + sum = (tension + curvature) * jacobianGaussWeight CASE(EQUATIONS_SET_FITTING_STRAIN_ENERGY_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) CASE DEFAULT @@ -2192,7 +2248,7 @@ SUBROUTINE Fitting_EquationsSetGaussSetup(equationsSet,equationsSetSetup,err,err TYPE(VARYING_STRING), INTENT(OUT) :: error !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..48cf3fa0 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 @@ -1840,6 +1840,7 @@ MODULE OpenCMISS_Iron !> \see OpenCMISS::Iron::DataProjection,OpenCMISS !>@{ INTEGER(INTG), PARAMETER :: CMFE_DATA_PROJECTION_CANCELLED = DATA_PROJECTION_CANCELLED !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 @@ -2112,25 +2120,49 @@ MODULE OpenCMISS_Iron MODULE PROCEDURE cmfe_DataProjection_ResultDistanceGetObj END INTERFACE cmfe_DataProjection_ResultDistanceGet - !>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 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 MODULE PROCEDURE cmfe_DataProjection_ResultExitTagGetNumber MODULE PROCEDURE cmfe_DataProjection_ResultExitTagGetObj @@ -2155,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 @@ -2175,8 +2207,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 @@ -2187,6 +2219,8 @@ MODULE OpenCMISS_Iron PUBLIC cmfe_DataProjection_Destroy + PUBLIC cmfe_DataProjection_DataPointFieldEvaluate + PUBLIC cmfe_DataProjection_DataPointsPositionEvaluate PUBLIC cmfe_DataProjection_ProjectionCancelByDataPoints @@ -2227,9 +2261,15 @@ 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_ResultElementXiGet PUBLIC cmfe_DataProjection_ResultExitTagGet @@ -2237,7 +2277,7 @@ MODULE OpenCMISS_Iron PUBLIC cmfe_DataProjection_ResultRMSErrorGet - PUBLIC cmfe_DataProjection_ResultXiGet,cmfe_DataProjection_ResultXiSet + PUBLIC cmfe_DataProjection_ResultProjectionXiGet,cmfe_DataProjection_ResultProjectionXiSet PUBLIC cmfe_DataProjection_ResultProjectionVectorGet @@ -2551,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 !@} @@ -3096,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, & @@ -3196,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 @@ -8554,6 +8606,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) @@ -20130,6 +20208,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) @@ -21167,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 @@ -21197,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 @@ -21227,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 @@ -21257,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 @@ -22794,6 +23005,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 +23145,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 +23289,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) @@ -23003,6 +23429,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) @@ -23206,9 +23700,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 !