From 9b4d9c087afdcc2b7330095bd2a1602bd799f64d Mon Sep 17 00:00:00 2001 From: Prasad Date: Sun, 3 Jun 2018 20:32:36 +1200 Subject: [PATCH 01/18] Fixes for dataprojections Distance vector sign was updated in previous commits by setting: distanceVector(1:numberOfCoordinates)=interpolatedPoint%values(:,NO_PART_DERIV)-dataPointLocation and by multiplying gradient and Hessian values by -1 compared with the original implementation. However, with the updated code, the gradient values are the same as what the old code produced, however the Hessian isn't. Fix was to revert to the original implementation for the gradient and Hessian and setting the distance vector as: distanceVector(1:numberOfCoordinates)=dataPointLocation-interpolatedPoint%values(:,NO_PART_DERIV) This has been done for all cases of the projection (line, face, element). --- src/data_projection_routines.f90 | 87 ++++++++++++++++++-------------- 1 file changed, 49 insertions(+), 38 deletions(-) diff --git a/src/data_projection_routines.f90 b/src/data_projection_routines.f90 index c53838dd..8792f40f 100644 --- a/src/data_projection_routines.f90 +++ b/src/data_projection_routines.f90 @@ -1184,6 +1184,16 @@ SUBROUTINE DataProjection_ErrorsCalculate(dataProjection,err,error,*) ELSE rmsError=SQRT(rmsError) ENDIF + IF(numberOfDataPoints==1) THEN + !Set maxDataPoint and maxError equal to minDataPoint and minError, respectively. + IF(maxDataPoint==0) THEN + maxDataPoint=minDataPoint + maxError=minError + ELSE + minDataPoint=maxDataPoint + minError=maxError + ENDIF + ENDIF dataProjection%rmsError=rmsError dataProjection%maximumError=maxError dataProjection%maximumErrorDataPoint=maxDataPoint @@ -1567,6 +1577,7 @@ SUBROUTINE DataProjection_DataPointsProjectionEvaluate(dataProjection,projection CALL FlagError(localError,err,error,*999) END SELECT ENDIF + IF(numberOfCandidates>dataProjection%maxNumberOfCandidates) dataProjection%maxNumberOfCandidates=numberOfCandidates !############################################################################################################## !find the clostest elements/faces/lines for each point in the current computational node base on starting xi !the clostest elements/faces/lines are required to shrink down on the list of possible projection candiates @@ -2260,7 +2271,7 @@ SUBROUTINE DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPo & 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)=interpolatedPoint%values(:,NO_PART_DERIV)-dataPointLocation + 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] @@ -2272,10 +2283,10 @@ SUBROUTINE DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPo bound=0 ENDIF !functionGradient - functionGradient=2.0_DP* & + functionGradient=-2.0_DP* & & (DOT_PRODUCT(distanceVector(1:numberOfCoordinates),interpolatedPoint%values(:,FIRST_PART_DERIV))) !functionHessian - functionHessian=2.0_DP* & + 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 @@ -2304,7 +2315,7 @@ SUBROUTINE DataProjection_NewtonElementsEvaluate_1(dataProjection,interpolatedPo newXi(1)=1.0_DP ENDIF CALL Field_InterpolateXi(SECOND_PART_DERIV,newXi,interpolatedPoint,err,error,*999,FIELD_GEOMETRIC_COMPONENTS_TYPE) - distanceVector=interpolatedPoint%values(:,NO_PART_DERIV)-dataPointLocation + 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) Date: Sun, 3 Jun 2018 21:17:06 +1200 Subject: [PATCH 02/18] Adding EQUATIONS_SET_FIBRE_FITTING_EQUATION_TYPE --- src/equations_set_constants.f90 | 1 + src/opencmiss_iron.f90 | 2 ++ 2 files changed, 3 insertions(+) diff --git a/src/equations_set_constants.f90 b/src/equations_set_constants.f90 index d250200d..40b3a983 100644 --- a/src/equations_set_constants.f90 +++ b/src/equations_set_constants.f90 @@ -103,6 +103,7 @@ MODULE EquationsSetConstants !Fitting class INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DATA_FITTING_EQUATION_TYPE=1 INTEGER(INTG), PARAMETER :: EQUATIONS_SET_GAUSS_FITTING_EQUATION_TYPE=2 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_FIBRE_FITTING_EQUATION_TYPE=3 !Multi physics class INTEGER(INTG), PARAMETER :: EQUATIONS_SET_FINITE_ELASTICITY_DARCY_TYPE=1 diff --git a/src/opencmiss_iron.f90 b/src/opencmiss_iron.f90 index ed0b281f..a3606760 100644 --- a/src/opencmiss_iron.f90 +++ b/src/opencmiss_iron.f90 @@ -2477,6 +2477,7 @@ MODULE OpenCMISS_Iron INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_LINEAR_ELASTIC_MODAL_TYPE = EQUATIONS_SET_LINEAR_ELASTIC_MODAL_TYPE ! Date: Sun, 3 Jun 2018 21:17:06 +1200 Subject: [PATCH 03/18] Revert "Adding EQUATIONS_SET_FIBRE_FITTING_EQUATION_TYPE" This reverts commit 3c7f4acdfac45fd4adfdf191c6b2e7fe9de46604. --- src/equations_set_constants.f90 | 1 - src/opencmiss_iron.f90 | 2 -- 2 files changed, 3 deletions(-) diff --git a/src/equations_set_constants.f90 b/src/equations_set_constants.f90 index 40b3a983..d250200d 100644 --- a/src/equations_set_constants.f90 +++ b/src/equations_set_constants.f90 @@ -103,7 +103,6 @@ MODULE EquationsSetConstants !Fitting class INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DATA_FITTING_EQUATION_TYPE=1 INTEGER(INTG), PARAMETER :: EQUATIONS_SET_GAUSS_FITTING_EQUATION_TYPE=2 - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_FIBRE_FITTING_EQUATION_TYPE=3 !Multi physics class INTEGER(INTG), PARAMETER :: EQUATIONS_SET_FINITE_ELASTICITY_DARCY_TYPE=1 diff --git a/src/opencmiss_iron.f90 b/src/opencmiss_iron.f90 index a3606760..ed0b281f 100644 --- a/src/opencmiss_iron.f90 +++ b/src/opencmiss_iron.f90 @@ -2477,7 +2477,6 @@ MODULE OpenCMISS_Iron INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_LINEAR_ELASTIC_MODAL_TYPE = EQUATIONS_SET_LINEAR_ELASTIC_MODAL_TYPE ! Date: Mon, 4 Jun 2018 14:25:21 +1200 Subject: [PATCH 04/18] Adding EQUATIONS_SET_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE --- src/equations_set_constants.f90 | 19 ++++++++++--------- src/opencmiss_iron.f90 | 5 ++++- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/equations_set_constants.f90 b/src/equations_set_constants.f90 index d250200d..15899b09 100644 --- a/src/equations_set_constants.f90 +++ b/src/equations_set_constants.f90 @@ -317,15 +317,16 @@ MODULE EquationsSetConstants !Fitting class !Data fitting INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE=1 - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_GENERALISED_DATA_FITTING_SUBTYPE=2 - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_MAT_PROPERTIES_DATA_FITTING_SUBTYPE=3 - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_MAT_PROPERTIES_INRIA_MODEL_DATA_FITTING_SUBTYPE=4 - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_VECTOR_DATA_FITTING_SUBTYPE=5 - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DIVFREE_VECTOR_DATA_FITTING_SUBTYPE=6 - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_VECTOR_DATA_PRE_FITTING_SUBTYPE=7 - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DIVFREE_VECTOR_DATA_PRE_FITTING_SUBTYPE=8 - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DATA_POINT_VECTOR_STATIC_FITTING_SUBTYPE=9 - INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DATA_POINT_VECTOR_QUASISTATIC_FITTING_SUBTYPE=10 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE=2 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_GENERALISED_DATA_FITTING_SUBTYPE=3 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_MAT_PROPERTIES_DATA_FITTING_SUBTYPE=4 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_MAT_PROPERTIES_INRIA_MODEL_DATA_FITTING_SUBTYPE=5 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_VECTOR_DATA_FITTING_SUBTYPE=6 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DIVFREE_VECTOR_DATA_FITTING_SUBTYPE=7 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_VECTOR_DATA_PRE_FITTING_SUBTYPE=8 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DIVFREE_VECTOR_DATA_PRE_FITTING_SUBTYPE=9 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DATA_POINT_VECTOR_STATIC_FITTING_SUBTYPE=10 + INTEGER(INTG), PARAMETER :: EQUATIONS_SET_DATA_POINT_VECTOR_QUASISTATIC_FITTING_SUBTYPE=11 !Gauss fitting INTEGER(INTG), PARAMETER :: EQUATIONS_SET_GAUSS_POINT_FITTING_SUBTYPE=1 !Fitting smoothing parameters diff --git a/src/opencmiss_iron.f90 b/src/opencmiss_iron.f90 index ed0b281f..13ae51c8 100644 --- a/src/opencmiss_iron.f90 +++ b/src/opencmiss_iron.f90 @@ -2750,6 +2750,8 @@ MODULE OpenCMISS_Iron INTEGER(INTG), PARAMETER :: CMFE_EQUATIONS_SET_SECOND_BIDOMAIN_SUBTYPE = EQUATIONS_SET_SECOND_BIDOMAIN_SUBTYPE ! Date: Mon, 4 Jun 2018 14:27:12 +1200 Subject: [PATCH 05/18] Adding EquationsSet setup for EQUATIONS_SET_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE Based this on those used for Finite Elasticity. --- src/fitting_routines.f90 | 567 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 565 insertions(+), 2 deletions(-) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 911cd9be..ca6d3807 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -1950,7 +1950,7 @@ SUBROUTINE Fitting_EquationsSetSetup(equationsSet,equationsSetSetup,err,error,*) SELECT CASE(equationsSet%specification(2)) CASE(EQUATIONS_SET_DATA_FITTING_EQUATION_TYPE) SELECT CASE(equationsSet%specification(3)) - CASE(EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE) + CASE(EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE,EQUATIONS_SET_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE) CALL Fitting_EquationsSetDataSetup(equationsSet,equationsSetSetup,err,error,*999) CASE(EQUATIONS_SET_VECTOR_DATA_FITTING_SUBTYPE,EQUATIONS_SET_VECTOR_DATA_PRE_FITTING_SUBTYPE) CALL FITTING_EQUATIONS_SET_VECTORDATA_SETUP(equationsSet,equationsSetSetup,err,error,*999) @@ -2022,7 +2022,8 @@ SUBROUTINE Fitting_EquationsSetSolutionMethodSet(equationsSet,solutionMethod,err SELECT CASE(equationsSet%specification(2)) CASE(EQUATIONS_SET_DATA_FITTING_EQUATION_TYPE) SELECT CASE(equationsSet%specification(3)) - CASE(EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE,EQUATIONS_SET_GENERALISED_DATA_FITTING_SUBTYPE, & + CASE(EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE,EQUATIONS_SET_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE, & + & EQUATIONS_SET_GENERALISED_DATA_FITTING_SUBTYPE, & & EQUATIONS_SET_MAT_PROPERTIES_DATA_FITTING_SUBTYPE,EQUATIONS_SET_MAT_PROPERTIES_INRIA_MODEL_DATA_FITTING_SUBTYPE, & & EQUATIONS_SET_VECTOR_DATA_FITTING_SUBTYPE,EQUATIONS_SET_VECTOR_DATA_PRE_FITTING_SUBTYPE, & & EQUATIONS_SET_DIVFREE_VECTOR_DATA_FITTING_SUBTYPE,EQUATIONS_SET_DIVFREE_VECTOR_DATA_PRE_FITTING_SUBTYPE, & @@ -2119,6 +2120,7 @@ SUBROUTINE Fitting_EquationsSetSpecificationSet(equationsSet,specification,err,e equationsSetSubtype=specification(3) SELECT CASE(equationsSetSubtype) CASE(EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE, & + & EQUATIONS_SET_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE, & & EQUATIONS_SET_VECTOR_DATA_FITTING_SUBTYPE, & & EQUATIONS_SET_DIVFREE_VECTOR_DATA_FITTING_SUBTYPE, & & EQUATIONS_SET_VECTOR_DATA_PRE_FITTING_SUBTYPE, & @@ -3357,6 +3359,567 @@ SUBROUTINE Fitting_EquationsSetDataSetup(equationsSet,equationsSetSetup,err,erro CALL FlagError(localError,err,error,*999) END SELECT + CASE DEFAULT + !----------------------------------------------------------------- + ! c a s e d e f a u l t + !----------------------------------------------------------------- + localError="The setup type of "//TRIM(NumberToVString(equationsSetSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a standard Galerkin projection." + CALL FlagError(localError,err,error,*999) + END SELECT + ELSEIF(equationsSet%specification(3)==EQUATIONS_SET_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE) THEN + SELECT CASE(equationsSetSetup%SETUP_TYPE) + + CASE(EQUATIONS_SET_SETUP_INITIAL_TYPE) + !----------------------------------------------------------------- + ! s o l u t i o n m e t h o d + !----------------------------------------------------------------- + SELECT CASE(equationsSetSetup%ACTION_TYPE) + CASE(EQUATIONS_SET_SETUP_START_ACTION) + CALL Fitting_EquationsSetSolutionMethodSet(equationsSet,EQUATIONS_SET_FEM_SOLUTION_METHOD,err,error,*999) + CALL EquationsSet_LabelSet(equationsSet,"Diffusion tensor fibre fitting equations set",err,error,*999) + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) + !Do nothing + CASE DEFAULT + localError="The action type of "//TRIM(NumberToVString(equationsSetSetup%ACTION_TYPE,"*",err,error))// & + & " for a setup type of "//TRIM(NumberToVString(equationsSetSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a standard Galerkin projection." + CALL FlagError(localError,err,error,*999) + END SELECT + + CASE(EQUATIONS_SET_SETUP_GEOMETRY_TYPE) + !----------------------------------------------------------------- + ! g e o m e t r y f i e l d + !----------------------------------------------------------------- + !Do nothing + + CASE(EQUATIONS_SET_SETUP_DEPENDENT_TYPE) + !----------------------------------------------------------------- + ! d e p e n d e n t f i e l d + !----------------------------------------------------------------- + SELECT CASE(equationsSetSetup%ACTION_TYPE) + CASE(EQUATIONS_SET_SETUP_START_ACTION) + IF(equationsSet%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN + !Create the auto created dependent field + CALL Field_CreateStart(equationsSetSetup%FIELD_USER_NUMBER,equationsSet%REGION,equationsSet%DEPENDENT% & + & DEPENDENT_FIELD,err,error,*999) + CALL Field_LabelSet(equationsSet%DEPENDENT%DEPENDENT_FIELD,"Dependent Field",err,error,*999) + CALL Field_TypeSetAndLock(equationsSet%DEPENDENT%DEPENDENT_FIELD,FIELD_GENERAL_TYPE,err,error,*999) + CALL Field_DependentTypeSetAndLock(equationsSet%DEPENDENT%DEPENDENT_FIELD,FIELD_DEPENDENT_TYPE,err,error,*999) + CALL Field_MeshDecompositionGet(equationsSet%GEOMETRY%GEOMETRIC_FIELD,geometricDecomposition,err,error,*999) + CALL Field_MeshDecompositionSetAndLock(equationsSet%DEPENDENT%DEPENDENT_FIELD,geometricDecomposition, & + & err,error,*999) + CALL Field_GeometricFieldSetAndLock(equationsSet%DEPENDENT%DEPENDENT_FIELD,equationsSet%GEOMETRY% & + & GEOMETRIC_FIELD,err,error,*999) + CALL Field_NumberOfVariablesSetAndLock(equationsSet%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999) + CALL Field_VariableTypesSetAndLock(equationsSet%DEPENDENT%DEPENDENT_FIELD,[FIELD_U_VARIABLE_TYPE, & + & FIELD_DELUDELN_VARIABLE_TYPE],err,error,*999) + CALL Field_VariableLabelSet(equationsSet%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE,"U",err,error,*999) + CALL Field_VariableLabelSet(equationsSet%DEPENDENT%DEPENDENT_FIELD,FIELD_DELUDELN_VARIABLE_TYPE,"del U/del n", & + & err,error,*999) + CALL Field_DimensionSetAndLock(equationsSet%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) + CALL Field_DimensionSetAndLock(equationsSet%DEPENDENT%DEPENDENT_FIELD,FIELD_DELUDELN_VARIABLE_TYPE, & + & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) + CALL Field_DataTypeSetAndLock(equationsSet%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_DP_TYPE,err,error,*999) + CALL Field_DataTypeSetAndLock(equationsSet%DEPENDENT%DEPENDENT_FIELD,FIELD_DELUDELN_VARIABLE_TYPE, & + & FIELD_DP_TYPE,err,error,*999) + !Set the number of components. + !If the independent field has been defined use that number of components + IF(ASSOCIATED(equationsSet%INDEPENDENT)) THEN + IF(ASSOCIATED(equationsSet%INDEPENDENT%INDEPENDENT_FIELD)) THEN + CALL Field_NumberOfComponentsGet(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & numberOfComponents,err,error,*999) + ELSE + numberOfComponents=1 + ENDIF + ELSE + numberOfComponents=1 + ENDIF + CALL Field_NumberOfComponentsSet(equationsSet%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & numberOfComponents,err,error,*999) + CALL Field_NumberOfComponentsSet(equationsSet%dependent%DEPENDENT_FIELD,FIELD_DELUDELN_VARIABLE_TYPE, & + & numberOfComponents,err,error,*999) + !Default to the geometric interpolation setup + CALL Field_ComponentMeshComponentGet(equationsSet%geometry%GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE,1, & + & geometricMeshComponent,err,error,*999) + DO componentIdx=1,numberOfComponents + CALL Field_ComponentMeshComponentSet(equationsSet%dependent%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & componentIdx,geometricMeshComponent,err,error,*999) + CALL Field_ComponentMeshComponentSet(equationsSet%dependent%DEPENDENT_FIELD,FIELD_DELUDELN_VARIABLE_TYPE, & + & componentIdx,geometricMeshComponent,err,error,*999) + ENDDO !componentIdx + SELECT CASE(equationsSet%SOLUTION_METHOD) + CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD) + DO componentIdx=1,numberOfComponents + CALL Field_ComponentInterpolationSetAndLock(equationsSet%dependent%DEPENDENT_FIELD, & + & FIELD_U_VARIABLE_TYPE,componentIdx,FIELD_NODE_BASED_INTERPOLATION,err,error,*999) + CALL Field_ComponentInterpolationSetAndLock(equationsSet%dependent%DEPENDENT_FIELD, & + & FIELD_DELUDELN_VARIABLE_TYPE,componentIdx,FIELD_NODE_BASED_INTERPOLATION,err,error,*999) + ENDDO !componentIdx + !Default the scaling to the geometric field scaling + CALL Field_ScalingTypeGet(equationsSet%geometry%GEOMETRIC_FIELD,geometricScalingType,err,error,*999) + CALL Field_ScalingTypeSet(equationsSet%dependent%DEPENDENT_FIELD,geometricScalingType,err,error,*999) + CASE(EQUATIONS_SET_BEM_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_FD_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_FV_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_GFEM_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_GFV_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE DEFAULT + localError="The solution method of "//TRIM(NumberToVString(equationsSet%SOLUTION_METHOD,"*",err,error))// & + & " is invalid." + CALL FlagError(localError,err,error,*999) + END SELECT + ELSE + !Check the user specified field + CALL Field_TypeCheck(equationsSetSetup%field,FIELD_GENERAL_TYPE,err,error,*999) + CALL Field_DependentTypeCheck(equationsSetSetup%field,FIELD_DEPENDENT_TYPE,err,error,*999) + CALL Field_NumberOfVariablesCheck(equationsSetSetup%field,2,err,error,*999) + CALL Field_VariableTypesCheck(equationsSetSetup%field,[FIELD_U_VARIABLE_TYPE,FIELD_DELUDELN_VARIABLE_TYPE], & + & err,error,*999) + CALL Field_DimensionCheck(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,FIELD_SCALAR_DIMENSION_TYPE,err,error,*999) + CALL Field_DimensionCheck(equationsSetSetup%field,FIELD_DELUDELN_VARIABLE_TYPE,FIELD_SCALAR_DIMENSION_TYPE, & + & err,error,*999) + CALL Field_DataTypeCheck(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,FIELD_DP_TYPE,err,error,*999) + CALL Field_DataTypeCheck(equationsSetSetup%field,FIELD_DELUDELN_VARIABLE_TYPE,FIELD_DP_TYPE,err,error,*999) + CALL Field_NumberOfComponentsGet(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,numberOfComponents,err,error,*999) + !If the independent field has been defined check the number of components is the same + IF(ASSOCIATED(equationsSet%INDEPENDENT)) THEN + IF(ASSOCIATED(equationsSet%INDEPENDENT%INDEPENDENT_FIELD)) THEN + CALL Field_NumberOfComponentsGet(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & numberOfIndependentComponents,err,error,*999) + IF(numberOfComponents /= numberOfIndependentComponents) THEN + localError="The number of components for the specified dependent field of "// & + & TRIM(NumberToVString(numberOfComponents,"*",err,error))// & + & " does not match the number of components for the independent field of "// & + & TRIM(NumberToVString(numberOfIndependentComponents,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + ENDIF + ENDIF + CALL Field_NumberOfComponentsGet(equationsSetSetup%field,FIELD_DELUDELN_VARIABLE_TYPE,numberOfComponents2, & + & err,error,*999) + IF(numberOfComponents2/=numberOfComponents) THEN + localError="The number of components in the independent field for variable type "// & + & TRIM(NumberToVstring(FIELD_DELUDELN_VARIABLE_TYPE,"*",err,error))//" of "// & + & TRIM(NumberToVstring(numberOfComponents2,"*",err,error))// & + & " does not match the number of components for variable type "// & + & TRIM(NumberToVstring(FIELD_U_VARIABLE_TYPE,"*",err,error))//" of "// & + & TRIM(NumberToVstring(numberOfComponents2,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + SELECT CASE(equationsSet%SOLUTION_METHOD) + CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD) + DO componentIdx=1,numberOfComponents + CALL Field_ComponentInterpolationCheck(equationsSetSetup%FIELD,FIELD_U_VARIABLE_TYPE,componentIdx, & + & FIELD_NODE_BASED_INTERPOLATION,err,error,*999) + CALL Field_ComponentInterpolationCheck(equationsSetSetup%FIELD,FIELD_DELUDELN_VARIABLE_TYPE,componentIdx, & + & FIELD_NODE_BASED_INTERPOLATION,err,error,*999) + ENDDO !numberOfComponents + CASE(EQUATIONS_SET_BEM_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_FD_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_FV_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_GFEM_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_GFV_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE DEFAULT + localError="The solution method of "//TRIM(NumberToVString(equationsSet%SOLUTION_METHOD,"*",err,error))// & + & " is invalid." + CALL FlagError(localError,err,error,*999) + END SELECT + ENDIF + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) + IF(equationsSet%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN + !Check that we have the same number of components as the independent field + IF(ASSOCIATED(equationsSet%dependent%DEPENDENT_FIELD)) THEN + CALL Field_NumberOfComponentsGet(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,numberOfComponents, & + & err,error,*999) + CALL Field_NumberOfComponentsGet(equationsSetSetup%field,FIELD_DELUDELN_VARIABLE_TYPE,numberOfComponents2, & + & err,error,*999) + IF(ASSOCIATED(equationsSet%INDEPENDENT)) THEN + IF(ASSOCIATED(equationsSet%INDEPENDENT%INDEPENDENT_FIELD)) THEN + CALL Field_NumberOfComponentsGet(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & numberOfIndependentComponents,err,error,*999) + IF(numberOfComponents /= numberOfIndependentComponents) THEN + localError="The number of components for the specified dependent field of "// & + & TRIM(NumberToVString(numberOfComponents,"*",err,error))// & + & " does not match the number of components for the independentt field of "// & + & TRIM(NumberToVString(numberOfIndependentComponents,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + ENDIF + ENDIF + IF(numberOfComponents2/=numberOfComponents) THEN + localError="The number of components in the dependent field for variable type "// & + & TRIM(NumberToVstring(FIELD_DELUDELN_VARIABLE_TYPE,"*",err,error))//" of "// & + & TRIM(NumberToVstring(numberOfComponents2,"*",err,error))// & + & " does not match the number of components for variable type "// & + & TRIM(NumberToVstring(FIELD_U_VARIABLE_TYPE,"*",err,error))//" of "// & + & TRIM(NumberToVstring(numberOfComponents2,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + !Finish creating the field + CALL Field_CreateFinish(equationsSet%DEPENDENT%DEPENDENT_FIELD,err,error,*999) + ELSE + CALL FlagError("The equations set dependent field is not associated.",err,error,*999) + ENDIF + ENDIF + CASE DEFAULT + localError="The action type of "//TRIM(NumberToVString(equationsSetSetup%ACTION_TYPE,"*",err,error))// & + & " for a setup type of "//TRIM(NumberToVString(equationsSetSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a standard Galerkin projection" + CALL FlagError(localError,err,error,*999) + END SELECT + + CASE(EQUATIONS_SET_SETUP_MATERIALS_TYPE) + !----------------------------------------------------------------- + ! m a t e r i a l f i e l d + !----------------------------------------------------------------- + SELECT CASE(equationsSet%specification(4)) + CASE(EQUATIONS_SET_FITTING_NO_SMOOTHING) + !Do nothing + CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_SMOOTHING,EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) + SELECT CASE(equationsSetSetup%ACTION_TYPE) + CASE(EQUATIONS_SET_SETUP_START_ACTION) + equationsMaterials=>equationsSet%materials + IF(ASSOCIATED(equationsMaterials)) THEN + IF(equationsMaterials%MATERIALS_FIELD_AUTO_CREATED) THEN + !Create the auto created materials field + CALL Field_CreateStart(equationsSetSetup%FIELD_USER_NUMBER,equationsSet%region,equationsSet% & + & materials%MATERIALS_FIELD,err,error,*999) + CALL Field_TypeSetAndLock(equationsMaterials%MATERIALS_FIELD,FIELD_MATERIAL_TYPE,err,error,*999) + CALL Field_DependentTypeSetAndLock(equationsMaterials%MATERIALS_FIELD,FIELD_INDEPENDENT_TYPE, & + & err,error,*999) + CALL Field_MeshDecompositionGet(equationsSet%geometry%GEOMETRIC_FIELD,geometricDecomposition, & + & err,error,*999) + !apply decomposition rule found on new created field + CALL Field_MeshDecompositionSetAndLock(equationsSet%materials%MATERIALS_FIELD,geometricDecomposition, & + & err,error,*999) + !point new field to geometric field + CALL Field_GeometricFieldSetAndLock(equationsMaterials%MATERIALS_FIELD,equationsSet%geometry% & + & GEOMETRIC_FIELD,err,error,*999) + CALL Field_NumberOfVariablesSet(equationsMaterials%MATERIALS_FIELD,1,err,error,*999) + CALL Field_VariableTypesSetAndLock(equationsMaterials%MATERIALS_FIELD,[FIELD_U_VARIABLE_TYPE], & + & err,error,*999) + CALL Field_DimensionSetAndLock(equationsMaterials%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) + CALL Field_DataTypeSetAndLock(equationsMaterials%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_DP_TYPE,err,error,*999) + !Sobolev smoothing material parameters- tau and kappa + CALL Field_NumberOfComponentsSetAndLock(equationsMaterials%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & 2,err,error,*999) + CALL Field_ComponentMeshComponentGet(equationsSet%geometry%GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE, & + & 1,geometricMeshComponent,err,error,*999) + CALL Field_ComponentMeshComponentSet(equationsMaterials%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & 1,geometricMeshComponent,err,error,*999) + CALL Field_ComponentMeshComponentSet(equationsMaterials%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & 2,geometricMeshComponent,err,error,*999) + CALL Field_ComponentInterpolationSet(equationsMaterials%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & 1,FIELD_CONSTANT_INTERPOLATION,err,error,*999) + CALL Field_ComponentInterpolationSet(equationsMaterials%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & 2,FIELD_CONSTANT_INTERPOLATION,err,error,*999) + !Default the field scaling to that of the geometric field + CALL Field_ScalingTypeGet(equationsSet%geometry%GEOMETRIC_FIELD,geometricScalingType,err,error,*999) + CALL Field_ScalingTypeSet(equationsMaterials%MATERIALS_FIELD,geometricScalingType,err,error,*999) + ELSE + !Check the user specified field + CALL Field_TypeCheck(equationsSetSetup%field,FIELD_MATERIAL_TYPE,err,error,*999) + CALL Field_DependentTypeCheck(equationsSetSetup%field,FIELD_INDEPENDENT_TYPE,err,error,*999) + CALL Field_NumberOfVariablesCheck(equationsSetSetup%field,1,err,error,*999) + CALL Field_VariableTypesCheck(equationsSetSetup%field,[FIELD_U_VARIABLE_TYPE],err,error,*999) + CALL Field_DimensionCheck(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,FIELD_VECTOR_DIMENSION_TYPE, & + & err,error,*999) + CALL Field_DataTypeCheck(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,FIELD_DP_TYPE,err,error,*999) + CALL Field_NumberOfComponentsCheck(equationsSetSetup%FIELD,FIELD_U_VARIABLE_TYPE,2,err,error,*999) + ENDIF + ELSE + CALL FlagError("Equations set materials is not associated.",err,error,*999) + END IF + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) + equationsMaterials=>equationsSet%materials + IF(ASSOCIATED(equationsMaterials)) THEN + IF(equationsMaterials%MATERIALS_FIELD_AUTO_CREATED) THEN + !Finish creating the materials field + CALL Field_CreateFinish(equationsMaterials%MATERIALS_FIELD,err,error,*999) + !Set the default values for the materials field + CALL Field_ComponentValuesInitialise(equationsMaterials%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_VALUES_SET_TYPE,1,0.0_DP,err,error,*999) + CALL Field_ComponentValuesInitialise(equationsMaterials%MATERIALS_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_VALUES_SET_TYPE,2,0.0_DP,err,error,*999) + ENDIF + ELSE + CALL FlagError("Equations set materials is not associated.",err,error,*999) + ENDIF + CASE DEFAULT + localError="The action type of "//TRIM(NumberToVString(equationsSetSetup%ACTION_TYPE,"*",err,error))// & + & " for a setup type of "//TRIM(NumberToVString(equationsSetSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for an update-materials Galerkin projection." + CALL FlagError(localError,err,error,*999) + END SELECT + CASE(EQUATIONS_SET_FITTING_STRAIN_ENERGY_SMOOTHING) + CALL FlagError("Not implemented.",err,error,*999) + CASE DEFAULT + localError="The fitting smoothing type of "//TRIM(NumberToVString(equationsSet%specification(4),"*",err,error))// & + & " is invalid." + CALL FlagError(localError,err,error,*999) + END SELECT + + CASE(EQUATIONS_SET_SETUP_INDEPENDENT_TYPE) + !----------------------------------------------------------------- + ! I n d e p e n d e n t t y p e + ! + ! (this field holds the data point based field of vectors to map to the dependent field) + !----------------------------------------------------------------- + SELECT CASE(equationsSetSetup%ACTION_TYPE) + !Set start action + CASE(EQUATIONS_SET_SETUP_START_ACTION) + IF(equationsSet%INDEPENDENT%INDEPENDENT_FIELD_AUTO_CREATED) THEN + !Create the auto created independent field + !start field creation with name 'INDEPENDENT_FIELD' + CALL Field_CreateStart(equationsSetSetup%FIELD_USER_NUMBER,equationsSet%region, & + & equationsSet%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999) + !start creation of a new field + CALL Field_TypeSetAndLock(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,FIELD_GENERAL_TYPE,err,error,*999) + !label the field + CALL Field_LabelSetAndLock(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,"Independent Field",err,error, & + & *999) + !define new created field to be independent + CALL Field_DependentTypeSetAndLock(equationsSet%independent%INDEPENDENT_FIELD, & + & FIELD_INDEPENDENT_TYPE,err,error,*999) + !look for decomposition rule already defined + CALL Field_MeshDecompositionGet(equationsSet%geometry%GEOMETRIC_FIELD,geometricDecomposition,err,error,*999) + !apply decomposition rule found on new created field + CALL Field_MeshDecompositionSetAndLock(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,geometricDecomposition, & + & err,error,*999) + !point new field to geometric field + CALL Field_GeometricFieldSetAndLock(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,equationsSet% & + & geometry%GEOMETRIC_FIELD,err,error,*999) + !Create two variables: U for data and V for weights + CALL Field_NumberOfVariablesSetAndLock(equationsSet%INDEPENDENT%INDEPENDENT_FIELD, & + & 2,err,error,*999) + CALL Field_VariableTypesSetAndLock(equationsSet%INDEPENDENT%INDEPENDENT_FIELD, & + & [FIELD_U_VARIABLE_TYPE,FIELD_V_VARIABLE_TYPE],err,error,*999) + CALL Field_ComponentMeshComponentGet(equationsSet%geometry%GEOMETRIC_FIELD,FIELD_U_VARIABLE_TYPE, & + & 1,geometricMeshComponent,err,error,*999) + !If the dependent field has been created then use that number of components + IF(ASSOCIATED(equationsSet%DEPENDENT%DEPENDENT_FIELD)) THEN + CALL Field_NumberOfComponentsGet(equationsSet%DEPENDENT%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & numberOfComponents,err,error,*999) + ELSE + numberOfComponents=1 + ENDIF + ! U Variable: data points + CALL Field_DimensionSetAndLock(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) + CALL Field_DataTypeSetAndLock(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & FIELD_DP_TYPE,err,error,*999) + CALL Field_NumberOfComponentsSet(equationsSet%INDEPENDENT%INDEPENDENT_FIELD, & + & FIELD_U_VARIABLE_TYPE,numberOfComponents,err,error,*999) + !Default to the geometric interpolation setup + ! V Variable: data point weights + CALL Field_DimensionSetAndLock(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,FIELD_V_VARIABLE_TYPE, & + & FIELD_VECTOR_DIMENSION_TYPE,err,error,*999) + CALL Field_DataTypeSetAndLock(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,FIELD_V_VARIABLE_TYPE, & + & FIELD_DP_TYPE,err,error,*999) + CALL Field_NumberOfComponentsSet(equationsSet%INDEPENDENT%INDEPENDENT_FIELD, & + & FIELD_V_VARIABLE_TYPE,numberOfComponents,err,error,*999) + !Default to the geometric interpolation setup + DO componentIdx=1,numberOfComponents + CALL Field_ComponentMeshComponentSet(equationsSet%INDEPENDENT%INDEPENDENT_FIELD, & + & FIELD_U_VARIABLE_TYPE,componentIdx,geometricMeshComponent,err,error,*999) + CALL Field_ComponentMeshComponentSet(equationsSet%INDEPENDENT%INDEPENDENT_FIELD, & + & FIELD_V_VARIABLE_TYPE,componentIdx,geometricMeshComponent,err,error,*999) + ENDDO !componentIdx + SELECT CASE(equationsSet%SOLUTION_METHOD) + !Specify fem solution method + CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD) + DO componentIdx = 1,numberOfComponents + CALL Field_ComponentInterpolationSetAndLock(equationsSet%independent%INDEPENDENT_FIELD, & + & FIELD_U_VARIABLE_TYPE,componentIdx,FIELD_DATA_POINT_BASED_INTERPOLATION,err,error,*999) + CALL Field_ComponentInterpolationSetAndLock(equationsSet%independent%INDEPENDENT_FIELD, & + & FIELD_V_VARIABLE_TYPE,componentIdx,FIELD_DATA_POINT_BASED_INTERPOLATION,err,error,*999) + ENDDO !componentIdx + CALL Field_ScalingTypeGet(equationsSet%geometry%GEOMETRIC_FIELD,geometricScalingType,err,error,*999) + CALL Field_ScalingTypeSet(equationsSet%independent%INDEPENDENT_FIELD,geometricScalingType,err,error,*999) + CASE DEFAULT + localError="The solution method of " & + & //TRIM(NumberToVString(equationsSet%SOLUTION_METHOD,"*",err,error))// " is invalid." + CALL FlagError(localError,err,error,*999) + END SELECT + ELSE + !Check the user specified field + CALL Field_TypeCheck(equationsSetSetup%field,FIELD_GENERAL_TYPE,err,error,*999) + CALL Field_DependentTypeCheck(equationsSetSetup%field,FIELD_INDEPENDENT_TYPE,err,error,*999) + ! U (vector) variable + CALL Field_DimensionCheck(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,FIELD_VECTOR_DIMENSION_TYPE, & + & err,error,*999) + CALL Field_DataTypeCheck(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,FIELD_DP_TYPE,err,error,*999) + CALL Field_NumberOfComponentsGet(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE, & + & numberOfComponents,err,error,*999) + !If the dependent field has been defined use that number of components + IF(ASSOCIATED(equationsSet%dependent%DEPENDENT_FIELD)) THEN + CALL Field_NumberOfComponentsGet(equationsSet%dependent%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & numberOfDependentComponents,err,error,*999) + IF(numberOfComponents /= numberOfDependentComponents) THEN + localError="The number of components for the specified independent field of "// & + & TRIM(NumberToVString(numberOfComponents,"*",err,error))// & + & " does not match the number of components for the dependent field of "// & + & TRIM(NumberToVString(numberOfDependentComponents,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + ENDIF + ! V (weight) variable + CALL Field_DimensionCheck(equationsSetSetup%FIELD,FIELD_V_VARIABLE_TYPE,FIELD_VECTOR_DIMENSION_TYPE, & + & err,error,*999) + CALL Field_DataTypeCheck(equationsSetSetup%FIELD,FIELD_V_VARIABLE_TYPE,FIELD_DP_TYPE,err,error,*999) + CALL Field_NumberOfComponentsCheck(equationsSetSetup%field,FIELD_V_VARIABLE_TYPE,numberOfComponents,err,error,*999) + SELECT CASE(equationsSet%SOLUTION_METHOD) + CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD) + DO componentIdx=1,numberOfComponents + CALL Field_ComponentInterpolationCheck(equationsSetSetup%FIELD,FIELD_U_VARIABLE_TYPE,componentIdx, & + & FIELD_GAUSS_POINT_BASED_INTERPOLATION,err,error,*999) + CALL Field_ComponentInterpolationCheck(equationsSetSetup%FIELD,FIELD_V_VARIABLE_TYPE,componentIdx, & + & FIELD_GAUSS_POINT_BASED_INTERPOLATION,err,error,*999) + ENDDO !componentIdx + CASE DEFAULT + localError="The solution method of "//TRIM(NumberToVString(equationsSet%SOLUTION_METHOD, & + &"*",err,error))//" is invalid." + CALL FlagError(localError,err,error,*999) + END SELECT + ENDIF + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) + IF(equationsSet%independent%INDEPENDENT_FIELD_AUTO_CREATED) THEN + !Check that we have the same number of components as the dependent field + IF(ASSOCIATED(equationsSet%dependent%DEPENDENT_FIELD)) THEN + CALL Field_NumberOfComponentsGet(equationsSet%dependent%DEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE, & + & numberOfDependentComponents,err,error,*999) + CALL Field_NumberOfComponentsGet(equationsSetSetup%field,FIELD_U_VARIABLE_TYPE,numberOfComponents,err,error,*999) + IF(numberOfComponents /= numberOfDependentComponents) THEN + localError="The number of components for the specified independent field of "// & + & TRIM(NumberToVString(numberOfComponents,"*",err,error))// & + & " does not match the number of components for the dependentt field of "// & + & TRIM(NumberToVString(numberOfDependentComponents,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + CALL Field_NumberOfComponentsGet(equationsSetSetup%field,FIELD_V_VARIABLE_TYPE,numberOfComponents2,err,error,*999) + IF(numberOfComponents /= numberOfComponents2) THEN + localError="The number of components for the U variable of the specified independent field of "// & + & TRIM(NumberToVString(numberOfComponents,"*",err,error))// & + & " does not match the number of components for the V variable of the field of "// & + & TRIM(NumberToVString(numberOfComponents2,"*",err,error))//"." + CALL FlagError(localError,err,error,*999) + ENDIF + ENDIF + !Specify finish action + CALL Field_CreateFinish(equationsSet%INDEPENDENT%INDEPENDENT_FIELD,err,error,*999) + !Initialise the weights to 1.0 + DO componentIdx=1,numberOfComponents + CALL Field_ComponentValuesInitialise(equationsSet%independent%INDEPENDENT_FIELD,FIELD_V_VARIABLE_TYPE, & + & FIELD_VALUES_SET_TYPE,componentIdx,1.0_DP,err,error,*999) + ENDDO !componentIdx + ENDIF + CASE DEFAULT + localError="The action type of "//TRIM(NumberToVString(equationsSetSetup%ACTION_TYPE,"*",err,error))// & + & " for a setup type of "//TRIM(NumberToVString(equationsSetSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a fitting equations set." + CALL FlagError(localError,err,error,*999) + END SELECT + CASE(EQUATIONS_SET_SETUP_SOURCE_TYPE) + !----------------------------------------------------------------- + ! s o u r c e t y p e + !----------------------------------------------------------------- + SELECT CASE(equationsSetSetup%ACTION_TYPE) + CASE(EQUATIONS_SET_SETUP_START_ACTION) + !Do nothing + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) + !Do nothing + CASE DEFAULT + localError="The action type of "//TRIM(NumberToVString(equationsSetSetup%ACTION_TYPE,"*",err,error))// & + & " for a setup type of "//TRIM(NumberToVString(equationsSetSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a standard Galerkin projection." + CALL FlagError(localError,err,error,*999) + END SELECT + + CASE(EQUATIONS_SET_SETUP_EQUATIONS_TYPE) + !----------------------------------------------------------------- + ! e q u a t i o n s t y p e + !----------------------------------------------------------------- + SELECT CASE(equationsSetSetup%ACTION_TYPE) + CASE(EQUATIONS_SET_SETUP_START_ACTION) + IF(equationsSet%DEPENDENT%DEPENDENT_FINISHED) THEN + CALL Equations_CreateStart(equationsSet,equations,err,error,*999) + CALL Equations_LinearityTypeSet(equations,EQUATIONS_NONLINEAR,err,error,*999) + CALL Equations_TimeDependenceTypeSet(equations,EQUATIONS_STATIC,err,error,*999) + ELSE + CALL FlagError("Equations set dependent field has not been finished.",err,error,*999) + ENDIF + CASE(EQUATIONS_SET_SETUP_FINISH_ACTION) + SELECT CASE(equationsSet%SOLUTION_METHOD) + CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD) + !Finish the equations creation + CALL EquationsSet_EquationsGet(equationsSet,equations,err,error,*999) + CALL Equations_CreateFinish(equations,err,error,*999) + NULLIFY(vectorEquations) + CALL Equations_VectorEquationsGet(equations,vectorEquations,err,error,*999) + !Create the equations mapping. + CALL EquationsMapping_VectorCreateStart(vectorEquations,FIELD_DELUDELN_VARIABLE_TYPE,vectorMapping,err,error,*999) + CALL EquationsMapping_ResidualVariableTypesSet(vectorMapping,[FIELD_U_VARIABLE_TYPE],err,error,*999) + CALL EquationsMapping_LinearMatricesNumberSet(vectorMapping,0,err,error,*999) + CALL EquationsMapping_RHSVariableTypeSet(vectorMapping,FIELD_DELUDELN_VARIABLE_TYPE,err,error,*999) + CALL EquationsMapping_VectorCreateFinish(vectorMapping,err,error,*999) + !Create the equations matrices + CALL EquationsMatrices_VectorCreateStart(vectorEquations,vectorMatrices,err,error,*999) + ! set structure and storage types + SELECT CASE(equations%sparsityType) + CASE(EQUATIONS_MATRICES_FULL_MATRICES) + CALL EquationsMatrices_NonlinearStorageTypeSet(vectorMatrices,MATRIX_BLOCK_STORAGE_TYPE, & + & err,error,*999) + CASE(EQUATIONS_MATRICES_SPARSE_MATRICES) + CALL EquationsMatrices_NonlinearStorageTypeSet(vectorMatrices,MATRIX_COMPRESSED_ROW_STORAGE_TYPE, & + & err,error,*999) + CALL EquationsMatrices_NonlinearStructureTypeSet(vectorMatrices,EQUATIONS_MATRIX_FEM_STRUCTURE, & + & err,error,*999) + CASE DEFAULT + localError="The equations matrices sparsity type of "// & + & TRIM(NumberToVString(equations%sparsityType,"*",err,error))//" is invalid." + CALL FlagError(localError,err,error,*999) + END SELECT + !Set Jacobian matrices calculation type to default finite difference. + CALL EquationsMatrices_JacobianTypesSet(vectorMatrices,[EQUATIONS_JACOBIAN_FINITE_DIFFERENCE_CALCULATED], & + & err,error,*999) + CALL EquationsMatrices_VectorCreateFinish(vectorMatrices,err,error,*999) + CASE(EQUATIONS_SET_BEM_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_FD_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_FV_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_GFEM_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE(EQUATIONS_SET_GFV_SOLUTION_METHOD) + CALL FlagError("Not implemented.",err,error,*999) + CASE DEFAULT + localError="The solution method of "//TRIM(NumberToVString(equationsSet%SOLUTION_METHOD,"*",err,error))// & + & " is invalid." + CALL FlagError(localError,err,error,*999) + END SELECT + CASE DEFAULT + localError="The action type of "//TRIM(NumberToVString(equationsSetSetup%ACTION_TYPE,"*",err,error))// & + & " for a setup type of "//TRIM(NumberToVString(equationsSetSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a standard Galerkin projection." + CALL FlagError(localError,err,error,*999) + END SELECT + CASE DEFAULT !----------------------------------------------------------------- ! c a s e d e f a u l t From b6e131fb5db7899fecb6bf49da2277c6089446a0 Mon Sep 17 00:00:00 2001 From: Prasad Date: Mon, 4 Jun 2018 14:41:26 +1200 Subject: [PATCH 06/18] Renaming existing Fitting_FiniteElementCalculate to LinearFitting_FiniteElementCalculate (linear) Fitting_FiniteElementCalculate was originally used for EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE --- src/equations_set_routines.f90 | 8 +++++++- src/fitting_routines.f90 | 12 ++++++------ 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/src/equations_set_routines.f90 b/src/equations_set_routines.f90 index c37c09c7..3bd27ba9 100644 --- a/src/equations_set_routines.f90 +++ b/src/equations_set_routines.f90 @@ -2598,7 +2598,13 @@ SUBROUTINE EquationsSet_FiniteElementCalculate(equationsSet,elementNumber,err,er CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) CALL CLASSICAL_FIELD_FINITE_ELEMENT_CALCULATE(equationsSet,elementNumber,err,error,*999) CASE(EQUATIONS_SET_FITTING_CLASS) - CALL Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,*999) + IF(equationsSet%specification(3)==EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE) THEN + CALL LinearFitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,*999) + ELSE + localError="The equations set subtype of "//TRIM(NumberToVString(equationsSet%specification(3),"*",err,error))// & + & " does not equal a standard Galerkin projection subtype." + CALL FlagError(localError,err,error,*999) + ENDIF CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) IF(SIZE(equationsSet%specification,1)<2) & & CALL FlagError("Equations set specification must have at least two entries for a bioelectrics equation class.", & diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index ca6d3807..32245388 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -97,7 +97,7 @@ MODULE FittingRoutines PUBLIC Fitting_ProblemSetup PUBLIC Fitting_ProblemSpecificationSet - PUBLIC Fitting_FiniteElementCalculate + PUBLIC LinearFitting_FiniteElementCalculate PUBLIC Fitting_PreSolve PUBLIC Fitting_PostSolve @@ -109,7 +109,7 @@ MODULE FittingRoutines ! !>Calculates the element stiffness matrices and RHS for a Galerkin projection finite element equations set. - SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,*) + SUBROUTINE LinearFitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,*) !Argument variables TYPE(EQUATIONS_SET_TYPE), POINTER :: equationsSet ! Date: Mon, 4 Jun 2018 15:20:04 +1200 Subject: [PATCH 07/18] Adding skeleton for NonlinearFitting_FiniteElementResidualEvaluate Setting FiniteElementJacobianEvaluate to not implemented and instead will be calculated with FD --- src/equations_set_routines.f90 | 12 ++++---- src/fitting_routines.f90 | 54 ++++++++++++++++++++++++++++++++++ 2 files changed, 59 insertions(+), 7 deletions(-) diff --git a/src/equations_set_routines.f90 b/src/equations_set_routines.f90 index 3bd27ba9..f99098f3 100644 --- a/src/equations_set_routines.f90 +++ b/src/equations_set_routines.f90 @@ -2598,13 +2598,7 @@ SUBROUTINE EquationsSet_FiniteElementCalculate(equationsSet,elementNumber,err,er CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) CALL CLASSICAL_FIELD_FINITE_ELEMENT_CALCULATE(equationsSet,elementNumber,err,error,*999) CASE(EQUATIONS_SET_FITTING_CLASS) - IF(equationsSet%specification(3)==EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE) THEN - CALL LinearFitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,*999) - ELSE - localError="The equations set subtype of "//TRIM(NumberToVString(equationsSet%specification(3),"*",err,error))// & - & " does not equal a standard Galerkin projection subtype." - CALL FlagError(localError,err,error,*999) - ENDIF + CALL LinearFitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,*999) CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) IF(SIZE(equationsSet%specification,1)<2) & & CALL FlagError("Equations set specification must have at least two entries for a bioelectrics equation class.", & @@ -2779,6 +2773,8 @@ SUBROUTINE EquationsSet_FiniteElementJacobianEvaluate(equationsSet,elementNumber CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) CALL ClassicalField_FiniteElementJacobianEvaluate(equationsSet,elementNumber,err,error,*999) + CASE(EQUATIONS_SET_FITTING_CLASS) + CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_MODAL_CLASS) @@ -3018,6 +3014,8 @@ SUBROUTINE EquationsSet_FiniteElementResidualEvaluate(equationsSet,elementNumber CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) CALL ClassicalField_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*999) + CASE(EQUATIONS_SET_FITTING_CLASS) + CALL NonlinearFitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*999) CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_MODAL_CLASS) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 32245388..0d344d13 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -98,6 +98,7 @@ MODULE FittingRoutines PUBLIC Fitting_ProblemSpecificationSet PUBLIC LinearFitting_FiniteElementCalculate + PUBLIC NonlinearFitting_FiniteElementResidualEvaluate PUBLIC Fitting_PreSolve PUBLIC Fitting_PostSolve @@ -1338,6 +1339,59 @@ END SUBROUTINE LinearFitting_FiniteElementCalculate !================================================================================================================================ ! + !>Evaluates the element residual and rhs vectors for the given element number for a fitting class finite element equation set. + SUBROUTINE NonlinearFitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*) + + !Argument variables + TYPE(EQUATIONS_SET_TYPE), POINTER :: equationsSet !equationsSet%EQUATIONS + IF(ASSOCIATED(EQUATIONS)) THEN + NULLIFY(vectorEquations) + CALL Equations_VectorEquationsGet(equations,vectorEquations,err,error,*999) + SELECT CASE(equationsSet%SPECIFICATION(3)) + CASE(EQUATIONS_SET_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE) + + CASE DEFAULT + localError="Equations set subtype "//TRIM(NumberToVString(equationsSet%SPECIFICATION(3),"*",err,error))// & + & " is not valid for a fitting equation set class." + CALL FlagError(localError,err,error,*999) + END SELECT + ELSE + CALL FlagError("Equations set equations is not associated.",err,error,*999) + ENDIF + ELSE + CALL FlagError("Equations set is not associated.",err,error,*999) + ENDIF + + EXITS("NonlinearFitting_FiniteElementResidualEvaluate") + RETURN +999 ERRORS("NonlinearFitting_FiniteElementResidualEvaluate",err,error) + EXITS("NonlinearFitting_FiniteElementResidualEvaluate") + RETURN 1 + + END SUBROUTINE NonlinearFitting_FiniteElementResidualEvaluate + + ! + !================================================================================================================================ + ! + !>Sets up the update-materials Galerkin projection. SUBROUTINE FITTING_EQUATIONS_SET_MAT_PROPERTIES_SETUP(equationsSet,equationsSetSetup,err,error,*) From bd50d784136da0c26fb149bdfc4c3908cc6c18f4 Mon Sep 17 00:00:00 2001 From: Prasad Date: Mon, 4 Jun 2018 15:36:39 +1200 Subject: [PATCH 08/18] Adding PROBLEM_STATIC_NONLINEAR_FITTING_SUBTYPE --- src/fitting_routines.f90 | 101 +++++++++++++++++++++++++++++++++++++- src/opencmiss_iron.f90 | 5 +- src/problem_constants.f90 | 19 +++---- 3 files changed, 114 insertions(+), 11 deletions(-) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 0d344d13..22786c57 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -5468,7 +5468,97 @@ SUBROUTINE Fitting_ProblemStaticSetup(problem,problemSetup,err,error,*) & " is invalid for a static fitting problem." CALL FlagError(localError,err,error,*999) END SELECT - CASE DEFAULT + CASE DEFAULT + localError="The setup type of "//TRIM(NumberToVString(problemSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a static fitting problem." + CALL FlagError(localError,err,error,*999) + END SELECT + ELSEIF(problem%specification(3)==PROBLEM_STATIC_NONLINEAR_FITTING_SUBTYPE) THEN + SELECT CASE(problemSetup%SETUP_TYPE) + CASE(PROBLEM_SETUP_INITIAL_TYPE) + SELECT CASE(problemSetup%ACTION_TYPE) + CASE(PROBLEM_SETUP_START_ACTION) + !Do nothing???? + CASE(PROBLEM_SETUP_FINISH_ACTION) + !Do nothing??? + CASE DEFAULT + localError="The action type of "//TRIM(NumberToVString(problemSetup%ACTION_TYPE,"*",err,error))// & + & " for a setup type of "//TRIM(NumberToVString(problemSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a static nonlinear fitting problem." + CALL FlagError(localError,err,error,*999) + END SELECT + CASE(PROBLEM_SETUP_CONTROL_TYPE) + SELECT CASE(problemSetup%ACTION_TYPE) + CASE(PROBLEM_SETUP_START_ACTION) + !Set up a simple control loop + CALL ControlLoop_CreateStart(problem,controlLoop,err,error,*999) + CASE(PROBLEM_SETUP_FINISH_ACTION) + !Finish the control loops + controlLoopRoot=>problem%CONTROL_LOOP + CALL ControlLoop_Get(controlLoopRoot,CONTROL_LOOP_NODE,controlLoop,err,error,*999) + CALL ControlLoop_CreateFinish(controlLoop,err,error,*999) + CASE DEFAULT + localError="The action type of "//TRIM(NumberToVString(problemSetup%ACTION_TYPE,"*",err,error))// & + & " for a setup type of "//TRIM(NumberToVString(problemSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a static fitting problem." + CALL FlagError(localError,err,error,*999) + END SELECT + CASE(PROBLEM_SETUP_SOLVERS_TYPE) + !Get the control loop + controlLoopRoot=>problem%CONTROL_LOOP + CALL ControlLoop_Get(controlLoopRoot,CONTROL_LOOP_NODE,controlLoop,err,error,*999) + SELECT CASE(problemSetup%ACTION_TYPE) + CASE(PROBLEM_SETUP_START_ACTION) + !Start the solvers creation + CALL Solvers_CreateStart(controlLoop,solvers,err,error,*999) + CALL Solvers_NumberSet(solvers,1,err,error,*999) + !Set the solver to be a nonlinear solver + CALL Solvers_SolverGet(solvers,1,solver,err,error,*999) + CALL Solver_TypeSet(solver,SOLVER_NONLINEAR_TYPE,err,error,*999) + !Set solver defaults + CALL Solver_LibraryTypeSet(solver,SOLVER_PETSC_LIBRARY,err,error,*999) + CASE(PROBLEM_SETUP_FINISH_ACTION) + !Get the solvers + CALL ControlLoop_SolversGet(controlLoop,solvers,err,error,*999) + !Finish the solvers creation + CALL Solvers_CreateFinish(solvers,err,error,*999) + CASE DEFAULT + localError="The action type of "//TRIM(NumberToVString(problemSetup%ACTION_TYPE,"*",err,error))// & + & " for a setup type of "//TRIM(NumberToVString(problemSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a static fitting problem." + CALL FlagError(localError,err,error,*999) + END SELECT + CASE(PROBLEM_SETUP_SOLVER_EQUATIONS_TYPE) + SELECT CASE(problemSetup%ACTION_TYPE) + CASE(PROBLEM_SETUP_START_ACTION) + !Get the control loop + controlLoopRoot=>problem%CONTROL_LOOP + CALL ControlLoop_Get(controlLoopRoot,CONTROL_LOOP_NODE,controlLoop,err,error,*999) + !Get the solver + CALL ControlLoop_SolversGet(controlLoop,solvers,err,error,*999) + CALL Solvers_SolverGet(solvers,1,solver,err,error,*999) + !Create the solver equations + CALL SolverEquations_CreateStart(solver,solverEquations,err,error,*999) + CALL SolverEquations_LinearityTypeSet(solverEquations,SOLVER_EQUATIONS_NONLINEAR,err,error,*999) + CALL SolverEquations_TimeDependenceTypeSet(solverEquations,SOLVER_EQUATIONS_STATIC,err,error,*999) + CALL SolverEquations_SparsityTypeSet(solverEquations,SOLVER_SPARSE_MATRICES,err,error,*999) + CASE(PROBLEM_SETUP_FINISH_ACTION) + !Get the control loop + controlLoopRoot=>problem%CONTROL_LOOP + CALL ControlLoop_Get(controlLoopRoot,CONTROL_LOOP_NODE,controlLoop,err,error,*999) + !Get the solver equations + CALL ControlLoop_SolversGet(controlLoop,solvers,err,error,*999) + CALL Solvers_SolverGet(solvers,1,solver,err,error,*999) + CALL Solver_SolverEquationsGet(solver,solverEquations,err,error,*999) + !Finish the solver equations creation + CALL SolverEquations_CreateFinish(solverEquations,err,error,*999) + CASE DEFAULT + localError="The action type of "//TRIM(NumberToVString(problemSetup%ACTION_TYPE,"*",err,error))// & + & " for a setup type of "//TRIM(NumberToVString(problemSetup%SETUP_TYPE,"*",err,error))// & + & " is invalid for a static fitting problem." + CALL FlagError(localError,err,error,*999) + END SELECT + CASE DEFAULT localError="The setup type of "//TRIM(NumberToVString(problemSetup%SETUP_TYPE,"*",err,error))// & & " is invalid for a static fitting problem." CALL FlagError(localError,err,error,*999) @@ -5663,6 +5753,7 @@ SUBROUTINE Fitting_ProblemSpecificationSet(problem,problemSpecification,err,erro CASE(PROBLEM_DATA_FITTING_TYPE) SELECT CASE(problemSubtype) CASE(PROBLEM_STATIC_FITTING_SUBTYPE, & + & PROBLEM_STATIC_NONLINEAR_FITTING_SUBTYPE, & & PROBLEM_STANDARD_DATA_FITTING_SUBTYPE, & & PROBLEM_VECTOR_DATA_FITTING_SUBTYPE, & & PROBLEM_DIV_FREE_VECTOR_DATA_FITTING_SUBTYPE, & @@ -5779,6 +5870,8 @@ SUBROUTINE Fitting_PreSolve(solver,err,error,*) SELECT CASE(problem%specification(3)) CASE(PROBLEM_STATIC_FITTING_SUBTYPE) !Do nothing + CASE(PROBLEM_STATIC_NONLINEAR_FITTING_SUBTYPE) + !Do nothing CASE(PROBLEM_STANDARD_DATA_FITTING_SUBTYPE) !Do nothing CASE(PROBLEM_GENERALISED_DATA_FITTING_SUBTYPE) @@ -5860,6 +5953,8 @@ SUBROUTINE Fitting_PostSolve(solver,err,error,*) SELECT CASE(problem%specification(3)) CASE(PROBLEM_STATIC_FITTING_SUBTYPE) !Do nothing + CASE(PROBLEM_STATIC_NONLINEAR_FITTING_SUBTYPE) + !Do nothing CASE(PROBLEM_STANDARD_DATA_FITTING_SUBTYPE,PROBLEM_GENERALISED_DATA_FITTING_SUBTYPE, & & PROBLEM_MAT_PROPERTIES_DATA_FITTING_SUBTYPE) !Do nothing @@ -5939,6 +6034,8 @@ SUBROUTINE Fitting_PostSolveOutputData(solver,err,error,*) SELECT CASE(problem%specification(3)) CASE(PROBLEM_STATIC_FITTING_SUBTYPE) !Do nothing + CASE(PROBLEM_STATIC_NONLINEAR_FITTING_SUBTYPE) + !Do nothing CASE(PROBLEM_STANDARD_DATA_FITTING_SUBTYPE,PROBLEM_GENERALISED_DATA_FITTING_SUBTYPE, & & PROBLEM_MAT_PROPERTIES_DATA_FITTING_SUBTYPE, & & PROBLEM_DATA_POINT_VECTOR_STATIC_FITTING_SUBTYPE) @@ -6054,6 +6151,8 @@ SUBROUTINE Fitting_PreSolveUpdateInputData(solver,err,error,*) SELECT CASE(problem%specification(3)) CASE(PROBLEM_STATIC_FITTING_SUBTYPE) !Do nothing + CASE(PROBLEM_STATIC_NONLINEAR_FITTING_SUBTYPE) + !Do nothing CASE(PROBLEM_STANDARD_DATA_FITTING_SUBTYPE) !Do nothing CASE(PROBLEM_GENERALISED_DATA_FITTING_SUBTYPE) diff --git a/src/opencmiss_iron.f90 b/src/opencmiss_iron.f90 index 13ae51c8..72b485ce 100644 --- a/src/opencmiss_iron.f90 +++ b/src/opencmiss_iron.f90 @@ -5917,6 +5917,8 @@ MODULE OpenCMISS_Iron INTEGER(INTG), PARAMETER :: CMFE_PROBLEM_STATIC_FITTING_SUBTYPE = & & PROBLEM_STATIC_FITTING_SUBTYPE ! Date: Tue, 5 Jun 2018 12:03:48 +1200 Subject: [PATCH 09/18] Removing linear and nonlinear prefix from fitting equation sets --- src/equations_set_routines.f90 | 4 ++-- src/fitting_routines.f90 | 26 +++++++++++++------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/equations_set_routines.f90 b/src/equations_set_routines.f90 index f99098f3..0d5a598a 100644 --- a/src/equations_set_routines.f90 +++ b/src/equations_set_routines.f90 @@ -2598,7 +2598,7 @@ SUBROUTINE EquationsSet_FiniteElementCalculate(equationsSet,elementNumber,err,er CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) CALL CLASSICAL_FIELD_FINITE_ELEMENT_CALCULATE(equationsSet,elementNumber,err,error,*999) CASE(EQUATIONS_SET_FITTING_CLASS) - CALL LinearFitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,*999) + CALL Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,*999) CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) IF(SIZE(equationsSet%specification,1)<2) & & CALL FlagError("Equations set specification must have at least two entries for a bioelectrics equation class.", & @@ -3015,7 +3015,7 @@ SUBROUTINE EquationsSet_FiniteElementResidualEvaluate(equationsSet,elementNumber CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) CALL ClassicalField_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*999) CASE(EQUATIONS_SET_FITTING_CLASS) - CALL NonlinearFitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*999) + CALL Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*999) CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) CALL FlagError("Not implemented.",err,error,*999) CASE(EQUATIONS_SET_MODAL_CLASS) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 22786c57..23c7ea27 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -97,8 +97,8 @@ MODULE FittingRoutines PUBLIC Fitting_ProblemSetup PUBLIC Fitting_ProblemSpecificationSet - PUBLIC LinearFitting_FiniteElementCalculate - PUBLIC NonlinearFitting_FiniteElementResidualEvaluate + PUBLIC Fitting_FiniteElementCalculate + PUBLIC Fitting_FiniteElementResidualEvaluate PUBLIC Fitting_PreSolve PUBLIC Fitting_PostSolve @@ -110,7 +110,7 @@ MODULE FittingRoutines ! !>Calculates the element stiffness matrices and RHS for a Galerkin projection finite element equations set. - SUBROUTINE LinearFitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,*) + SUBROUTINE Fitting_FiniteElementCalculate(equationsSet,elementNumber,err,error,*) !Argument variables TYPE(EQUATIONS_SET_TYPE), POINTER :: equationsSet !Evaluates the element residual and rhs vectors for the given element number for a fitting class finite element equation set. - SUBROUTINE NonlinearFitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*) + SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err,error,*) !Argument variables TYPE(EQUATIONS_SET_TYPE), POINTER :: equationsSet ! Date: Sat, 16 Jun 2018 20:45:51 +1200 Subject: [PATCH 10/18] Adding code use for linear fitting as a place holder to Fitting_FiniteElementResidualEvaluate --- src/fitting_routines.f90 | 336 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 323 insertions(+), 13 deletions(-) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 23c7ea27..75a12751 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -1348,31 +1348,341 @@ SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err, INTEGER(INTG), INTENT(OUT) :: err !equationsSet%EQUATIONS - IF(ASSOCIATED(EQUATIONS)) THEN + equations=>equationsSet%equations + IF(ASSOCIATED(equations)) THEN NULLIFY(vectorEquations) CALL Equations_VectorEquationsGet(equations,vectorEquations,err,error,*999) - SELECT CASE(equationsSet%SPECIFICATION(3)) - CASE(EQUATIONS_SET_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE) + IF(.NOT.ALLOCATED(equationsSet%specification)) & + & CALL FlagError("Equations set specification is not allocated.",err,error,*999) + IF(SIZE(equationsSet%specification,1)/=4) & + & CALL FlagError("Equations set specification must have four entries for a fitting type equations set.",err,error,*999) + smoothingType=equationsSet%specification(4) + SELECT CASE(equationsSet%specification(2)) + CASE(EQUATIONS_SET_DATA_FITTING_EQUATION_TYPE) + + SELECT CASE(equationsSet%specification(3)) + CASE(EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE) + geometricField=>equations%interpolation%geometricField + dependentField=>equations%interpolation%dependentField + independentField=>equations%interpolation%independentField + vectorMatrices=>vectorEquations%vectorMatrices + linearMatrices=>vectorMatrices%linearMatrices + equationsMatrix=>linearMatrices%matrices(1)%ptr + rhsVector=>vectorMatrices%rhsVector + vectorMapping=>vectorEquations%vectorMapping + linearMapping=>vectorMapping%linearMapping + dependentVariable=>linearMapping%equationsMatrixToVarMaps(1)%VARIABLE + dependentVariableType=dependentVariable%VARIABLE_TYPE + dataVariable=>independentField%VARIABLE_TYPE_MAP(FIELD_U_VARIABLE_TYPE)%ptr + dataWeightVariable=>independentField%VARIABLE_TYPE_MAP(FIELD_V_VARIABLE_TYPE)%ptr + dependentBasis=>dependentField%decomposition%domain(dependentField%decomposition%MESH_COMPONENT_NUMBER)%ptr% & + & topology%elements%elements(elementNumber)%basis + geometricBasis=>geometricField%decomposition%domain(geometricField%decomposition%MESH_COMPONENT_NUMBER)%ptr% & + & topology%elements%elements(elementNumber)%basis + quadratureScheme=>dependentBasis%quadrature%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr + CALL Field_InterpolationParametersElementGet(FIELD_VALUES_SET_TYPE,elementNumber,equations%interpolation% & + & geometricInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + CALL Field_InterpolationParametersElementGet(FIELD_VALUES_SET_TYPE,elementNumber,equations%interpolation% & + & dependentInterpParameters(dependentVariableType)%ptr,err,error,*999) + CALL Field_NumberOfComponentsGet(geometricField,FIELD_U_VARIABLE_TYPE,numberOfDimensions,err,error,*999) + CALL Field_NumberOfComponentsGet(independentField,FIELD_U_VARIABLE_TYPE,numberOfDataComponents,err,error,*999) + IF(numberOfDataComponents>99) CALL FlagError("Increase the size of the data point vectors.",err,error,*999) + numberOfXi = dependentBasis%NUMBER_OF_XI + + !Get data point vector parameters + CALL Field_ParameterSetDataGet(independentField,FIELD_U_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE, & + & independentVectorParameters,err,error,*999) + !Get data point weight parameters + CALL Field_ParameterSetDataGet(independentField,FIELD_V_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE, & + & independentWeightParameters,err,error,*999) + + !=============================== + ! D a t a P o i n t F i t + !=============================== + dataProjection=>independentField%dataProjection + IF(.NOT.ASSOCIATED(dataProjection)) & + & CALL FlagError("Data projection is not associated on independent field.",err,error,*999) + decompositionTopology=>independentField%decomposition%topology + IF(ASSOCIATED(decompositionTopology)) THEN + dataPoints=>decompositionTopology%dataPoints + IF(.NOT.ASSOCIATED(dataPoints)) & + & CALL FlagError("Data points are not associated on the decomposition topology of the independent field.", & + & err,error,*999) + ELSE + CALL FlagError("Decomposition topology is not associated on the independent field.",err,error,*999) + ENDIF + !Loop over data points + DO dataPointIdx=1,dataPoints%elementDataPoint(elementNumber)%numberOfProjectedData + dataPointUserNumber = dataPoints%elementDataPoint(elementNumber)%dataIndices(dataPointIdx)%userNumber + dataPointLocalNumber = dataPoints%elementDataPoint(elementNumber)%dataIndices(dataPointIdx)%localNumber + 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) + !Get data point vector value and weight + DO componentIdx=1,numberOfDataComponents + localDof=dataVariable%components(componentIdx)%PARAM_TO_DOF_MAP% & + & DATA_POINT_PARAM2DOF_MAP%DATA_POINTS(dataPointLocalNumber) + dataPointVector(componentIdx)=independentVectorParameters(localDof) + localDof=dataWeightVariable%components(componentIdx)%PARAM_TO_DOF_MAP% & + & DATA_POINT_PARAM2DOF_MAP%DATA_POINTS(dataPointLocalNumber) + dataPointWeight(componentIdx)=independentWeightParameters(localDof) + ENDDO !componentIdx + + dependentParameterRowIdx=0 + !Loop over element rows + DO dependentComponentRowIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + meshComponentRow=dependentVariable%components(dependentComponentRowIdx)%MESH_COMPONENT_NUMBER + dependentBasisRow=>dependentField%decomposition%domain(meshComponentRow)%ptr%topology%elements% & + & elements(elementNumber)%basis + DO dependentElementParameterRowIdx=1,dependentBasisRow%NUMBER_OF_ELEMENT_PARAMETERS + dependentParameterRowIdx=dependentParameterRowIdx+1 + dependentParameterColumnIdx=0 + basisFunctionRow=Basis_EvaluateXi(dependentBasisRow,dependentElementParameterRowIdx,NO_PART_DERIV, & + & projectionXi,err,error) + IF(equationsMatrix%updateMatrix) THEN + !Loop over element columns + DO dependentComponentColumnIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + meshComponentColumn=dependentVariable%components(dependentComponentColumnIdx)%MESH_COMPONENT_NUMBER + dependentBasisColumn=>dependentField%decomposition%domain(meshComponentColumn)%ptr% & + & topology%elements%elements(elementNumber)%basis + DO dependentElementParameterColumnIdx=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS + dependentParameterColumnIdx=dependentParameterColumnIdx+1 + !Treat each component as separate and independent so only calculate the diagonal blocks + IF(dependentComponentColumnIdx==dependentComponentRowIdx) THEN + basisFunctionColumn=Basis_EvaluateXi(dependentBasisColumn,dependentElementParameterColumnIdx, & + & NO_PART_DERIV,projectionXi,err,error) + sum = basisFunctionRow*basisFunctionColumn*dataPointWeight(dependentComponentRowIdx) + equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)= & + & equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)+sum + ENDIF + ENDDO !dependentElementParameterColumnIdx + ENDDO !dependentComponentColumnIdx + ENDIF + IF(rhsVector%updateVector) THEN + sum = basisFunctionRow*dataPointVector(dependentComponentRowIdx)*dataPointWeight(dependentComponentRowIdx) + rhsVector%elementVector%vector(dependentParameterRowIdx)= & + & rhsVector%elementVector%vector(dependentParameterRowIdx)+sum + ENDIF + ENDDO !dependentElementParameterRowIdx + ENDDO !dependentComponentRowIdx + ENDDO !dataPointIdx + + !Restore data point vector parameters + CALL Field_ParameterSetDataRestore(independentField,FIELD_U_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE, & + & independentVectorParameters,err,error,*999) + !Restore data point weight parameters + CALL Field_ParameterSetDataRestore(independentField,FIELD_V_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE, & + & independentWeightParameters,err,error,*999) + + SELECT CASE(smoothingType) + CASE(EQUATIONS_SET_FITTING_NO_SMOOTHING) + !Do nothing + CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_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) + !Loop over Gauss points + DO gaussPointIdx=1,quadratureScheme%NUMBER_OF_GAUSS + !Interpolate fields + CALL Field_InterpolateGauss(FIRST_PART_DERIV,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) + CALL Field_InterpolatedPointMetricsCalculate(geometricBasis%NUMBER_OF_XI,equations%interpolation% & + & geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + !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) + jacobianGaussWeight=equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr%jacobian* & + & quadratureScheme%GAUSS_WEIGHTS(gaussPointIdx) + + !Loop over field components + dependentParameterRowIdx=0 + DO dependentComponentRowIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + !Loop over element rows + meshComponentRow=dependentVariable%components(dependentComponentRowIdx)%MESH_COMPONENT_NUMBER + dependentBasisRow=>dependentField%decomposition%domain(meshComponentRow)%ptr% & + & topology%elements%elements(elementNumber)%basis + quadratureSchemeRow=>dependentBasisRow%quadrature%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr + DO dependentElementParameterRowIdx=1,dependentBasisRow%NUMBER_OF_ELEMENT_PARAMETERS + dependentParameterRowIdx=dependentParameterRowIdx+1 + dependentParameterColumnIdx=0 + !Loop over element columns + DO dependentComponentColumnIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + meshComponentColumn=dependentVariable%components(dependentComponentColumnIdx)%MESH_COMPONENT_NUMBER + dependentBasisColumn=>dependentField%decomposition%domain(meshComponentColumn)%ptr% & + & topology%elements%elements(elementNumber)%basis + quadratureSchemeColumn=>dependentBasisColumn%quadrature%QUADRATURE_SCHEME_MAP( & + & BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr + DO dependentElementParameterColumnIdx=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS + 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)) + 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)) + 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)) + ENDIF ! 3D + ENDIF ! 2 or 3D + sum = (tension + curvature) * jacobianGaussWeight + + equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)= & + equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)+sum + + ENDDO !dependentElementParameterColumnIdx + ENDDO !dependentComponentColumnIdx + 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 + localError="The fitting smoothing type of "//TRIM(NumberToVString(smoothingType,"*",err,error))// & + & " is invalid." + CALL FlagError(localError,err,error,*999) + END SELECT + + CASE DEFAULT + localError="Equations set subtype "//TRIM(NumberToVString(equationsSet%specification(3),"*",err,error))// & + & " is not valid for a Gauss fitting equations set class." + CALL FlagError(localError,err,error,*999) + END SELECT CASE DEFAULT - localError="Equations set subtype "//TRIM(NumberToVString(equationsSet%SPECIFICATION(3),"*",err,error))// & - & " is not valid for a fitting equation set class." + localError="Equations set type "//TRIM(NumberToVString(equationsSet%specification(2),"*",err,error))// & + & " is not valid for a fitting equations set class." CALL FlagError(localError,err,error,*999) END SELECT + + !Scale factor adjustment + IF(dependentField%SCALINGS%SCALING_TYPE/=FIELD_NO_SCALING) THEN + CALL Field_InterpolationParametersScaleFactorsElementGet(elementNumber,equations%interpolation% & + & dependentInterpParameters(dependentVariableType)%ptr,err,error,*999) + dependentParameterRowIdx=0 + !Loop over element rows + DO dependentComponentRowIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + meshComponentRow=dependentVariable%components(dependentComponentRowIdx)%MESH_COMPONENT_NUMBER + dependentBasisRow=>dependentField%decomposition%domain(meshComponentRow)%ptr% & + & topology%elements%elements(elementNumber)%basis + DO dependentElementParameterRowIdx=1,dependentBasisRow%NUMBER_OF_ELEMENT_PARAMETERS + dependentParameterRowIdx=dependentParameterRowIdx+1 + dependentParameterColumnIdx=0 + IF(equationsMatrix%updateMatrix) THEN + !Loop over element columns + DO dependentComponentColumnIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + meshComponentColumn=dependentVariable%components(dependentComponentColumnIdx)%MESH_COMPONENT_NUMBER + dependentBasisColumn=>dependentField%decomposition%domain(meshComponentColumn)%ptr% & + & topology%elements%elements(elementNumber)%basis + DO dependentElementParameterColumnIdx=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS + dependentParameterColumnIdx=dependentParameterColumnIdx+1 + equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)= & + & equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)* & + & equations%interpolation%dependentInterpParameters(dependentVariableType)%ptr% & + & SCALE_FACTORS(dependentElementParameterRowIdx,dependentComponentRowIdx)* & + & equations%interpolation%dependentInterpParameters(dependentVariableType)%ptr% & + & SCALE_FACTORS(dependentElementParameterColumnIdx,dependentComponentColumnIdx) + ENDDO !dependentElementParameterColumnIdx + ENDDO !dependentComponentColumnIdx + ENDIF + IF(rhsVector%updateVector) THEN + rhsVector%elementVector%vector(dependentParameterRowIdx)= & + & rhsVector%elementVector%vector(dependentParameterRowIdx)* & + & equations%interpolation%dependentInterpParameters(dependentVariableType)%ptr% & + & SCALE_FACTORS(dependentElementParameterRowIdx,dependentComponentRowIdx) + ENDIF + ENDDO !dependentElementParameterRowIdx + ENDDO !dependentComponentRowIdx + ENDIF ELSE CALL FlagError("Equations set equations is not associated.",err,error,*999) ENDIF From f7a7f9158938a03abe56c8a008a77f2d16d27b69 Mon Sep 17 00:00:00 2001 From: Prasad Date: Sat, 16 Jun 2018 21:39:35 +1200 Subject: [PATCH 11/18] Adding PROBLEM_FIBRE_FITTING_TYPE --- src/fitting_routines.f90 | 9 +++++++++ src/opencmiss_iron.f90 | 3 ++- src/problem_constants.f90 | 3 ++- 3 files changed, 13 insertions(+), 2 deletions(-) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 75a12751..2099c4f7 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -6077,6 +6077,15 @@ SUBROUTINE Fitting_ProblemSpecificationSet(problem,problemSpecification,err,erro & " is not valid for a Galerkin projection type of a data fitting problem." CALL FlagError(localError,err,error,*999) END SELECT + CASE(PROBLEM_FIBRE_FITTING_TYPE) + SELECT CASE(problemSubtype) + CASE(PROBLEM_STATIC_NONLINEAR_FITTING_SUBTYPE) + !ok + CASE DEFAULT + localError="The third problem specification of "//TRIM(NumberToVstring(problemSubtype,"*",err,error))// & + & " is not valid for a Galerkin projection type of a fibre fitting problem." + CALL FlagError(localError,err,error,*999) + END SELECT CASE DEFAULT localError="The second problem specification of "//TRIM(NumberToVstring(problemType,"*",err,error))// & & " is not valid for a data fitting problem." diff --git a/src/opencmiss_iron.f90 b/src/opencmiss_iron.f90 index 72b485ce..8d15daca 100644 --- a/src/opencmiss_iron.f90 +++ b/src/opencmiss_iron.f90 @@ -5824,6 +5824,7 @@ MODULE OpenCMISS_Iron INTEGER(INTG), PARAMETER :: CMFE_PROBLEM_BIDOMAIN_EQUATION_TYPE = PROBLEM_BIDOMAIN_EQUATION_TYPE ! Date: Sat, 16 Jun 2018 22:08:23 +1200 Subject: [PATCH 12/18] Bug fixes for diffusion tesnor fibre fitting problem and equation set setup --- src/fitting_routines.f90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 2099c4f7..055d2e7f 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -1398,7 +1398,7 @@ SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err, CASE(EQUATIONS_SET_DATA_FITTING_EQUATION_TYPE) SELECT CASE(equationsSet%specification(3)) - CASE(EQUATIONS_SET_DATA_POINT_FITTING_SUBTYPE) + CASE(EQUATIONS_SET_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE) geometricField=>equations%interpolation%geometricField dependentField=>equations%interpolation%dependentField independentField=>equations%interpolation%independentField @@ -1634,7 +1634,7 @@ SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err, CASE DEFAULT localError="Equations set subtype "//TRIM(NumberToVString(equationsSet%specification(3),"*",err,error))// & - & " is not valid for a Gauss fitting equations set class." + & " is not valid for a data fitting equations set class." CALL FlagError(localError,err,error,*999) END SELECT @@ -5499,6 +5499,8 @@ SUBROUTINE Fitting_ProblemSetup(problem,problemSetup,err,error,*) SELECT CASE(problem%specification(3)) CASE(PROBLEM_STATIC_FITTING_SUBTYPE) CALL Fitting_ProblemStaticSetup(problem,problemSetup,err,error,*999) + CASE(PROBLEM_STATIC_NONLINEAR_FITTING_SUBTYPE) + CALL Fitting_ProblemStaticSetup(problem,problemSetup,err,error,*999) CASE(PROBLEM_STANDARD_DATA_FITTING_SUBTYPE) CALL FITTING_PROBLEM_STANDARD_SETUP(problem,problemSetup,err,error,*999) CASE(PROBLEM_VECTOR_DATA_FITTING_SUBTYPE) @@ -6063,7 +6065,6 @@ SUBROUTINE Fitting_ProblemSpecificationSet(problem,problemSpecification,err,erro CASE(PROBLEM_DATA_FITTING_TYPE) SELECT CASE(problemSubtype) CASE(PROBLEM_STATIC_FITTING_SUBTYPE, & - & PROBLEM_STATIC_NONLINEAR_FITTING_SUBTYPE, & & PROBLEM_STANDARD_DATA_FITTING_SUBTYPE, & & PROBLEM_VECTOR_DATA_FITTING_SUBTYPE, & & PROBLEM_DIV_FREE_VECTOR_DATA_FITTING_SUBTYPE, & From 0c572d06b5cc1075ff10792d88dac86639248d58 Mon Sep 17 00:00:00 2001 From: Prasad Date: Thu, 21 Jun 2018 14:54:30 +1200 Subject: [PATCH 13/18] Nullifying indpendent field parameters pointers --- src/fitting_routines.f90 | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 055d2e7f..dc24200b 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -1355,8 +1355,10 @@ SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err, TYPE(EquationsType), POINTER :: equations TYPE(EquationsMappingVectorType), POINTER :: vectorMapping TYPE(EquationsMappingLinearType), POINTER :: linearMapping + TYPE(EquationsMappingNonlinearType), POINTER :: nonlinearMapping TYPE(EquationsMatricesVectorType), POINTER :: vectorMatrices TYPE(EquationsMatricesLinearType), POINTER :: linearMatrices + TYPE(EquationsMatricesNonlinearType), POINTER :: nonlinearMatrices TYPE(EquationsMatricesRHSType), POINTER :: rhsVector TYPE(EquationsMatrixType), POINTER :: equationsMatrix TYPE(EquationsVectorType), POINTER :: vectorEquations @@ -1403,12 +1405,15 @@ SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err, dependentField=>equations%interpolation%dependentField independentField=>equations%interpolation%independentField vectorMatrices=>vectorEquations%vectorMatrices - linearMatrices=>vectorMatrices%linearMatrices - equationsMatrix=>linearMatrices%matrices(1)%ptr + nonlinearMatrices=>vectorMatrices%nonlinearMatrices + !!!linearMatrices=>vectorMatrices%linearMatrices + !!!equationsMatrix=>linearMatrices%matrices(1)%ptr rhsVector=>vectorMatrices%rhsVector vectorMapping=>vectorEquations%vectorMapping - linearMapping=>vectorMapping%linearMapping - dependentVariable=>linearMapping%equationsMatrixToVarMaps(1)%VARIABLE + !!!linearMapping=>vectorMapping%linearMapping + nonlinearMapping=>vectorMapping%nonlinearMapping + dependentVariable=>nonlinearMapping%residualVariables(1)%ptr + !!!dependentVariable=>linearMapping%equationsMatrixToVarMaps(1)%VARIABLE dependentVariableType=dependentVariable%VARIABLE_TYPE dataVariable=>independentField%VARIABLE_TYPE_MAP(FIELD_U_VARIABLE_TYPE)%ptr dataWeightVariable=>independentField%VARIABLE_TYPE_MAP(FIELD_V_VARIABLE_TYPE)%ptr @@ -1427,9 +1432,11 @@ SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err, numberOfXi = dependentBasis%NUMBER_OF_XI !Get data point vector parameters + NULLIFY(independentVectorParameters) CALL Field_ParameterSetDataGet(independentField,FIELD_U_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE, & & independentVectorParameters,err,error,*999) !Get data point weight parameters + NULLIFY(independentWeightParameters) CALL Field_ParameterSetDataGet(independentField,FIELD_V_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE, & & independentWeightParameters,err,error,*999) From 322fdbda235d4ca983877441fa613e63ca6bd09a Mon Sep 17 00:00:00 2001 From: Prasad Date: Thu, 21 Jun 2018 14:55:22 +1200 Subject: [PATCH 14/18] Updating fibre fitting equation set setup --- src/fitting_routines.f90 | 156 +++++++++++++++++++++++++-------------- 1 file changed, 100 insertions(+), 56 deletions(-) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index dc24200b..19032dc8 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -3161,7 +3161,7 @@ SUBROUTINE Fitting_EquationsSetDataSetup(equationsSet,equationsSetSetup,err,erro TYPE(VARYING_STRING), INTENT(OUT) :: error ! Date: Thu, 21 Jun 2018 14:55:35 +1200 Subject: [PATCH 15/18] Updating fibre fitting residual evaluate --- src/fitting_routines.f90 | 41 +++++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 19 deletions(-) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 19032dc8..1dbaae9f 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -1489,25 +1489,28 @@ SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err, dependentParameterColumnIdx=0 basisFunctionRow=Basis_EvaluateXi(dependentBasisRow,dependentElementParameterRowIdx,NO_PART_DERIV, & & projectionXi,err,error) - IF(equationsMatrix%updateMatrix) THEN - !Loop over element columns - DO dependentComponentColumnIdx=1,dependentVariable%NUMBER_OF_COMPONENTS - meshComponentColumn=dependentVariable%components(dependentComponentColumnIdx)%MESH_COMPONENT_NUMBER - dependentBasisColumn=>dependentField%decomposition%domain(meshComponentColumn)%ptr% & - & topology%elements%elements(elementNumber)%basis - DO dependentElementParameterColumnIdx=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS - dependentParameterColumnIdx=dependentParameterColumnIdx+1 - !Treat each component as separate and independent so only calculate the diagonal blocks - IF(dependentComponentColumnIdx==dependentComponentRowIdx) THEN - basisFunctionColumn=Basis_EvaluateXi(dependentBasisColumn,dependentElementParameterColumnIdx, & - & NO_PART_DERIV,projectionXi,err,error) - sum = basisFunctionRow*basisFunctionColumn*dataPointWeight(dependentComponentRowIdx) - equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)= & - & equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)+sum - ENDIF - ENDDO !dependentElementParameterColumnIdx - ENDDO !dependentComponentColumnIdx - ENDIF + !!!IF(equationsMatrix%updateMatrix) THEN + !Loop over element columns + sum = 0.0_DP + DO dependentComponentColumnIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + meshComponentColumn=dependentVariable%components(dependentComponentColumnIdx)%MESH_COMPONENT_NUMBER + dependentBasisColumn=>dependentField%decomposition%domain(meshComponentColumn)%ptr% & + & topology%elements%elements(elementNumber)%basis + DO dependentElementParameterColumnIdx=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS + dependentParameterColumnIdx=dependentParameterColumnIdx+1 + !Treat each component as separate and independent so only calculate the diagonal blocks + IF(dependentComponentColumnIdx==dependentComponentRowIdx) THEN + basisFunctionColumn=Basis_EvaluateXi(dependentBasisColumn,dependentElementParameterColumnIdx, & + & NO_PART_DERIV,projectionXi,err,error) + sum = basisFunctionRow*basisFunctionColumn*dataPointWeight(dependentComponentRowIdx)* & + & equations%interpolation%dependentinterpparameters(1)%ptr%parameters( & + & dependentComponentRowIdx,dependentElementParameterRowIdx) + ENDIF + ENDDO !dependentElementParameterColumnIdx + ENDDO !dependentComponentColumnIdx + nonlinearMatrices%elementResidual%vector(dependentParameterRowIdx) = & + & nonlinearMatrices%elementResidual%vector(dependentParameterRowIdx) + sum + !!!ENDIF IF(rhsVector%updateVector) THEN sum = basisFunctionRow*dataPointVector(dependentComponentRowIdx)*dataPointWeight(dependentComponentRowIdx) rhsVector%elementVector%vector(dependentParameterRowIdx)= & From 44fe4891225f17216570dd7ce3c74d2c87a535c8 Mon Sep 17 00:00:00 2001 From: Prasad Date: Thu, 21 Jun 2018 20:28:53 +1200 Subject: [PATCH 16/18] Setting fibre fitting to use data_point_based_interpolation --- src/fitting_routines.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 1dbaae9f..6f975fef 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -4217,9 +4217,9 @@ SUBROUTINE Fitting_EquationsSetDataSetup(equationsSet,equationsSetSetup,err,erro CASE(EQUATIONS_SET_FEM_SOLUTION_METHOD) DO componentIdx=1,numberOfComponents CALL Field_ComponentInterpolationCheck(equationsSetSetup%FIELD,FIELD_U_VARIABLE_TYPE,componentIdx, & - & FIELD_GAUSS_POINT_BASED_INTERPOLATION,err,error,*999) + & FIELD_DATA_POINT_BASED_INTERPOLATION,err,error,*999) CALL Field_ComponentInterpolationCheck(equationsSetSetup%FIELD,FIELD_V_VARIABLE_TYPE,componentIdx, & - & FIELD_GAUSS_POINT_BASED_INTERPOLATION,err,error,*999) + & FIELD_DATA_POINT_BASED_INTERPOLATION,err,error,*999) ENDDO !componentIdx CASE DEFAULT localError="The solution method of "//TRIM(NumberToVString(equationsSet%SOLUTION_METHOD, & From 3f46a99e01b54885e91919abd5f48be41775d078 Mon Sep 17 00:00:00 2001 From: Prasad Date: Thu, 21 Jun 2018 21:46:11 +1200 Subject: [PATCH 17/18] Removing references to equationsMatrix/updateMatrix in fitting residual evaluate --- src/fitting_routines.f90 | 210 +++++++++++++++++++-------------------- 1 file changed, 104 insertions(+), 106 deletions(-) diff --git a/src/fitting_routines.f90 b/src/fitting_routines.f90 index 6f975fef..4e4d3f94 100644 --- a/src/fitting_routines.f90 +++ b/src/fitting_routines.f90 @@ -1502,9 +1502,9 @@ SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err, IF(dependentComponentColumnIdx==dependentComponentRowIdx) THEN basisFunctionColumn=Basis_EvaluateXi(dependentBasisColumn,dependentElementParameterColumnIdx, & & NO_PART_DERIV,projectionXi,err,error) - sum = basisFunctionRow*basisFunctionColumn*dataPointWeight(dependentComponentRowIdx)* & - & equations%interpolation%dependentinterpparameters(1)%ptr%parameters( & - & dependentComponentRowIdx,dependentElementParameterRowIdx) + !sum = basisFunctionRow*basisFunctionColumn*dataPointWeight(dependentComponentRowIdx)* & + !& equations%interpolation%dependentinterpparameters(1)%ptr%parameters( & + !& dependentComponentRowIdx,dependentElementParameterRowIdx) ENDIF ENDDO !dependentElementParameterColumnIdx ENDDO !dependentComponentColumnIdx @@ -1531,106 +1531,105 @@ SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err, CASE(EQUATIONS_SET_FITTING_NO_SMOOTHING) !Do nothing CASE(EQUATIONS_SET_FITTING_SOBOLEV_VALUE_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) - !Loop over Gauss points - DO gaussPointIdx=1,quadratureScheme%NUMBER_OF_GAUSS - !Interpolate fields - CALL Field_InterpolateGauss(FIRST_PART_DERIV,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) - CALL Field_InterpolatedPointMetricsCalculate(geometricBasis%NUMBER_OF_XI,equations%interpolation% & - & geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) - !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) - jacobianGaussWeight=equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr%jacobian* & - & quadratureScheme%GAUSS_WEIGHTS(gaussPointIdx) - - !Loop over field components - dependentParameterRowIdx=0 - DO dependentComponentRowIdx=1,dependentVariable%NUMBER_OF_COMPONENTS - !Loop over element rows - meshComponentRow=dependentVariable%components(dependentComponentRowIdx)%MESH_COMPONENT_NUMBER - dependentBasisRow=>dependentField%decomposition%domain(meshComponentRow)%ptr% & - & topology%elements%elements(elementNumber)%basis - quadratureSchemeRow=>dependentBasisRow%quadrature%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr - DO dependentElementParameterRowIdx=1,dependentBasisRow%NUMBER_OF_ELEMENT_PARAMETERS - dependentParameterRowIdx=dependentParameterRowIdx+1 - dependentParameterColumnIdx=0 - !Loop over element columns - DO dependentComponentColumnIdx=1,dependentVariable%NUMBER_OF_COMPONENTS - meshComponentColumn=dependentVariable%components(dependentComponentColumnIdx)%MESH_COMPONENT_NUMBER - dependentBasisColumn=>dependentField%decomposition%domain(meshComponentColumn)%ptr% & - & topology%elements%elements(elementNumber)%basis - quadratureSchemeColumn=>dependentBasisColumn%quadrature%QUADRATURE_SCHEME_MAP( & - & BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr - DO dependentElementParameterColumnIdx=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS - dependentParameterColumnIdx=dependentParameterColumnIdx+1 + materialsField=>equations%interpolation%materialsField + CALL Field_InterpolationParametersElementGet(FIELD_VALUES_SET_TYPE,elementNumber,equations%interpolation% & + & materialsInterpParameters(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + !Loop over Gauss points + DO gaussPointIdx=1,quadratureScheme%NUMBER_OF_GAUSS + !Interpolate fields + CALL Field_InterpolateGauss(FIRST_PART_DERIV,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) + CALL Field_InterpolatedPointMetricsCalculate(geometricBasis%NUMBER_OF_XI,equations%interpolation% & + & geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr,err,error,*999) + !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) + jacobianGaussWeight=equations%interpolation%geometricInterpPointMetrics(FIELD_U_VARIABLE_TYPE)%ptr%jacobian* & + & quadratureScheme%GAUSS_WEIGHTS(gaussPointIdx) + + !Loop over field components + dependentParameterRowIdx=0 + DO dependentComponentRowIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + !Loop over element rows + meshComponentRow=dependentVariable%components(dependentComponentRowIdx)%MESH_COMPONENT_NUMBER + dependentBasisRow=>dependentField%decomposition%domain(meshComponentRow)%ptr% & + & topology%elements%elements(elementNumber)%basis + quadratureSchemeRow=>dependentBasisRow%quadrature%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr + DO dependentElementParameterRowIdx=1,dependentBasisRow%NUMBER_OF_ELEMENT_PARAMETERS + dependentParameterRowIdx=dependentParameterRowIdx+1 + dependentParameterColumnIdx=0 + !Loop over element columns + DO dependentComponentColumnIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + meshComponentColumn=dependentVariable%components(dependentComponentColumnIdx)%MESH_COMPONENT_NUMBER + dependentBasisColumn=>dependentField%decomposition%domain(meshComponentColumn)%ptr% & + & topology%elements%elements(elementNumber)%basis + quadratureSchemeColumn=>dependentBasisColumn%quadrature%QUADRATURE_SCHEME_MAP( & + & BASIS_DEFAULT_QUADRATURE_SCHEME)%ptr + DO dependentElementParameterColumnIdx=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS + dependentParameterColumnIdx=dependentParameterColumnIdx+1 - !Calculate Sobolev surface tension and curvature smoothing terms - tension = tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1, & + !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)) + 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_S1, & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S2, & & gaussPointIdx)) - curvature = kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S1, & + curvature = curvature + kappaParam*2.0_DP* ( & + & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2_S2, & & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S1_S1, & + & 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)) - IF(numberOfXi > 1) THEN + IF(numberOfXi > 2) THEN tension = tension + tauParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2, & + & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S3, & & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S2, & + & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S3, & & gaussPointIdx)) curvature = curvature + kappaParam*2.0_DP* ( & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S2_S2, & + & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S3_S3, & & gaussPointIdx)* & - & quadratureSchemeColumn%GAUSS_BASIS_FNS(dependentElementParameterColumnIdx,PART_DERIV_S2_S2, & - & gaussPointIdx) + & - & quadratureSchemeRow%GAUSS_BASIS_FNS(dependentElementParameterRowIdx,PART_DERIV_S1_S2, & + & 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_S2, & + & 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)) - 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)) - ENDIF ! 3D - ENDIF ! 2 or 3D - sum = (tension + curvature) * jacobianGaussWeight + ENDIF ! 3D + ENDIF ! 2 or 3D + sum = (tension + curvature) * jacobianGaussWeight - equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)= & - equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)+sum + !!! Update residual vector + !!!equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)= & + !!! equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)+sum - ENDDO !dependentElementParameterColumnIdx - ENDDO !dependentComponentColumnIdx - ENDDO !dependentElementParameterRowIdx - ENDDO !dependentComponentRowIdx - ENDDO !gaussPointIdx - ENDIF + ENDDO !dependentElementParameterColumnIdx + ENDDO !dependentComponentColumnIdx + ENDDO !dependentElementParameterRowIdx + ENDDO !dependentComponentRowIdx + ENDDO !gaussPointIdx CASE(EQUATIONS_SET_FITTING_SOBOLEV_DIFFERENCE_SMOOTHING) CALL FlagError("Not implemented.",err,error,*999) @@ -1667,23 +1666,22 @@ SUBROUTINE Fitting_FiniteElementResidualEvaluate(equationsSet,elementNumber,err, DO dependentElementParameterRowIdx=1,dependentBasisRow%NUMBER_OF_ELEMENT_PARAMETERS dependentParameterRowIdx=dependentParameterRowIdx+1 dependentParameterColumnIdx=0 - IF(equationsMatrix%updateMatrix) THEN - !Loop over element columns - DO dependentComponentColumnIdx=1,dependentVariable%NUMBER_OF_COMPONENTS - meshComponentColumn=dependentVariable%components(dependentComponentColumnIdx)%MESH_COMPONENT_NUMBER - dependentBasisColumn=>dependentField%decomposition%domain(meshComponentColumn)%ptr% & - & topology%elements%elements(elementNumber)%basis - DO dependentElementParameterColumnIdx=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS - dependentParameterColumnIdx=dependentParameterColumnIdx+1 - equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)= & - & equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)* & - & equations%interpolation%dependentInterpParameters(dependentVariableType)%ptr% & - & SCALE_FACTORS(dependentElementParameterRowIdx,dependentComponentRowIdx)* & - & equations%interpolation%dependentInterpParameters(dependentVariableType)%ptr% & - & SCALE_FACTORS(dependentElementParameterColumnIdx,dependentComponentColumnIdx) - ENDDO !dependentElementParameterColumnIdx - ENDDO !dependentComponentColumnIdx - ENDIF + !Loop over element columns + DO dependentComponentColumnIdx=1,dependentVariable%NUMBER_OF_COMPONENTS + meshComponentColumn=dependentVariable%components(dependentComponentColumnIdx)%MESH_COMPONENT_NUMBER + dependentBasisColumn=>dependentField%decomposition%domain(meshComponentColumn)%ptr% & + & topology%elements%elements(elementNumber)%basis + DO dependentElementParameterColumnIdx=1,dependentBasisColumn%NUMBER_OF_ELEMENT_PARAMETERS + dependentParameterColumnIdx=dependentParameterColumnIdx+1 + !!! Add scale factors to residual vector + !equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)= & + ! & equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)* & + ! & equations%interpolation%dependentInterpParameters(dependentVariableType)%ptr% & + ! & SCALE_FACTORS(dependentElementParameterRowIdx,dependentComponentRowIdx)* & + ! & equations%interpolation%dependentInterpParameters(dependentVariableType)%ptr% & + ! & SCALE_FACTORS(dependentElementParameterColumnIdx,dependentComponentColumnIdx) + ENDDO !dependentElementParameterColumnIdx + ENDDO !dependentComponentColumnIdx IF(rhsVector%updateVector) THEN rhsVector%elementVector%vector(dependentParameterRowIdx)= & & rhsVector%elementVector%vector(dependentParameterRowIdx)* & From 61ff08c57ca6b01b880d7e6e58f178e541c73c6c Mon Sep 17 00:00:00 2001 From: Prasad Date: Thu, 21 Jun 2018 21:53:38 +1200 Subject: [PATCH 18/18] Updating pre and post solve residual evaluate for nonlinear fitting problems --- src/problem_routines.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/problem_routines.f90 b/src/problem_routines.f90 index 5121a2ea..6bddc7b3 100644 --- a/src/problem_routines.f90 +++ b/src/problem_routines.f90 @@ -1474,6 +1474,8 @@ SUBROUTINE PROBLEM_PRE_RESIDUAL_EVALUATE(SOLVER,err,error,*) !Pre residual evaluate not used CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) !Pre residual evaluate not used + CASE(PROBLEM_FITTING_CLASS) + !do nothing CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) !Pre residual evaluate not used CASE(EQUATIONS_SET_MODAL_CLASS) @@ -1615,6 +1617,8 @@ SUBROUTINE PROBLEM_POST_RESIDUAL_EVALUATE(SOLVER,err,error,*) !Post residual evaluate not used CASE(EQUATIONS_SET_CLASSICAL_FIELD_CLASS) !Post residual evaluate not used + CASE(PROBLEM_FITTING_CLASS) + !do nothing CASE(EQUATIONS_SET_BIOELECTRICS_CLASS) !Post residual evaluate not used CASE(EQUATIONS_SET_MODAL_CLASS)