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)Evaluates the element residual and rhs vectors for the given element number for a fitting class finite element equation set. + SUBROUTINE Fitting_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) + 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_DIFFUSION_TENSOR_FIBRE_FITTING_SUBTYPE) + geometricField=>equations%interpolation%geometricField + dependentField=>equations%interpolation%dependentField + independentField=>equations%interpolation%independentField + vectorMatrices=>vectorEquations%vectorMatrices + nonlinearMatrices=>vectorMatrices%nonlinearMatrices + !!!linearMatrices=>vectorMatrices%linearMatrices + !!!equationsMatrix=>linearMatrices%matrices(1)%ptr + rhsVector=>vectorMatrices%rhsVector + vectorMapping=>vectorEquations%vectorMapping + !!!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 + 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 + 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) + + !=============================== + ! 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 + 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)= & + & 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) + 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 + + !!! Update residual vector + !!!equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)= & + !!! equationsMatrix%elementMatrix%matrix(dependentParameterRowIdx,dependentParameterColumnIdx)+sum + + ENDDO !dependentElementParameterColumnIdx + ENDDO !dependentComponentColumnIdx + ENDDO !dependentElementParameterRowIdx + ENDDO !dependentComponentRowIdx + ENDDO !gaussPointIdx + + 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 data fitting equations set class." + CALL FlagError(localError,err,error,*999) + END SELECT + + CASE DEFAULT + 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 + !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)* & + & 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 + ELSE + CALL FlagError("Equations set is not associated.",err,error,*999) + ENDIF + + EXITS("Fitting_FiniteElementResidualEvaluate") + RETURN +999 ERRORS("Fitting_FiniteElementResidualEvaluate",err,error) + EXITS("Fitting_FiniteElementResidualEvaluate") + RETURN 1 + + END SUBROUTINE Fitting_FiniteElementResidualEvaluate + + ! + !================================================================================================================================ + ! + !>Sets up the update-materials Galerkin projection. SUBROUTINE FITTING_EQUATIONS_SET_MAT_PROPERTIES_SETUP(equationsSet,equationsSetSetup,err,error,*) @@ -1950,7 +2322,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 +2394,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 +2492,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, & @@ -2749,67 +3123,623 @@ SUBROUTINE Fitting_EquationsSetGaussSetup(equationsSet,equationsSetSetup,err,err CALL FlagError(localError,err,error,*999) END SELECT CASE DEFAULT - localError="The setup type of "//TRIM(NumberToVString(equationsSetSetup%SETUP_TYPE,"*",err,error))// & - & " is invalid for a fitting equations set." + localError="The 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 DEFAULT + localError="The equations set subtype of "//TRIM(NumberToVString(equationsSet%specification(3),"*",err,error))// & + & " is invalid." + CALL FlagError(localError,err,error,*999) + END SELECT + ELSE + localError="The equations set type of "//TRIM(NumberToVString(equationsSet%specification(2),"*",err,error))// & + & " does not equal a Gauss fitting type." + CALL FlagError(localError,err,error,*999) + ENDIF + ELSE + CALL FlagError("Equations set is not associated.",err,error,*999) + ENDIF + + EXITS("Fitting_EquationsSetGaussSetup") + RETURN +999 ERRORSEXITS("Fitting_EquationsSetGaussSetup",err,error) + RETURN 1 + + END SUBROUTINE Fitting_EquationsSetGaussSetup + + ! + !================================================================================================================================ + ! + + !>Sets up the standard data point fitting equations set. + SUBROUTINE Fitting_EquationsSetDataSetup(equationsSet,equationsSetSetup,err,error,*) + + !Argument variables + TYPE(EQUATIONS_SET_TYPE), POINTER :: equationsSet !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_LINEAR,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_U_VARIABLE_TYPE,vectorMapping,err,error,*999) + CALL EquationsMapping_LinearMatricesNumberSet(vectorMapping,1,err,error,*999) + CALL EquationsMapping_LinearMatricesVariableTypesSet(vectorMapping,[FIELD_U_VARIABLE_TYPE], & + & 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) + SELECT CASE(equations%sparsityType) + CASE(EQUATIONS_MATRICES_FULL_MATRICES) + CALL EquationsMatrices_LinearStorageTypeSet(vectorMatrices,[MATRIX_BLOCK_STORAGE_TYPE],err,error,*999) + CASE(EQUATIONS_MATRICES_SPARSE_MATRICES) + CALL EquationsMatrices_LinearStorageTypeSet(vectorMatrices,[MATRIX_COMPRESSED_ROW_STORAGE_TYPE], & + & err,error,*999) + CALL EquationsMatrices_LinearStructureTypeSet(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 + 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 - localError="The equations set subtype of "//TRIM(NumberToVString(equationsSet%specification(3),"*",err,error))// & - & " is invalid." + !----------------------------------------------------------------- + ! 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 - ELSE - localError="The equations set type of "//TRIM(NumberToVString(equationsSet%specification(2),"*",err,error))// & - & " does not equal a Gauss fitting type." - CALL FlagError(localError,err,error,*999) - ENDIF - ELSE - CALL FlagError("Equations set is not associated.",err,error,*999) - ENDIF - - EXITS("Fitting_EquationsSetGaussSetup") - RETURN -999 ERRORSEXITS("Fitting_EquationsSetGaussSetup",err,error) - RETURN 1 - - END SUBROUTINE Fitting_EquationsSetGaussSetup - - ! - !================================================================================================================================ - ! - - !>Sets up the standard data point fitting equations set. - SUBROUTINE Fitting_EquationsSetDataSetup(equationsSet,equationsSetSetup,err,error,*) - - !Argument variables - TYPE(EQUATIONS_SET_TYPE), POINTER :: equationsSet !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) @@ -5059,6 +6130,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." @@ -5162,6 +6242,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) @@ -5243,6 +6325,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 @@ -5322,6 +6406,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) @@ -5437,6 +6523,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 ed0b281f..8d15daca 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 !