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..303b5928 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 @@ -893,6 +901,35 @@ END SUBROUTINE DataProjection_DataProjectionCandidatesInitialise !================================================================================================================================ ! + !>Cancels a data projection result + SUBROUTINE DataProjection_DataProjectionResultCancel(dataProjectionResult,err,error,*) + + !Argument variables + TYPE(DataProjectionResultType) :: dataProjectionResult !Finalises the data projection result and deallocates all memory SUBROUTINE DataProjection_DataProjectionResultFinalise(dataProjectionResult,err,error,*) @@ -906,7 +943,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 +976,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 @@ -1290,6 +1329,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 +1484,7 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection INTEGER(INTG), ALLOCATABLE :: globalNumberOfProjectedPoints(:) INTEGER(INTG) :: MPIClosestDistances,dataProjectionGlobalNumber INTEGER(INTG) :: MPIIError - INTEGER(INTG), ALLOCATABLE :: projectedElement(:),projectedLineFace(:),projectionExitTag(:) + INTEGER(INTG), ALLOCATABLE :: projectedGlobalElement(:),projectedLineFace(:),projectedLocalElement(:),projectionExitTag(:) REAL(DP), ALLOCATABLE :: closestDistances(:,:),globalClosestDistances(:,:) REAL(DP), ALLOCATABLE :: projectedDistance(:,:),projectedXi(:,:),projectionVectors(:,:) INTEGER(INTG) :: elementIdx,elementNumber,lineFaceIdx,lineFaceNumber,computationalNodeIdx,xiIdx,localElementNumber, & @@ -1673,8 +1785,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 +1838,8 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection projectedDistance(2,:)=myComputationalNode !Find the globally closest distance in the current domain DO dataPointIdx=1,numberOfDataPoints + projectionExitTag(dataPointIdx)=dataProjection%dataProjectionResults(dataPointIdx)%exitTag + projectedLocalElement(dataPointIdx)=dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber CALL Sorting_BubbleIndexSort(globalClosestDistances(dataPointIdx,:),sortingIndices1,err,error,*999) sortingIndices1(1:totalNumberOfClosestCandidates)=sortingIndices1(1:totalNumberOfClosestCandidates)- & & globalMPIDisplacements(myComputationalNode+1) !shift the index to current computational node @@ -1747,11 +1863,11 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonLinesEvaluate(dataProjection,interpolatedPoint, & & dataPoints%dataPoints(dataPointIdx)%position,closestElements( & & dataPointIdx,1:numberOfClosestCandidates),closestLinesFaces(dataPointIdx,1: & - & numberOfClosestCandidates),projectionExitTag(dataPointIdx),projectedElement(dataPointIdx), & + & numberOfClosestCandidates),projectionExitTag(dataPointIdx),projectedLocalElement(dataPointIdx), & & projectedLineFace(dataPointIdx),projectedDistance(1,dataPointIdx),projectedXi(:,dataPointIdx), & & projectionVectors(:,dataPointIdx),err,error,*999) !Map the element number to global number - projectedElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedElement(dataPointIdx)) + projectedGlobalElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedLocalElement(dataPointIdx)) ENDIF ENDDO !dataPointIdx CASE(DATA_PROJECTION_BOUNDARY_FACES_PROJECTION_TYPE) @@ -1762,11 +1878,11 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonFacesEvaluate(dataProjection,interpolatedPoint, & & dataPoints%dataPoints(dataPointIdx)%position,closestElements( & & dataPointIdx,1:numberOfClosestCandidates),closestLinesFaces(dataPointIdx, & - & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx),projectedElement(dataPointIdx), & + & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx),projectedLocalElement(dataPointIdx), & & projectedLineFace(dataPointIdx),projectedDistance(1,dataPointIdx),projectedXi(:,dataPointIdx), & & projectionVectors(:,dataPointIdx),err,error,*999) !Map the element number to global number - projectedElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedElement(dataPointIdx)) + projectedGlobalElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedLocalElement(dataPointIdx)) ENDIF ENDDO !dataPointIdx CASE(DATA_PROJECTION_ALL_ELEMENTS_PROJECTION_TYPE) @@ -1779,10 +1895,10 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx, & & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx), & - & projectedElement(dataPointIdx),projectedDistance(1,dataPointIdx), & + & projectedLocalElement(dataPointIdx),projectedDistance(1,dataPointIdx), & & projectedXi(:,dataPointIdx),projectionVectors(:,dataPointIdx),err,error,*999) !Map the element number to global number - projectedElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedElement(dataPointIdx)) + projectedGlobalElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedLocalElement(dataPointIdx)) ENDIF ENDDO !dataPointIdx CASE(2) !2D element @@ -1792,11 +1908,11 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonElementsEvaluate_2(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx, & & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx), & - & projectedElement(dataPointIdx),projectedDistance(1,dataPointIdx), & + & projectedLocalElement(dataPointIdx),projectedDistance(1,dataPointIdx), & & projectedXi(:,dataPointIdx),projectionVectors(:,dataPointIdx), & & err,error,*999) !Map the element number to global number - projectedElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedElement(dataPointIdx)) + projectedGlobalElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedLocalElement(dataPointIdx)) ENDIF ENDDO !dataPointIdx CASE(3) !3D element @@ -1806,10 +1922,10 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonElementsEvaluate_3(dataProjection,interpolatedPoint,dataPoints% & & dataPoints(dataPointIdx)%position,closestElements(dataPointIdx, & & 1:numberOfClosestCandidates),projectionExitTag(dataPointIdx), & - & projectedElement(dataPointIdx),projectedDistance(1,dataPointIdx), & + & projectedLocalElement(dataPointIdx),projectedDistance(1,dataPointIdx), & & projectedXi(:,dataPointIdx),projectionVectors(:,dataPointIdx),err,error,*999) !Map the element number to global number - projectedElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedElement(dataPointIdx)) + projectedGlobalElement(dataPointIdx)=domainMappingElements%LOCAL_TO_GLOBAL_MAP(projectedLocalElement(dataPointIdx)) ENDIF ENDDO !dataPointIdx CASE DEFAULT @@ -1841,8 +1957,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 +1991,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 +2003,8 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection ENDDO !dataPointIdx projectedXi(:,sortingIndices2)=projectedXi projectionVectors(:, sortingIndices2)=projectionVectors - projectedElement(sortingIndices2)=projectedElement + projectedLocalElement(sortingIndices2)=projectedLocalElement + projectedGlobalElement(sortingIndices2)=projectedGlobalElement IF(dataProjection%projectionType==DATA_PROJECTION_BOUNDARY_LINES_PROJECTION_TYPE) THEN DO dataPointIdx=1,numberOfDataPoints dataProjection%dataProjectionResults(sortingIndices2(dataPointIdx))%elementLineFaceNumber=projectedLineFace(dataPointIdx) @@ -1899,10 +2023,12 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonLinesEvaluate(dataProjection,interpolatedPoint,dataPoints%dataPoints( & & dataPointIdx)%position,closestElements(dataPointIdx,:),closestLinesFaces(dataPointIdx,:), & & dataProjection%dataProjectionResults(dataPointIdx)%exitTag,dataProjection% & - & dataProjectionResults(dataPointIdx)%elementNumber,dataProjection% & + & dataProjectionResults(dataPointIdx)%elementLocalNumber,dataProjection% & & dataProjectionResults(dataPointIdx)%elementLineFaceNumber,dataProjection%dataProjectionResults( & - & dataPointIdx)%DISTANCE,dataProjection%dataProjectionResults(dataPointIdx)%xi, & + & dataPointIdx)%distance,dataProjection%dataProjectionResults(dataPointIdx)%xi, & & dataProjection%dataProjectionResults(dataPointIdx)%projectionVector,err,error,*999) + dataProjection%dataProjectionResults(dataPointIdx)%elementGlobalNumber= & + & dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber ENDDO !dataPointIdx CASE(DATA_PROJECTION_BOUNDARY_FACES_PROJECTION_TYPE) !Newton project to closest faces, and find miminum projection @@ -1910,10 +2036,12 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL DataProjection_NewtonFacesEvaluate(dataProjection,interpolatedPoint,dataPoints%dataPoints( & & dataPointIdx)%position,closestElements(dataPointIdx,:),closestLinesFaces(dataPointIdx,:), & & dataProjection%dataProjectionResults(dataPointIdx)%exitTag,dataProjection% & - & dataProjectionResults(dataPointIdx)%elementNumber,dataProjection% & + & dataProjectionResults(dataPointIdx)%elementLocalNumber,dataProjection% & & dataProjectionResults(dataPointIdx)%elementLineFaceNumber,dataProjection%dataProjectionResults( & & dataPointIdx)%DISTANCE,dataProjection%dataProjectionResults(dataPointIdx)%xi, & & dataProjection%dataProjectionResults(dataPointIdx)%projectionVector,err,error,*999) + dataProjection%dataProjectionResults(dataPointIdx)%elementGlobalNumber= & + & dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber ENDDO !dataPointIdx CASE(DATA_PROJECTION_ALL_ELEMENTS_PROJECTION_TYPE) !Newton project to closest elements, and find miminum projection @@ -1923,30 +2051,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 +2101,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 +2125,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 +2147,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 +2208,7 @@ SUBROUTINE DataProjection_DataPointsPositionEvaluate(dataProjection,field,fieldV interpolatedPoint=>interpolatedPoints(fieldVariableType)%ptr !Loop through data points DO dataPointIdx=1,dataPoints%numberOfDataPoints - elementNumber=dataProjection%dataProjectionResults(dataPointIdx)%elementNumber + elementNumber=dataProjection%dataProjectionResults(dataPointIdx)%elementLocalNumber CALL Field_InterpolationParametersElementGet(fieldParameterSetType,elementNumber, & & interpolationParameters(fieldVariableType)%ptr,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) CALL Field_InterpolateXi(NO_PART_DERIV,dataProjection%dataProjectionResults(dataPointIdx)%xi, & @@ -2253,118 +2389,128 @@ SUBROUTINE DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPo IF(.NOT.ASSOCIATED(dataProjection)) CALL FlagError("Data projection is not associated.",err,error,*999) IF(.NOT.dataProjection%dataProjectionFinished) CALL FlagError("Data projection has not been finished.",err,error,*999) - projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT numberOfCoordinates=dataProjection%numberOfCoordinates - relativeTolerance=dataProjection%relativeTolerance - absoluteTolerance=dataProjection%absoluteTolerance - maximumDelta=dataProjection%maximumIterationUpdate - minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small - !Project on each candidate elements - DO elementIdx=1,SIZE(candidateElements,1) - elementNumber=candidateElements(elementIdx) - IF(elementNumber>0) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT - converged=.FALSE. - !Start at half the maximumDelta as we do not know if quadratic model is a good approximation yet - delta=0.5_DP*maximumDelta - CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,elementNumber,interpolatedPoint% & - & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - xi=dataProjection%startingXi - CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector(1:numberOfCoordinates)=dataPointLocation-interpolatedPoint%values(:,NO_PART_DERIV) - functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations !(outer loop) - !Check for bounds [0,1] - IF(ABS(xi(1))absoluteTolerance) THEN !positive: minimum exists - updateXi(1)=-functionGradient/functionHessian - updateXiNorm=DABS(updateXi(1)) - insideRegion=updateXiNorm<=delta - ENDIF !positive - IF(.NOT.insideRegion) THEN !minimum not in the region - updateXi(1)=-DSIGN(delta,functionGradient) - updateXiNorm=delta - ENDIF - IF((bound/=0).AND.(bound>0.EQV.updateXi(1)>0.0_DP)) THEN !projection go out of element bound - exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS - EXIT main_loop - ENDIF - converged=updateXiNorm1.0_DP) THEN - newXi(1)=1.0_DP + IF(projectionExitTag==DATA_PROJECTION_USER_SPECIFIED) THEN + IF(projectionElementNumber==0) CALL FlagError("The projection element number has not been set.",err,error,*999) + CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,projectionElementNumber, & + & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + CALL Field_InterpolateXi(NO_PART_DERIV,projectionXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + projectionVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + projectionDistance=DSQRT(DOT_PRODUCT(projectionVector(1:numberOfCoordinates),projectionVector(1:numberOfCoordinates))) + ELSE + projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + relativeTolerance=dataProjection%relativeTolerance + absoluteTolerance=dataProjection%absoluteTolerance + maximumDelta=dataProjection%maximumIterationUpdate + minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small + !Project on each candidate elements + DO elementIdx=1,SIZE(candidateElements,1) + elementNumber=candidateElements(elementIdx) + IF(elementNumber>0) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + converged=.FALSE. + !Start at half the maximumDelta as we do not know if quadratic model is a good approximation yet + delta=0.5_DP*maximumDelta + CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,elementNumber,interpolatedPoint% & + & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + xi=dataProjection%startingXi + CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation-interpolatedPoint%values(:,NO_PART_DERIV) + functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations !(outer loop) + !Check for bounds [0,1] + IF(ABS(xi(1))absoluteTolerance) THEN !bad model: reduce step size - IF(delta<=minimumDelta) THEN - !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + !functionGradient + functionGradient=-2.0_DP* & + & (DOT_PRODUCT(distanceVector(1:numberOfCoordinates),interpolatedPoint%values(:,FIRST_PART_DERIV))) + !functionHessian + functionHessian=-2.0_DP* & + & (DOT_PRODUCT(distanceVector(1:numberOfCoordinates),interpolatedPoint%values(:,SECOND_PART_DERIV))- & + & DOT_PRODUCT(interpolatedPoint%values(:,FIRST_PART_DERIV),interpolatedPoint%values(:,FIRST_PART_DERIV))) + !A model trust region approach, directly taken from CMISS CLOS22: V = -(H + EIGEN_SHIFT*I)g + !The calculation of EIGEN_SHIFT are only approximated as opposed to the common trust region approach + !(inner loop: adjust region size) usually EXIT at 1 or 2 iterations + DO iterationIdx2=1,dataProjection%maximumNumberOfIterations + insideRegion=.FALSE. + IF(functionHessian>absoluteTolerance) THEN !positive: minimum exists + updateXi(1)=-functionGradient/functionHessian + updateXiNorm=DABS(updateXi(1)) + insideRegion=updateXiNorm<=delta + ENDIF !positive + IF(.NOT.insideRegion) THEN !minimum not in the region + updateXi(1)=-DSIGN(delta,functionGradient) + updateXiNorm=delta + ENDIF + IF((bound/=0).AND.(bound>0.EQV.updateXi(1)>0.0_DP)) THEN !projection go out of element bound + exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.25_DP*delta) - ELSE - predictedReduction=updateXi(1)*(functionGradient+0.5_DP*functionHessian*updateXi(1)) - predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction - IF(predictionAccuracy<0.01_DP) THEN !bad model: reduce region size + converged=updateXiNorm1.0_DP) THEN + newXi(1)=1.0_DP + ENDIF + CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector=dataPointLocation-interpolatedPoint%values(:,NO_PART_DERIV) + newFunctionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !second half of the convergence test + converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN !bad model: reduce step size IF(delta<=minimumDelta) THEN !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.5_DP*delta) - ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN - !good model: increase region size - delta=DMIN1(maximumDelta,2.0_DP*delta) - EXIT + delta=DMAX1(minimumDelta,0.25_DP*delta) ELSE - !ok model: keep the current region size - EXIT + predictedReduction=updateXi(1)*(functionGradient+0.5_DP*functionHessian*updateXi(1)) + predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction + IF(predictionAccuracy<0.01_DP) THEN !bad model: reduce region size + IF(delta<=minimumDelta) THEN + !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + EXIT main_loop + ENDIF + delta=DMAX1(minimumDelta,0.5_DP*delta) + ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN + !good model: increase region size + delta=DMIN1(maximumDelta,2.0_DP*delta) + EXIT + ELSE + !ok model: keep the current region size + EXIT + ENDIF ENDIF + ENDDO !iterationIdx2 (inner loop: adjust region size) + functionValue=newFunctionValue + xi=newXi + IF(converged) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED + EXIT ENDIF - ENDDO !iterationIdx2 (inner loop: adjust region size) - functionValue=newFunctionValue - xi=newXi - IF(converged) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED - EXIT + ENDDO main_loop !iterationIdx1 (outer loop) + IF(exitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT.AND.iterationIdx1>=dataProjection%maximumNumberOfIterations) & + & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION + IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)=dataProjection%maximumNumberOfIterations) & - & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION - IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%domain(meshComponentNumber)%ptr% & & mappings%elements numberOfCoordinates=dataProjection%numberOfCoordinates - relativeTolerance=dataProjection%relativeTolerance - absoluteTolerance=dataProjection%absoluteTolerance - maximumDelta=dataProjection%maximumIterationUpdate - minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small - DO elementIdx=1,SIZE(candidateElements,1) !project on each candidate elements - elementNumber=candidateElements(elementIdx) - IF(elementNumber>0) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT - converged=.FALSE. - !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet - delta=0.5_DP*maximumDelta - CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,elementNumber, & - & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - xi=dataProjection%startingXi - CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & - & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - !Outer loop - main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations - !Check for bounds [0,1] - DO xiIdx=1,2 - IF(ABS(xi(xiIdx))=temp1).AND.(minEigen>absoluteTolerance) - IF(insideRegion) THEN - determinant=minEigen*maxEigen !det(H) - hessianDiagonal(1)=functionHessian(1,1) - hessianDiagonal(2)=functionHessian(2,2) - ELSE - eigenShift=MAX(temp1-minEigen,absoluteTolerance) !shift towards steepest decent - determinant=temp1*(maxEigen+eigenShift) !det(H) - hessianDiagonal(1)=functionHessian(1,1)+eigenShift - hessianDiagonal(2)=functionHessian(2,2)+eigenShift - ENDIF - updateXi(1)=-(hessianDiagonal(2)*functionGradient(1)-functionHessian(1,2)*functionGradient(2))/determinant - updateXi(2)=(functionHessian(1,2)*functionGradient(1)-hessianDiagonal(1)*functionGradient(2))/determinant - updateXiNorm=DSQRT(DOT_PRODUCT(updateXi,updateXi)) - free=.TRUE. - DO xiIdx=1,2 - IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN - !Projection will go out of element bounds - IF(.NOT.free) THEN !both xi are fixed - exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS - EXIT main_loop - ENDIF - free=.FALSE. - fixedXiIdx=xiIdx - ENDIF - ENDDO !xiIdx - IF(free) THEN - !Both xi are free - IF(.NOT.insideRegion) THEN - IF(updateXiNorm>0.0_DP) THEN - updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound - ENDIF - ENDIF - ELSE - !xi are not free - updateXi(fixedXiIdx)=0.0_DP - xiIdx=3-fixedXiIdx - insideRegion=.FALSE. - IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN - !Positive: minimum exists in the unbounded direction - updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) - updateXiNorm=DABS(updateXi(xiIdx)) - insideRegion=updateXiNorm<=delta - ENDIF - IF(.NOT.insideRegion) THEN - !Minimum not in the region - updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) - updateXiNorm=delta + IF(projectionExitTag==DATA_PROJECTION_USER_SPECIFIED) THEN + IF(projectionElementNumber==0) CALL FlagError("The projection element number has not been set.",err,error,*999) + CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,projectionElementNumber, & + & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + CALL Field_InterpolateXi(NO_PART_DERIV,projectionXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + projectionVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + projectionDistance=DSQRT(DOT_PRODUCT(projectionVector(1:numberOfCoordinates),projectionVector(1:numberOfCoordinates))) + ELSE + projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + relativeTolerance=dataProjection%relativeTolerance + absoluteTolerance=dataProjection%absoluteTolerance + maximumDelta=dataProjection%maximumIterationUpdate + minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small + DO elementIdx=1,SIZE(candidateElements,1) !project on each candidate elements + elementNumber=candidateElements(elementIdx) + IF(elementNumber>0) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + converged=.FALSE. + !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet + delta=0.5_DP*maximumDelta + CALL Field_InterpolationParametersElementGet(dataProjection%projectionSetType,elementNumber, & + & interpolatedPoint%INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + xi=dataProjection%startingXi + CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Outer loop + main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations + !Check for bounds [0,1] + DO xiIdx=1,2 + IF(ABS(xi(xiIdx))1.0_DP) THEN - newXi(xiIdx)=1.0_DP - newXi(3-xiIdx)=xi(3-xiIdx)+updateXi(3-xiIdx)*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) + ENDDO !xiIdx + !functionGradient + functionGradient(1)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1))) + functionGradient(2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + !functionHessian + functionHessian(1,1)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1_S1))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1))) + functionHessian(1,2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1_S2))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + functionHessian(2,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=dataProjection%startingXi + CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + functionValue=DOT_PRODUCT(distanceVector,distanceVector) + !Outer loop + main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations + !Check for bounds [0,1] DO xiIdx=1,3 - IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN - !Projection will go out of element bounds - nBound=nBound+1 - free=.FALSE. - IF(nBound<=2) THEN - fixedXiIdx2(nBound)=xiIdx - ELSE - !All xi are fixed - exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS - EXIT main_loop - ENDIF - ENDIF - ENDDO !xiIdx - IF(free) THEN - !All xi are free - IF(.NOT.insideRegion) THEN - IF(updateXiNorm>0.0_DP) THEN - updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound - ENDIF + IF(ABS(xi(xiIdx))-1.0E-5_DP) THEN !include some negatives for numerical errors + minEigen=-temp1 !all eigenvalues are the same ELSE - !At least one of the xi are not free - !Try 2D projection - free=.TRUE. - fixedXiIdx=fixedXiIdx2(1) - IF(nBound==2) THEN - !only fix the direction that is most strongly suggesting leaving the element - IF(updateXi(fixedXiIdx2(2))>updateXi(fixedXiIdx2(1))) fixedXiIdx=fixedXiIdx2(2) - ENDIF - updateXi(fixedXiIdx)=0.0_DP - faceXiIdxs(1)=1+MOD(fixedXiIdx,3) - faceXiIdxs(2)=1+MOD(fixedXiIdx+1,3) - !functionGradient2 - functionGradient2(1)=functionGradient(faceXiIdxs(1)) - functionGradient2(2)=functionGradient(faceXiIdxs(2)) - !functionHessian2 - functionHessian2(1,1)=functionHessian(faceXiIdxs(1),faceXiIdxs(1)) - functionHessian2(1,2)=functionHessian(faceXiIdxs(1),faceXiIdxs(2)) - functionHessian2(2,2)=functionHessian(faceXiIdxs(2),faceXiIdxs(2)) - !re-estimate the trust solution in 2D - temp1=0.5_DP*(functionHessian2(1,1)+functionHessian2(2,2)) - temp2=DSQRT((0.5_DP*(functionHessian2(1,1)-functionHessian2(2,2)))**2+functionHessian2(1,2)**2) - minEigen=temp1-temp2 - maxEigen=temp1+temp2 - temp3=DSQRT(DOT_PRODUCT(functionGradient2,functionGradient2))/delta + temp3=DSQRT(-temp3) + temp4=(determinant+3.0_DP*(temp1*temp2)-2.0_DP*temp1**3)/(2.0_DP*temp3**3) + minEigen=2.0_DP*temp3*DCOS((DACOS(temp4)+TWOPI)/3.0_DP)-temp1 + ENDIF + functionGradientNorm=DSQRT(DOT_PRODUCT(functionGradient,functionGradient)) + !Inner loop: adjust region size, usually EXIT at 1 or 2 iterations + DO iterationIdx2=1,dataProjection%maximumNumberOfIterations + temp1=functionGradientNorm/delta !Estimate if the solution is inside the trust region without calculating a Newton step, this also guarantees !the Hessian matrix is positive definite - insideRegion=(minEigen>=temp3).AND.(minEigen>absoluteTolerance) + insideRegion=(minEigen>=temp1).AND.(minEigen>absoluteTolerance) IF(insideRegion) THEN - determinant=minEigen*maxEigen !determinant of functionHessian - hessianDiagonal2(1)=functionHessian2(1,1) - hessianDiagonal2(2)=functionHessian2(2,2) + hessianDiagonal(1)=functionHessian(1,1) + hessianDiagonal(2)=functionHessian(2,2) + hessianDiagonal(3)=functionHessian(3,3) ELSE - eigenShift=MAX(temp3-minEigen,absoluteTolerance) !shift towards steepest decent - determinant=temp3*(maxEigen+eigenShift) !determinant of shifted functionHessian - hessianDiagonal2(1)=functionHessian2(1,1)+eigenShift - hessianDiagonal2(2)=functionHessian2(2,2)+eigenShift + eigenShift=MAX(temp1-minEigen,absoluteTolerance) !shift towards steepest decent + determinant=determinant+eigenShift*(trace2+eigenShift*(trace+eigenShift)) !shift the determinant + hessianDiagonal(1)=functionHessian(1,1)+eigenShift + hessianDiagonal(2)=functionHessian(2,2)+eigenShift + hessianDiagonal(3)=functionHessian(3,3)+eigenShift ENDIF - updateXi(faceXiIdxs(1))=-(hessianDiagonal2(2)*functionGradient2(1)- & - & functionHessian2(1,2)*functionGradient2(2))/determinant - updateXi(faceXiIdxs(2))=(functionHessian2(1,2)*functionGradient2(1)- & - & hessianDiagonal2(1)*functionGradient2(2))/determinant + temp2=functionHessian(1,3)*functionHessian(2,3)-functionHessian(1,2)*hessianDiagonal(3) + temp3=functionHessian(1,2)*functionHessian(2,3)-functionHessian(1,3)*hessianDiagonal(2) + temp4=functionHessian(1,2)*functionHessian(1,3)-functionHessian(2,3)*hessianDiagonal(1) + updateXi(1)=((functionHessian(2,3)**2-hessianDiagonal(2)*hessianDiagonal(3))*functionGradient(1)- & + & temp2*functionGradient(2)-temp3*functionGradient(3))/determinant + updateXi(2)=((functionHessian(1,3)**2-hessianDiagonal(1)*hessianDiagonal(3))*functionGradient(2)- & + & temp2*functionGradient(1)-temp4*functionGradient(3))/determinant + updateXi(3)=((functionHessian(1,2)**2-hessianDiagonal(1)*hessianDiagonal(2))*functionGradient(3)- & + & temp3*functionGradient(1)-temp4*functionGradient(2))/determinant updateXiNorm=DSQRT(DOT_PRODUCT(updateXi,updateXi)) - !Check again for bounds - DO boundXiIdx=1,2 - IF((bound(faceXiIdxs(boundXiIdx))/=0).AND.(bound(faceXiIdxs(boundXiIdx))>0.EQV. & - & updateXi(faceXiIdxs(boundXiIdx))>0.0_DP)) THEN + free=.TRUE. + nBound=0 + DO xiIdx=1,3 + IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN !Projection will go out of element bounds - IF(.NOT.free) THEN - !Both xi are fixed + nBound=nBound+1 + free=.FALSE. + IF(nBound<=2) THEN + fixedXiIdx2(nBound)=xiIdx + ELSE + !All xi are fixed exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS EXIT main_loop ENDIF - free=.FALSE. - fixedBoundXiIdx=boundXiIdx ENDIF ENDDO !xiIdx IF(free) THEN - !Both xi are free + !All xi are free IF(.NOT.insideRegion) THEN IF(updateXiNorm>0.0_DP) THEN updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound ENDIF ENDIF ELSE - !xi are not free - updateXi(faceXiIdxs(fixedBoundXiIdx))=0.0_DP - xiIdx=faceXiIdxs(3-fixedBoundXiIdx) - insideRegion=.FALSE. - IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN - !positive: minimum exists in the unbounded direction - updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) - updateXiNorm=DABS(updateXi(xiIdx)) - insideRegion=updateXiNorm<=delta + !At least one of the xi are not free + !Try 2D projection + free=.TRUE. + fixedXiIdx=fixedXiIdx2(1) + IF(nBound==2) THEN + !only fix the direction that is most strongly suggesting leaving the element + IF(updateXi(fixedXiIdx2(2))>updateXi(fixedXiIdx2(1))) fixedXiIdx=fixedXiIdx2(2) ENDIF - IF(.NOT.insideRegion) THEN - !minimum not in the region - updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) - updateXiNorm=delta + updateXi(fixedXiIdx)=0.0_DP + faceXiIdxs(1)=1+MOD(fixedXiIdx,3) + faceXiIdxs(2)=1+MOD(fixedXiIdx+1,3) + !functionGradient2 + functionGradient2(1)=functionGradient(faceXiIdxs(1)) + functionGradient2(2)=functionGradient(faceXiIdxs(2)) + !functionHessian2 + functionHessian2(1,1)=functionHessian(faceXiIdxs(1),faceXiIdxs(1)) + functionHessian2(1,2)=functionHessian(faceXiIdxs(1),faceXiIdxs(2)) + functionHessian2(2,2)=functionHessian(faceXiIdxs(2),faceXiIdxs(2)) + !re-estimate the trust solution in 2D + temp1=0.5_DP*(functionHessian2(1,1)+functionHessian2(2,2)) + temp2=DSQRT((0.5_DP*(functionHessian2(1,1)-functionHessian2(2,2)))**2+functionHessian2(1,2)**2) + minEigen=temp1-temp2 + maxEigen=temp1+temp2 + temp3=DSQRT(DOT_PRODUCT(functionGradient2,functionGradient2))/delta + !Estimate if the solution is inside the trust region without calculating a Newton step, this also guarantees + !the Hessian matrix is positive definite + insideRegion=(minEigen>=temp3).AND.(minEigen>absoluteTolerance) + IF(insideRegion) THEN + determinant=minEigen*maxEigen !determinant of functionHessian + hessianDiagonal2(1)=functionHessian2(1,1) + hessianDiagonal2(2)=functionHessian2(2,2) + ELSE + eigenShift=MAX(temp3-minEigen,absoluteTolerance) !shift towards steepest decent + determinant=temp3*(maxEigen+eigenShift) !determinant of shifted functionHessian + hessianDiagonal2(1)=functionHessian2(1,1)+eigenShift + hessianDiagonal2(2)=functionHessian2(2,2)+eigenShift ENDIF - ENDIF !if xi are free (2D) - ENDIF !if xi are free (3D) - !First half of the convergence test - converged=updateXiNorm0.EQV. & + & updateXi(faceXiIdxs(boundXiIdx))>0.0_DP)) THEN + !Projection will go out of element bounds + IF(.NOT.free) THEN + !Both xi are fixed + exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS + EXIT main_loop ENDIF - ENDDO - ELSEIF(newXi(xiIdx)>1.0_DP) THEN - newXi(xiIdx)=1.0_DP - DO xiIdx2 = 1,3 - IF(xiIdx2 /= xiIdx) THEN - newXi(xiIdx2)=xi(xiIdx2)+updateXi(xiIdx2) + free=.FALSE. + fixedBoundXiIdx=boundXiIdx + ENDIF + ENDDO !xiIdx + IF(free) THEN + !Both xi are free + IF(.NOT.insideRegion) THEN + IF(updateXiNorm>0.0_DP) THEN + updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound ENDIF - ENDDO + ENDIF + ELSE + !xi are not free + updateXi(faceXiIdxs(fixedBoundXiIdx))=0.0_DP + xiIdx=faceXiIdxs(3-fixedBoundXiIdx) + insideRegion=.FALSE. + IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN + !positive: minimum exists in the unbounded direction + updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) + updateXiNorm=DABS(updateXi(xiIdx)) + insideRegion=updateXiNorm<=delta + ENDIF + IF(.NOT.insideRegion) THEN + !minimum not in the region + updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) + updateXiNorm=delta + ENDIF + ENDIF !if xi are free (2D) + ENDIF !if xi are free (3D) + !First half of the convergence test + converged=updateXiNorm1.0_DP) THEN + newXi(xiIdx)=1.0_DP + DO xiIdx2 = 1,3 + IF(xiIdx2 /= xiIdx) THEN + newXi(xiIdx2)=xi(xiIdx2)+updateXi(xiIdx2) + ENDIF + ENDDO + ENDIF + ELSE IF(newXi(xiIdx)<0.0_DP) THEN + !boundary collision check + newXi(xiIdx)=0.0_DP + newXi=xi-updateXi*xi(xiIdx)/updateXi(xiIdx) + ELSE IF(newXi(xiIdx)>1.0_DP) THEN + newXi(xiIdx)=1.0_DP + newXi=xi+updateXi*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) ENDIF - ELSE IF(newXi(xiIdx)<0.0_DP) THEN - !boundary collision check - newXi(xiIdx)=0.0_DP - newXi=xi-updateXi*xi(xiIdx)/updateXi(xiIdx) - ELSE IF(newXi(xiIdx)>1.0_DP) THEN - newXi(xiIdx)=1.0_DP - newXi=xi+updateXi*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) - ENDIF - ENDDO !xiIdx - CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector=dataPointLocation-interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - newFunctionValue=DOT_PRODUCT(distanceVector,distanceVector) - !Second half of the convergence test (before collision detection) - converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN - !bad model: reduce step size - IF(delta<=minimumDelta) THEN - !something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION ! \it will get stucked!! - EXIT main_loop - ENDIF - delta=DMAX1(minimumDelta,0.25_DP*delta) - ELSE - predictedReduction=DOT_PRODUCT(functionGradient,updateXi)+ & - & 0.5_DP*(updateXi(1)*(updateXi(1)*functionHessian(1,1)+2.0_DP*updateXi(2)*functionHessian(1,2)+ & - & 2.0_DP*updateXi(3)*functionHessian(1,3))+updateXi(2)*(updateXi(2)*functionHessian(2,2)+ & - & 2.0_DP*updateXi(3)*functionHessian(2,3))+updateXi(2)**2*functionHessian(2,2)) - IF (ABS(predictedReduction) < ZERO_TOLERANCE) THEN - converged = .TRUE. !prediction reduction converged - EXIT - ENDIF - predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction - IF(predictionAccuracy<0.01_DP) THEN - !bad model: reduce region size + ENDDO !xiIdx + CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector=dataPointLocation-interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + newFunctionValue=DOT_PRODUCT(distanceVector,distanceVector) + !Second half of the convergence test (before collision detection) + converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN + !bad model: reduce step size IF(delta<=minimumDelta) THEN !something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION ! \it will get stucked!! EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.5_DP*delta) - ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN - !good model: increase region size - delta=DMIN1(maximumDelta,2.0_DP*delta) - EXIT + delta=DMAX1(minimumDelta,0.25_DP*delta) ELSE - !ok model: keep the current region size - EXIT + predictedReduction=DOT_PRODUCT(functionGradient,updateXi)+ & + & 0.5_DP*(updateXi(1)*(updateXi(1)*functionHessian(1,1)+2.0_DP*updateXi(2)*functionHessian(1,2)+ & + & 2.0_DP*updateXi(3)*functionHessian(1,3))+updateXi(2)*(updateXi(2)*functionHessian(2,2)+ & + & 2.0_DP*updateXi(3)*functionHessian(2,3))+updateXi(2)**2*functionHessian(2,2)) + IF (ABS(predictedReduction) < ZERO_TOLERANCE) THEN + converged = .TRUE. !prediction reduction converged + EXIT + ENDIF + predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction + IF(predictionAccuracy<0.01_DP) THEN + !bad model: reduce region size + IF(delta<=minimumDelta) THEN + !something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + EXIT main_loop + ENDIF + delta=DMAX1(minimumDelta,0.5_DP*delta) + ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN + !good model: increase region size + delta=DMIN1(maximumDelta,2.0_DP*delta) + EXIT + ELSE + !ok model: keep the current region size + EXIT + ENDIF ENDIF + ENDDO !iterationIdx2 (inner loop: adjust region size) + functionValue=newFunctionValue + xi=newXi + IF(converged) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED + EXIT ENDIF - ENDDO !iterationIdx2 (inner loop: adjust region size) - functionValue=newFunctionValue - xi=newXi - IF(converged) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED - EXIT + ENDDO main_loop !iterationIdx1 (outer loop) + IF(exitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT.AND.iterationIdx1>=dataProjection%maximumNumberOfIterations) & + & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION + IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)=dataProjection%maximumNumberOfIterations) & - & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION - IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)interpolatedPoint%INTERPOLATION_PARAMETERS%FIELD%decomposition%domain(meshComponentNumber)% & & PTR%MAPPINGS%ELEMENTS numberOfCoordinates=dataProjection%numberOfCoordinates - relativeTolerance=dataProjection%relativeTolerance - absoluteTolerance=dataProjection%absoluteTolerance - maximumDelta=dataProjection%maximumIterationUpdate - minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small - !Project on each candidate elements - DO elementIdx=1,SIZE(candidateElements,1) - elementNumber=candidateElements(elementIdx) - IF(elementNumber>0) THEN - elementFaceNumber=candidateElementFaces(elementIdx) - faceNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & - & elements(elementNumber)%ELEMENT_FACES(elementFaceNumber) - exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT - converged=.FALSE. - !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet - delta=0.5_DP*maximumDelta - CALL Field_InterpolationParametersFaceGet(dataProjection%projectionSetType,faceNumber,interpolatedPoint% & - & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - xi=dataProjection%startingXi - CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & - & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - !Outer loop - main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations - !Check for bounds [0,1] - DO xiIdx=1,2 - IF(ABS(xi(xiIdx))=temp1).AND.(minEigen>absoluteTolerance) - IF(insideRegion) THEN - determinant=minEigen*maxEigen !det(H) - hessianDiagonal(1)=functionHessian(1,1) - hessianDiagonal(2)=functionHessian(2,2) - ELSE - eigenShift=MAX(temp1-minEigen,absoluteTolerance) !shift towards steepest decent - determinant=temp1*(maxEigen+eigenShift) !det(H) - hessianDiagonal(1)=functionHessian(1,1)+eigenShift - hessianDiagonal(2)=functionHessian(2,2)+eigenShift - ENDIF - updateXi(1)=-(hessianDiagonal(2)*functionGradient(1)-functionHessian(1,2)*functionGradient(2))/determinant - updateXi(2)=(functionHessian(1,2)*functionGradient(1)-hessianDiagonal(1)*functionGradient(2))/determinant - updateXiNorm=DSQRT(DOT_PRODUCT(updateXi,updateXi)) - free=.TRUE. - DO xiIdx=1,2 - IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN - !Projection will go out of element bounds - IF(.NOT.free) THEN !both xi are fixed - exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS - EXIT main_loop - ENDIF - free=.FALSE. - fixedXiIdx=xiIdx - ENDIF - ENDDO !xiIdx - IF(free) THEN - !Both xi are free - IF(.NOT.insideRegion) THEN - IF(updateXiNorm>0.0_DP) THEN - updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound - ENDIF - ENDIF - ELSE - !xi are not free - updateXi(fixedXiIdx)=0.0_DP - xiIdx=3-fixedXiIdx - insideRegion=.FALSE. - IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN - !Positive: minimum exists in the unbounded direction - updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) - updateXiNorm=DABS(updateXi(xiIdx)) - insideRegion=updateXiNorm<=delta - ENDIF - IF(.NOT.insideRegion) THEN - !Minimum not in the region - updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) - updateXiNorm=delta + IF(projectionExitTag==DATA_PROJECTION_USER_SPECIFIED) THEN + IF(projectionElementNumber==0) CALL FlagError("The projection element number has not been set.",err,error,*999) + IF(projectionElementFaceNumber==0) CALL FlagError("The projection element face number has not been set.",err,error,*999) + faceNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & + & elements(projectionElementNumber)%ELEMENT_FACES(projectionElementFaceNumber) + CALL Field_InterpolationParametersFaceGet(dataProjection%projectionSetType,faceNumber,interpolatedPoint% & + & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + CALL Field_InterpolateXi(NO_PART_DERIV,projectionXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + projectionVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + projectionDistance=DSQRT(DOT_PRODUCT(projectionVector(1:numberOfCoordinates),projectionVector(1:numberOfCoordinates))) + ELSE + projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + relativeTolerance=dataProjection%relativeTolerance + absoluteTolerance=dataProjection%absoluteTolerance + maximumDelta=dataProjection%maximumIterationUpdate + minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small + !Project on each candidate elements + DO elementIdx=1,SIZE(candidateElements,1) + elementNumber=candidateElements(elementIdx) + IF(elementNumber>0) THEN + elementFaceNumber=candidateElementFaces(elementIdx) + faceNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & + & elements(elementNumber)%ELEMENT_FACES(elementFaceNumber) + exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + converged=.FALSE. + !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet + delta=0.5_DP*maximumDelta + CALL Field_InterpolationParametersFaceGet(dataProjection%projectionSetType,faceNumber,interpolatedPoint% & + & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + xi=dataProjection%startingXi + CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Outer loop + main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations + !Check for bounds [0,1] + DO xiIdx=1,2 + IF(ABS(xi(xiIdx))1.0_DP) THEN - newXi(xiIdx)=1.0_DP - newXi(3-xiIdx)=xi(3-xiIdx)+updateXi(3-xiIdx)*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) + ENDDO !xiIdx + !functionGradient + functionGradient(1)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1))) + functionGradient(2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + !functionHessian + functionHessian(1,1)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1_S1))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1))) + functionHessian(1,2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1_S2))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S1), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + functionHessian(2,2)=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2_S2))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2), & + & interpolatedPoint%values(1:numberOfCoordinates,PART_DERIV_S2))) + !A model trust region approach, a Newton step is taken if the minimum lies inside the trust region (delta), + !if not, shift the step towards the steepest descent + temp1=0.5_DP*(functionHessian(1,1)+functionHessian(2,2)) + temp2=DSQRT((0.5_DP*(functionHessian(1,1)-functionHessian(2,2)))**2+functionHessian(1,2)**2) + minEigen=temp1-temp2 + maxEigen=temp1+temp2 + functionGradientNorm=DSQRT(DOT_PRODUCT(functionGradient,functionGradient)) + !Inner loop: adjust region size, usually EXIT at 1 or 2 iterations + DO iterationIdx2=1,dataProjection%maximumNumberOfIterations + temp1=functionGradientNorm/delta + !Estimate if the solution is inside the trust region without calculating a Newton step, this also guarantees + !the Hessian matrix is positive definite + insideRegion=(minEigen>=temp1).AND.(minEigen>absoluteTolerance) + IF(insideRegion) THEN + determinant=minEigen*maxEigen !det(H) + hessianDiagonal(1)=functionHessian(1,1) + hessianDiagonal(2)=functionHessian(2,2) + ELSE + eigenShift=MAX(temp1-minEigen,absoluteTolerance) !shift towards steepest decent + determinant=temp1*(maxEigen+eigenShift) !det(H) + hessianDiagonal(1)=functionHessian(1,1)+eigenShift + hessianDiagonal(2)=functionHessian(2,2)+eigenShift ENDIF - ENDDO - CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector=dataPointLocation-interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - newFunctionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - !Second half of the convergence test (before collision detection) - converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN - !Bad model: reduce step size - IF(delta<=minimumDelta) THEN - !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + updateXi(1)=-(hessianDiagonal(2)*functionGradient(1)-functionHessian(1,2)*functionGradient(2))/determinant + updateXi(2)=(functionHessian(1,2)*functionGradient(1)-hessianDiagonal(1)*functionGradient(2))/determinant + updateXiNorm=DSQRT(DOT_PRODUCT(updateXi,updateXi)) + free=.TRUE. + DO xiIdx=1,2 + IF((bound(xiIdx)/=0).AND.(bound(xiIdx)>0.EQV.updateXi(xiIdx)>0.0_DP)) THEN + !Projection will go out of element bounds + IF(.NOT.free) THEN !both xi are fixed + exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS + EXIT main_loop + ENDIF + free=.FALSE. + fixedXiIdx=xiIdx + ENDIF + ENDDO !xiIdx + IF(free) THEN + !Both xi are free + IF(.NOT.insideRegion) THEN + IF(updateXiNorm>0.0_DP) THEN + updateXi=delta/updateXiNorm*updateXi !readjust updateXi to lie on the region bound + ENDIF + ENDIF + ELSE + !xi are not free + updateXi(fixedXiIdx)=0.0_DP + xiIdx=3-fixedXiIdx + insideRegion=.FALSE. + IF(functionHessian(xiIdx,xiIdx)>0.0_DP) THEN + !Positive: minimum exists in the unbounded direction + updateXi(xiIdx)=-functionGradient(xiIdx)/functionHessian(xiIdx,xiIdx) + updateXiNorm=DABS(updateXi(xiIdx)) + insideRegion=updateXiNorm<=delta + ENDIF + IF(.NOT.insideRegion) THEN + !Minimum not in the region + updateXi(xiIdx)=-DSIGN(delta,functionGradient(xiIdx)) + updateXiNorm=delta + ENDIF + ENDIF !if xi is free + !First half of the convergence test + converged=updateXiNorm1.0_DP) THEN + newXi(xiIdx)=1.0_DP + newXi(3-xiIdx)=xi(3-xiIdx)+updateXi(3-xiIdx)*(1.0_DP-xi(xiIdx))/updateXi(xiIdx) + ENDIF + ENDDO + CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector=dataPointLocation-interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + newFunctionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Second half of the convergence test (before collision detection) + converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN + !Bad model: reduce step size IF(delta<=minimumDelta) THEN !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! functionValue=newFunctionValue EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.5_DP*delta) - ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN - !Good model: increase region size - delta=DMIN1(maximumDelta,2.0_DP*delta) - functionValue=newFunctionValue - EXIT + delta=DMAX1(minimumDelta,0.25_DP*delta) ELSE - !OK model: keep the current region size - functionValue=newFunctionValue - EXIT + predictedReduction=DOT_PRODUCT(functionGradient,updateXi)+ & + & 0.5_DP*(updateXi(1)*(updateXi(1)*functionHessian(1,1)+2.0_DP*updateXi(2)*functionHessian(1,2))+ & + & updateXi(2)**2*functionHessian(2,2)) + predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction + IF(predictionAccuracy<0.01_DP) THEN + !Bad model: reduce region size + IF(delta<=minimumDelta) THEN + !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + functionValue=newFunctionValue + EXIT main_loop + ENDIF + delta=DMAX1(minimumDelta,0.5_DP*delta) + ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN + !Good model: increase region size + delta=DMIN1(maximumDelta,2.0_DP*delta) + functionValue=newFunctionValue + EXIT + ELSE + !OK model: keep the current region size + functionValue=newFunctionValue + EXIT + ENDIF ENDIF + ENDDO !iterationIdx2 (inner loop: adjust region size) + functionValue=newFunctionValue + xi=newXi + IF(converged) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED + EXIT ENDIF - ENDDO !iterationIdx2 (inner loop: adjust region size) - functionValue=newFunctionValue - xi=newXi - IF(converged) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED - EXIT + ENDDO main_loop !iterationIdx1 (outer loop) + IF(exitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT.AND.iterationIdx1>=dataProjection%maximumNumberOfIterations) & + & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION + IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)=dataProjection%maximumNumberOfIterations) & - & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION - IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)0) THEN - elementLineNumber=candidateElementLines(elementIdx) - lineNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & - & elements(elementNumber)%ELEMENT_LINES(elementLineNumber) - exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT - converged=.FALSE. - !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet - delta=0.5_DP*maximumDelta - CALL Field_InterpolationParametersLineGet(dataProjection%projectionSetType,lineNumber,interpolatedPoint% & - & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - xi=dataProjection%startingXi - CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & - & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) - functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) - !Outer loop - main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations - !Check for bounds [0,1] - IF(ABS(xi(1))absoluteTolerance) THEN - !Positive: minimum exists - updateXi(1)=-functionGradient/functionHessian - updateXiNorm=DABS(updateXi(1)) - insideRegion=updateXiNorm<=delta - ENDIF - IF(.NOT.insideRegion) THEN - !minimum not in the region - updateXi(1)=-DSIGN(delta,functionGradient) - updateXiNorm=delta - ENDIF - IF((bound/=0).AND.(bound>0.EQV.updateXi(1)>0.0_DP)) THEN - !Projection will go out of element bounds - exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS - EXIT main_loop - ENDIF - !First half of the convergence test (before collision detection) - converged=updateXiNorm1.0_DP) THEN - newXi(1)=1.0_DP + IF(projectionExitTag==DATA_PROJECTION_USER_SPECIFIED) THEN + IF(projectionElementNumber==0) CALL FlagError("The projection element number has not been set.",err,error,*999) + IF(projectionElementLineNumber==0) CALL FlagError("The projection element line number has not been set.",err,error,*999) + lineNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & + & elements(projectionElementNumber)%ELEMENT_LINES(projectionElementLineNumber) + CALL Field_InterpolationParametersLineGet(dataProjection%projectionSetType,lineNumber,interpolatedPoint% & + & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + CALL Field_InterpolateXi(NO_PART_DERIV,projectionXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + projectionVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + projectionDistance=DSQRT(DOT_PRODUCT(projectionVector(1:numberOfCoordinates),projectionVector(1:numberOfCoordinates))) + ELSE + projectionExitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + relativeTolerance=dataProjection%relativeTolerance + absoluteTolerance=dataProjection%absoluteTolerance + maximumDelta=dataProjection%maximumIterationUpdate + minimumDelta=0.025_DP*maximumDelta !need to set a minimum, in case if it gets too small + !Project on each candidate elements + DO elementIdx=1,SIZE(candidateElements,1) + elementNumber=candidateElements(elementIdx) + IF(elementNumber>0) THEN + elementLineNumber=candidateElementLines(elementIdx) + lineNumber=interpolatedPoint%INTERPOLATION_PARAMETERS%field%decomposition%topology%elements% & + & elements(elementNumber)%ELEMENT_LINES(elementLineNumber) + exitTag=DATA_PROJECTION_EXIT_TAG_NO_ELEMENT + converged=.FALSE. + !start at half the maximumDelta as we do not know if quadratic model is a good approximation yet + delta=0.5_DP*maximumDelta + CALL Field_InterpolationParametersLineGet(dataProjection%projectionSetType,lineNumber,interpolatedPoint% & + & INTERPOLATION_PARAMETERS,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + xi=dataProjection%startingXi + CALL Field_InterpolateXi(SECOND_PART_DERIV,xi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + functionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Outer loop + main_loop: DO iterationIdx1=1,dataProjection%maximumNumberOfIterations + !Check for bounds [0,1] + IF(ABS(xi(1))absoluteTolerance) THEN - !Bad model: reduce step size - IF(delta<=minimumDelta) THEN - !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small - exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + !functionGradient + functionGradient=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,FIRST_PART_DERIV))) + !functionHessian + functionHessian=-2.0_DP*(DOT_PRODUCT(distanceVector(1:numberOfCoordinates), & + & interpolatedPoint%values(1:numberOfCoordinates,SECOND_PART_DERIV))- & + & DOT_PRODUCT(interpolatedPoint%values(1:numberOfCoordinates,FIRST_PART_DERIV), & + & interpolatedPoint%values(1:numberOfCoordinates,FIRST_PART_DERIV))) + !A model trust region approach, directly taken from CMISS CLOS22: V = -(H + eigenShift*I)g + !The calculation of eigenShift are only approximated as oppose to the common trust region approach + !Inner loop: adjust region size, usually EXIT at 1 or 2 iterations + DO iterationIdx2=1,dataProjection%maximumNumberOfIterations + insideRegion=.FALSE. + IF(functionHessian>absoluteTolerance) THEN + !Positive: minimum exists + updateXi(1)=-functionGradient/functionHessian + updateXiNorm=DABS(updateXi(1)) + insideRegion=updateXiNorm<=delta + ENDIF + IF(.NOT.insideRegion) THEN + !minimum not in the region + updateXi(1)=-DSIGN(delta,functionGradient) + updateXiNorm=delta + ENDIF + IF((bound/=0).AND.(bound>0.EQV.updateXi(1)>0.0_DP)) THEN + !Projection will go out of element bounds + exitTag=DATA_PROJECTION_EXIT_TAG_BOUNDS EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.25_DP*delta) - ELSE - predictedReduction=updateXi(1)*(functionGradient+0.5_DP*functionHessian*updateXi(1)) - predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction - IF(predictionAccuracy<0.01_DP) THEN - !Bad model: reduce region size + !First half of the convergence test (before collision detection) + converged=updateXiNorm1.0_DP) THEN + newXi(1)=1.0_DP + ENDIF + CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) + distanceVector(1:numberOfCoordinates)=dataPointLocation(1:numberOfCoordinates)- & + & interpolatedPoint%values(1:numberOfCoordinates,NO_PART_DERIV) + newFunctionValue=DOT_PRODUCT(distanceVector(1:numberOfCoordinates),distanceVector(1:numberOfCoordinates)) + !Second half of the convergence test + converged=converged.AND.(DABS(newFunctionValue-functionValue)/(1.0_DP+functionValue)absoluteTolerance) THEN + !Bad model: reduce step size IF(delta<=minimumDelta) THEN !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! EXIT main_loop ENDIF - delta=DMAX1(minimumDelta,0.5_DP*delta) - ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN - !Good model: increase region size - delta=DMIN1(maximumDelta,2.0_DP*delta) - EXIT + delta=DMAX1(minimumDelta,0.25_DP*delta) ELSE - !OK model: keep the current region size - EXIT + predictedReduction=updateXi(1)*(functionGradient+0.5_DP*functionHessian*updateXi(1)) + predictionAccuracy=(newFunctionValue-functionValue)/predictedReduction + IF(predictionAccuracy<0.01_DP) THEN + !Bad model: reduce region size + IF(delta<=minimumDelta) THEN + !Something went wrong, minimumDelta too large? not likely to happen if minimumDelta is small + exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION !it will get stucked!! + EXIT main_loop + ENDIF + delta=DMAX1(minimumDelta,0.5_DP*delta) + ELSEIF(predictionAccuracy>0.9_DP.AND.predictionAccuracy<1.1_DP) THEN + !Good model: increase region size + delta=DMIN1(maximumDelta,2.0_DP*delta) + EXIT + ELSE + !OK model: keep the current region size + EXIT + ENDIF ENDIF + ENDDO !iterationIdx2 (inner loop: adjust region size) + functionValue=newFunctionValue + xi=newXi + IF(converged) THEN + exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED + EXIT ENDIF - ENDDO !iterationIdx2 (inner loop: adjust region size) - functionValue=newFunctionValue - xi=newXi - IF(converged) THEN - exitTag=DATA_PROJECTION_EXIT_TAG_CONVERGED - EXIT + ENDDO main_loop !iterationIdx1 (outer loop) + IF(exitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT.AND.iterationIdx1>=dataProjection%maximumNumberOfIterations) & + & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION + IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)=dataProjection%maximumNumberOfIterations) & - & exitTag=DATA_PROJECTION_EXIT_TAG_MAX_ITERATION - IF((projectionExitTag==DATA_PROJECTION_EXIT_TAG_NO_ELEMENT).OR.(DSQRT(functionValue)distance) THEN - dataProjection%dataProjectionResults(dataPointIdx)%exitTag=DATA_PROJECTION_CANCELLED - dataProjection%dataProjectionResults(dataPointIdx)%elementNumber=0 - dataProjection%dataProjectionResults(dataPointIdx)%elementLineFaceNumber=0 - dataProjection%dataProjectionResults(dataPointIdx)%xi=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%elementXi=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%distance=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%projectionVector=0.0_DP + CALL DataProjection_DataProjectionResultCancel(dataProjection%dataProjectionResults(dataPointIdx),err,error,*999) ENDIF ENDDO !dataPointIdx CASE(DATA_PROJECTION_DISTANCE_GREATER_EQUAL) DO dataPointIdx=1,dataPoints%numberOfDataPoints IF(dataProjection%dataProjectionResults(dataPointIdx)%distance>=distance) THEN - dataProjection%dataProjectionResults(dataPointIdx)%exitTag=DATA_PROJECTION_CANCELLED - dataProjection%dataProjectionResults(dataPointIdx)%elementNumber=0 - dataProjection%dataProjectionResults(dataPointIdx)%elementLineFaceNumber=0 - dataProjection%dataProjectionResults(dataPointIdx)%xi=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%elementXi=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%distance=0.0_DP - dataProjection%dataProjectionResults(dataPointIdx)%projectionVector=0.0_DP + CALL DataProjection_DataProjectionResultCancel(dataProjection%dataProjectionResults(dataPointIdx),err,error,*999) ENDIF ENDDO !dataPointIdx CASE(DATA_PROJECTION_DISTANCE_LESS) DO dataPointIdx=1,dataPoints%numberOfDataPoints IF(dataProjection%dataProjectionResults(dataPointIdx)%distance=1) THEN @@ -4889,7 +5058,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 +5072,75 @@ 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 +5179,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 +5297,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 +5454,7 @@ END SUBROUTINE DataProjection_ResultExitTagGet ! !>Gets the projection xi for a data point identified by a given global number. - SUBROUTINE DataProjection_ResultXiGet(dataProjection,dataPointUserNumber,projectionXi,err,error,*) + SUBROUTINE DataProjection_ResultProjectionXiGet(dataProjection,dataPointUserNumber,projectionXi,err,error,*) !Argument variables TYPE(DataProjectionType), POINTER :: dataProjection !Sets the projection xi for a data point identified by a given global number. - SUBROUTINE DataProjection_ResultXiSet(dataProjection,dataPointUserNumber,projectionXi,err,error,*) + SUBROUTINE DataProjection_ResultProjectionXiSet(dataProjection,dataPointUserNumber,projectionXi,err,error,*) !Argument variables TYPE(DataProjectionType), POINTER :: dataProjection !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..46525954 100644 --- a/src/finite_elasticity_routines.F90 +++ b/src/finite_elasticity_routines.F90 @@ -3700,9 +3700,9 @@ SUBROUTINE FiniteElasticity_StressStrainCalculate(equationsSet,derivedType,field TYPE(VARYING_STRING), INTENT(OUT) :: error ! 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..510774e2 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 @@ -8554,6 +8594,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 +20196,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) @@ -22794,6 +22992,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 +23132,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 +23276,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 +23416,74 @@ END SUBROUTINE cmfe_DataProjection_ResultExitTagGetObj !================================================================================================================================ ! + !>Returns the element xi for a data point in a set of data points identified by user number. + SUBROUTINE cmfe_DataProjection_ResultElementXiGetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & + & dataPointUserNumber,elementXi,err) + !DLLEXPORT(cmfe_DataProjection_ResultElementXiGetNumber) + + !Argument variables + INTEGER(INTG), INTENT(IN) :: regionUserNumber !Returns the element xi for a data point in a set of data points identified by an object. + SUBROUTINE cmfe_DataProjection_ResultElementXiGetObj(dataProjection,dataPointUserNumber,elementXi,err) + !DLLEXPORT(cmfe_DataProjection_ResultElementXiGetObj) + + !Argument variables + TYPE(cmfe_DataProjectionType), INTENT(IN) :: dataProjection !Returns the maximum error for a data projection given by numbers. SUBROUTINE cmfe_DataProjection_ResultMaximumErrorGetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & & maximumDataPointUserNumber,maximumError,err) @@ -23206,9 +23687,9 @@ END SUBROUTINE cmfe_DataProjection_ResultRMSErrorGetObj ! !>Returns the projection xi for a data point in a set of data points identified by user number. - SUBROUTINE cmfe_DataProjection_ResultXiGetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & + SUBROUTINE cmfe_DataProjection_ResultProjectionXiGetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & & dataPointUserNumber,projectionXi,err) - !DLLEXPORT(cmfe_DataProjection_ResultXiGetNumber) + !DLLEXPORT(cmfe_DataProjection_ResultProjectionXiGetNumber) !Argument variables INTEGER(INTG), INTENT(IN) :: regionUserNumber !Returns the projection xi for a data point in a set of data points identified by an object. - SUBROUTINE cmfe_DataProjection_ResultXiGetObj(dataProjection,dataPointUserNumber,projectionXi,err) - !DLLEXPORT(cmfe_DataProjection_ResultXiGetObj) + SUBROUTINE cmfe_DataProjection_ResultProjectionXiGetObj(dataProjection,dataPointUserNumber,projectionXi,err) + !DLLEXPORT(cmfe_DataProjection_ResultProjectionXiGetObj) !Argument variables TYPE(cmfe_DataProjectionType), INTENT(IN) :: dataProjection !Sets the projection xi for a data point in a set of data points identified by user number. - SUBROUTINE cmfe_DataProjection_ResultXiSetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & + SUBROUTINE cmfe_DataProjection_ResultProjectionXiSetNumber(regionUserNumber,dataPointsUserNumber,dataProjectionUserNumber, & & dataPointUserNumber,projectionXi,err) - !DLLEXPORT(cmfe_DataProjection_ResultXiSetNumber) + !DLLEXPORT(cmfe_DataProjection_ResultProjectionXiSetNumber) !Argument variables INTEGER(INTG), INTENT(IN) :: regionUserNumber !Sets the projection xi for a data point in a set of data points identified by an object. - SUBROUTINE cmfe_DataProjection_ResultXiSetObj(dataProjection,dataPointUserNumber,projectionXi,err) - !DLLEXPORT(cmfe_DataProjection_ResultXiSetObj) + SUBROUTINE cmfe_DataProjection_ResultProjectionXiSetObj(dataProjection,dataPointUserNumber,projectionXi,err) + !DLLEXPORT(cmfe_DataProjection_ResultProjectionXiSetObj) !Argument variables TYPE(cmfe_DataProjectionType), INTENT(IN) :: dataProjection !