OpenCMISS-Iron Internal API Documentation
interface_operators_routines.f90
Go to the documentation of this file.
1 
43 
45 MODULE interface_operators_routines
46 
47  USE base_routines
48  USE basis_routines
49  USE constants
50  USE field_routines
51  USE input_output
52  USE interface_conditions_constants
53  USE interface_equations_routines
54  USE interface_mapping_routines
55  USE interface_matrices_routines
56  USE iso_varying_string
57  USE kinds
58  USE matrix_vector
59  USE strings
60  USE timer
61  USE types
62 
63 #include "macros.h"
64 
65  IMPLICIT NONE
66 
67  !Module types
68 
69  !Module variables
70 
71  !Interfaces
72 
73  PUBLIC fieldcontinuity_finiteelementcalculate
74 
75  PUBLIC frictionlesscontact_finiteelementcalculate
76 
77  PUBLIC solidfluidoperator_finiteelementcalculate
78 
79 CONTAINS
80 
81  !
82  !================================================================================================================================
83  !
84 
86  SUBROUTINE fieldcontinuity_finiteelementcalculate(interfaceCondition,elementNumber,err,error,*)
87 
88  !Argument variables
89  TYPE(interface_condition_type), POINTER :: interfacecondition
90  TYPE(interface_equations_type), POINTER :: interfaceequations
91  INTEGER(INTG), INTENT(IN) :: elementnumber
92  INTEGER(INTG), INTENT(OUT) :: err
93  TYPE(varying_string), INTENT(OUT) :: error
94  !Local Variables
95  INTEGER(INTG) :: gausspoint, rowcomponentidx, rowidx, rowparameteridx, colcomponentidx, colidx, colparameteridx
96  INTEGER(INTG) :: rowmeshcomponentnumber,derivativeidx,derivative,localelementnode,interfacenode,interfacederivative
97  INTEGER(INTG) :: coupledmeshelementnumber,coupledmeshidx,coupledmeshvariabletype,lagrangevariabletype
98  INTEGER(INTG) :: connectedline,decompositionlinenumber,locallinenodeidx,connectedface,decompositionfacenumber,localfacenodeidx
99  REAL(DP) :: xi(3),rwg,pgmsi,pgnsi,matrixcoefficient
100  TYPE(basis_type), POINTER :: interfacedependentbasis,coupledmeshbasis,interfacegeometricbasis, &
101  & interfacePenaltyBasis,interfaceConnectivityBasis
102  TYPE(quadrature_scheme_type), POINTER :: interfacequadraturescheme
103  TYPE(field_type), POINTER :: coupledmeshdependentfield,interfacedependentfield,interfacegeometricfield, &
104  & interfacePenaltyField
105  TYPE(field_variable_type), POINTER :: interfacematrixvariable,lagrangevariable
106  TYPE(element_matrix_type), POINTER :: interfaceelementmatrix
107  TYPE(interface_equations_domain_interpolation_type), POINTER :: interfaceinterpolation
108  TYPE(interface_element_connectivity_type), POINTER :: elementconnectivity
109  TYPE(domain_line_type), POINTER :: coupledmeshdomainline
110  TYPE(domain_face_type), POINTER :: coupledmeshdomainface
111  TYPE(varying_string) :: localerror
112 
113  TYPE(interface_type), POINTER :: interface
114  TYPE(interfacepointsconnectivitytype), POINTER :: pointsconnectivity
115  TYPE(decompositionelementdatapointstype), POINTER :: decompositionelementdata
116  TYPE(basis_type), POINTER :: coupledmeshdependentbasis
117  TYPE(field_type), POINTER :: coupledmeshgeometricfield
118  INTEGER(INTG) :: meshcomponentnumber,numberofcoupledmeshgeocomp,numberofinterfacemeshxi,numberofcoupledmeshxi, &
119  & numberOfMatrixCoupledElements
120  INTEGER(INTG) :: datapointidx,localelementnumber,matrixelementidx
121  INTEGER(INTG) :: matrixcoefficients(2),interfaceelementnumber
122 
123  enters("FieldContinuity_FiniteElementCalculate",err,error,*999)
124 
125  IF(.NOT.ASSOCIATED(interfacecondition)) CALL flag_error("Interface condition is not associated.",err,error,*999)
126  IF(.NOT.ASSOCIATED(interfacecondition%INTERFACE_EQUATIONS)) CALL flag_error("Interface equations is not associated." &
127  & ,err,error,*999)
128  IF(.NOT.ASSOCIATED(interfacecondition%INTERFACE)) CALL flag_error("Interface is not associated.",err,error,*999)
129 
130  interfaceequations=>interfacecondition%INTERFACE_EQUATIONS
131 
132  SELECT CASE(interfacecondition%METHOD)
133 
134  CASE(interface_condition_point_to_point_method)
135  CALL flag_error("Not implemented.",err,error,*999)
136 
137  CASE(interface_condition_lagrange_multipliers_method,interface_condition_penalty_method)
138 
139  SELECT CASE(interfacecondition%integrationType)
140 
141  CASE(interface_condition_gauss_integration)
142  !Pointers to interface variables (columns of interface element matrix)
143  interfaceinterpolation=>interfaceequations%INTERPOLATION%INTERFACE_INTERPOLATION
144  interfacegeometricfield=>interfaceinterpolation%GEOMETRIC_FIELD
145  interfacedependentfield=>interfaceinterpolation%DEPENDENT_FIELD
146  interfacegeometricbasis=>interfacegeometricfield%DECOMPOSITION%DOMAIN(interfacegeometricfield% &
147  & decomposition%MESH_COMPONENT_NUMBER)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(elementnumber)%BASIS
148  interfacedependentbasis=>interfacedependentfield%DECOMPOSITION%DOMAIN(interfacedependentfield% &
149  & decomposition%MESH_COMPONENT_NUMBER)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(elementnumber)%BASIS
150  SELECT CASE(interfacecondition%METHOD)
151  CASE(interface_condition_penalty_method)
152  interfacepenaltyfield=>interfaceinterpolation%PENALTY_FIELD
153  interfacepenaltybasis=>interfacepenaltyfield%DECOMPOSITION%DOMAIN(interfacepenaltyfield% &
154  & decomposition%MESH_COMPONENT_NUMBER)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(elementnumber)%BASIS
155  CALL field_interpolation_parameters_element_get(field_values_set_type,elementnumber,interfaceinterpolation% &
156  & penalty_interpolation(1)%INTERPOLATION_PARAMETERS(field_u_variable_type)%PTR,err,error,*999)
157  ENDSELECT
158  !Integrate using the interface quadrature scheme
159  interfacequadraturescheme=>interfacegeometricbasis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
160  lagrangevariable=>interfaceequations%INTERFACE_MAPPING%LAGRANGE_VARIABLE
161  lagrangevariabletype=lagrangevariable%VARIABLE_TYPE
162  !Get element interpolation parameters from the first geometric interpolation set (to get Jacobian for interface surface integral)
163  CALL field_interpolation_parameters_element_get(field_values_set_type,elementnumber,interfaceinterpolation% &
164  & geometric_interpolation(1)%INTERPOLATION_PARAMETERS(field_u_variable_type)%PTR,err,error,*999)
165  !\todo Use matrix coefficients in solver routines when assembling solver matrices instead of multiplying them here
166  matrixcoefficient=1.0_dp
167  DO coupledmeshidx=1,interfaceequations%INTERFACE_MATRICES%NUMBER_OF_INTERFACE_MATRICES
168  IF(interfaceequations%INTERFACE_MATRICES%MATRICES(coupledmeshidx)%PTR%UPDATE_MATRIX) THEN
169  !\todo Use matrix coefficients in solver routines when assembling solver matrices instead of multiplying them here
170  IF(coupledmeshidx>1) THEN
171  matrixcoefficient=-1.0_dp
172  ENDIF
173  !Pointers to the coupledMeshIdx'th coupled mesh variables (rows of interface element matrix)
174  coupledmeshdependentfield=>interfaceequations%INTERPOLATION%VARIABLE_INTERPOLATION(coupledmeshidx)%DEPENDENT_FIELD
175  elementconnectivity=>interfacecondition%INTERFACE%MESH_CONNECTIVITY%ELEMENT_CONNECTIVITY(elementnumber,coupledmeshidx)
176  coupledmeshelementnumber=elementconnectivity%COUPLED_MESH_ELEMENT_NUMBER
177  interfacematrixvariable=> &
178  & interfaceequations%INTERFACE_MAPPING%INTERFACE_MATRIX_ROWS_TO_VAR_MAPS(coupledmeshidx)%VARIABLE
179  coupledmeshvariabletype=interfacematrixvariable%VARIABLE_TYPE
180  interfaceelementmatrix=>interfaceequations%INTERFACE_MATRICES%MATRICES(coupledmeshidx)%PTR%ELEMENT_MATRIX
181  interfaceconnectivitybasis=>interfacecondition%INTERFACE%MESH_CONNECTIVITY%BASIS
182 
183  !coupledMeshDependentInterpolation=>interfaceEquations%INTERPOLATION%VARIABLE_INTERPOLATION(coupledMeshIdx)% &
184  ! & DEPENDENT_INTERPOLATION
185 
186  !Loop over gauss points
187  DO gausspoint=1,interfacequadraturescheme%NUMBER_OF_GAUSS
188  !CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,GaussPoint, &
189  ! & coupledMeshDependentInterpolation%GEOMETRIC_INTERPOLATION(1)%INTERPOLATED_POINT(FIELD_U_VARIABLE_TYPE)%PTR, &
190  ! & err,error,*999)
191  !CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(interfaceGeometricBasis%NUMBER_OF_XI,interfaceInterpolation% &
192  ! & GEOMETRIC_INTERPOLATION(1)%INTERPOLATED_POINT_METRICS(FIELD_U_VARIABLE_TYPE)%PTR,err,error,*999)
193 
194  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,gausspoint,interfaceinterpolation% &
195  & geometric_interpolation(1)%INTERPOLATED_POINT(field_u_variable_type)%PTR,err,error,*999)
196  CALL field_interpolated_point_metrics_calculate(interfacegeometricbasis%NUMBER_OF_XI,interfaceinterpolation% &
197  & geometric_interpolation(1)%INTERPOLATED_POINT_METRICS(field_u_variable_type)%PTR,err,error,*999)
198  rwg=interfaceinterpolation%GEOMETRIC_INTERPOLATION(1)%INTERPOLATED_POINT_METRICS(field_u_variable_type)%PTR% &
199  & jacobian*interfacequadraturescheme%GAUSS_WEIGHTS(gausspoint)
200  IF(interfacecondition%METHOD==interface_condition_penalty_method .AND. &
201  & coupledmeshidx==interfaceequations%INTERFACE_MATRICES%NUMBER_OF_INTERFACE_MATRICES) THEN
202  CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,gausspoint,interfaceinterpolation% &
203  & penalty_interpolation(1)%INTERPOLATED_POINT(field_u_variable_type)%PTR,err,error,*999)
204  rowidx=0
205  DO rowcomponentidx=1,lagrangevariable%NUMBER_OF_COMPONENTS
206  !Loop over the Lagrange variable matrix rows
207  DO rowparameteridx=1,interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS
208  pgnsi=interfacequadraturescheme%GAUSS_BASIS_FNS(rowparameteridx,no_part_deriv,gausspoint)
209  rowidx=rowidx+1
210  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,rowidx)- &
211  & (1.0_dp/interfaceinterpolation%PENALTY_INTERPOLATION(1)% &
212  & interpolated_point(field_u_variable_type)%PTR%VALUES(1,1))*pgnsi**2.0_dp*rwg
213  ENDDO !rowParameterIdx
214  ENDDO !rowComponentIdx
215  ELSE
216  !\todo defaults to first mesh component, generalise
217  !XI=InterfaceOperators_InterfToCoupledMeshGaussTransform( &
218  ! & elementConnectivity,interfaceConnectivityBasis,GaussPoint,err,error)
219  xi(1:interfacedependentbasis%NUMBER_OF_XI)=interfaceoperators_interftocoupledmeshgausstransform( &
220  & elementconnectivity,interfaceconnectivitybasis,gausspoint,err,error)
221  !XI=interfaceCondition%interface%pointsConnectivity%pointsConnectivity(GaussPoint,coupledMeshIdx)%xi
222  ! Loop over number of Lagrange variable components as not all components in the dependent field variable may be coupled
223  !\todo Currently Lagrange field variable component numbers must match each coupled dependent field variable component numbers. Generalise ordering
224  DO rowcomponentidx=1,lagrangevariable%NUMBER_OF_COMPONENTS
225  rowmeshcomponentnumber=interfacematrixvariable%COMPONENTS(rowcomponentidx)%MESH_COMPONENT_NUMBER
226  coupledmeshbasis=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(rowmeshcomponentnumber)%PTR%TOPOLOGY% &
227  & elements%ELEMENTS(coupledmeshelementnumber)%BASIS
228 
229  SELECT CASE(interfacedependentbasis%NUMBER_OF_XI)
230 
231  CASE(1) !1D interface (line)
232  connectedline = elementconnectivity%CONNECTED_LINE
233  decompositionlinenumber=coupledmeshdependentfield%DECOMPOSITION%TOPOLOGY% &
234  & elements%ELEMENTS(coupledmeshelementnumber)%ELEMENT_LINES(connectedline)
235  coupledmeshdomainline=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(rowmeshcomponentnumber)%PTR%TOPOLOGY% &
236  & lines%LINES(decompositionlinenumber)
237  DO locallinenodeidx=1,coupledmeshbasis%NUMBER_OF_NODES_IN_LOCAL_LINE(connectedline)
238  localelementnode=coupledmeshbasis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx,connectedline)
239  DO derivativeidx=1,coupledmeshdomainline%BASIS%NUMBER_OF_DERIVATIVES(locallinenodeidx)
240  derivative=coupledmeshbasis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx,connectedline)
241  derivative=coupledmeshdomainline%DERIVATIVES_IN_LINE(1,derivativeidx,locallinenodeidx)
242  rowparameteridx=coupledmeshbasis%ELEMENT_PARAMETER_INDEX(derivative,localelementnode)
243  pgmsi=basis_evaluate_xi(coupledmeshbasis,rowparameteridx,no_part_deriv,xi,err,error)
244  rowidx=rowparameteridx+coupledmeshbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
245  DO interfacenode=1,interfacedependentbasis%NUMBER_OF_NODES
246  DO interfacederivative=1,interfacedependentbasis%NUMBER_OF_DERIVATIVES(interfacenode)
247  !\todo requires equal number of nodes between interface mesh and coupled mesh. Generalize
248  colparameteridx=interfacedependentbasis%ELEMENT_PARAMETER_INDEX(interfacederivative,interfacenode)
249  pgnsi=interfacequadraturescheme%GAUSS_BASIS_FNS(colparameteridx,no_part_deriv,gausspoint)
250  colidx=colparameteridx+interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
251  !\todo Use matrix coefficients in solver routines when assembling solver matrices instead of multiplying them here
252  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx)+ &
253  & pgnsi*pgmsi*rwg*matrixcoefficient
254  ENDDO !interfaceDerivative
255  ENDDO !interfaceNode
256  ENDDO !derivativeIdx
257  ENDDO !localLineNodeIdx
258 
259  CASE(2) !2D interface (face)
260 
261  SELECT CASE(coupledmeshbasis%NUMBER_OF_XI)
262 
263  CASE(2) !Coupled Mesh has 2 xi directions
264  DO localelementnode=1,coupledmeshbasis%NUMBER_OF_NODES
265  DO derivative=1,coupledmeshbasis%NUMBER_OF_DERIVATIVES(localelementnode)
266  rowparameteridx=coupledmeshbasis%ELEMENT_PARAMETER_INDEX(derivative,localelementnode)
267  pgmsi=basis_evaluate_xi(coupledmeshbasis,rowparameteridx,no_part_deriv, &
268  & xi(1:coupledmeshbasis%NUMBER_OF_XI),err,error)
269  rowidx=rowparameteridx+coupledmeshbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
270  DO interfacenode=1,interfacedependentbasis%NUMBER_OF_NODES
271  DO interfacederivative=1,interfacedependentbasis%NUMBER_OF_DERIVATIVES(interfacenode)
272  !\todo requires equal number of nodes between interface mesh and coupled mesh. Generalize
273  colparameteridx=interfacedependentbasis%ELEMENT_PARAMETER_INDEX(interfacederivative,interfacenode)
274  pgnsi=interfacequadraturescheme%GAUSS_BASIS_FNS(colparameteridx,no_part_deriv,gausspoint)
275  colidx=colparameteridx+interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
276  !\todo Use matrix coefficients in solver routines when assembling solver matrices instead of multiplying them here
277  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx)+ &
278  & pgnsi*pgmsi*rwg*matrixcoefficient
279  ENDDO !interfaceDerivative
280  ENDDO !interfaceNode
281  ENDDO !derivative
282  ENDDO !localElementNode
283 
284  CASE(3) !Coupled Mesh has 3 xi directions
285  connectedface = elementconnectivity%CONNECTED_FACE
286  decompositionfacenumber=coupledmeshdependentfield%DECOMPOSITION%TOPOLOGY% &
287  & elements%ELEMENTS(coupledmeshelementnumber)%ELEMENT_FACES(connectedface)
288  coupledmeshdomainface=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(rowmeshcomponentnumber)%PTR%TOPOLOGY% &
289  & faces%FACES(decompositionfacenumber)
290  DO localfacenodeidx=1,coupledmeshbasis%NUMBER_OF_NODES_IN_LOCAL_FACE(connectedface)
291  localelementnode=coupledmeshbasis%NODE_NUMBERS_IN_LOCAL_FACE(localfacenodeidx,connectedface)
292  DO derivativeidx=1,coupledmeshdomainface%BASIS%NUMBER_OF_DERIVATIVES(localfacenodeidx)
293  derivative=coupledmeshbasis% &
294  & derivative_numbers_in_local_face(derivativeidx,localfacenodeidx,connectedface)
295  rowparameteridx=coupledmeshbasis%ELEMENT_PARAMETER_INDEX(derivative,localelementnode)
296  pgmsi=basis_evaluate_xi(coupledmeshbasis,rowparameteridx,no_part_deriv, &
297  & xi(1:coupledmeshbasis%NUMBER_OF_XI),err,error)
298  rowidx=rowparameteridx+coupledmeshbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
299  DO interfacenode=1,interfacedependentbasis%NUMBER_OF_NODES
300  DO interfacederivative=1,interfacedependentbasis%NUMBER_OF_DERIVATIVES(interfacenode)
301  !\todo requires equal number of nodes between interface mesh and coupled mesh. Generalize
302  colparameteridx=interfacedependentbasis%ELEMENT_PARAMETER_INDEX(interfacederivative,interfacenode)
303  pgnsi=interfacequadraturescheme%GAUSS_BASIS_FNS(colparameteridx,no_part_deriv,gausspoint)
304  colidx=colparameteridx+interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
305  !\todo Use matrix coefficients in solver routines when assembling solver matrices instead of multiplying them here
306  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx)+ &
307  & pgnsi*pgmsi*rwg*matrixcoefficient
308  ENDDO !interfaceDerivative
309  ENDDO !interfaceNode
310  ENDDO !derivativeIdx
311  ENDDO !FaceNodeIdx
312 
313  END SELECT !coupledMeshBasis%NUMBER_OF_XI
314 
315  END SELECT !interfaceDependentBasis%NUMBER_OF_XI
316 
317  ENDDO !rowComponentIdx
318  ENDIF
319  ENDDO !GaussPoint
320 
321  !Scale factor adjustment
322  !\todo check if scale factor adjustments are already made elsewhere eg when calculating the interface matrix contribution to the residual for non-linear problems
323  !\todo update looping of variables/components for non-zero matrix elements as done above
324  IF(interfacecondition%METHOD==interface_condition_penalty_method .AND. &
325  & coupledmeshidx==interfaceequations%INTERFACE_MATRICES%NUMBER_OF_INTERFACE_MATRICES) THEN
326  !Scale factor adjustment for the Lagrange Variable (columns)
327  IF(interfacedependentfield%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
328  CALL field_interpolationparametersscalefactorselementget(elementnumber, &
329  & interfaceinterpolation%DEPENDENT_INTERPOLATION(1)%INTERPOLATION_PARAMETERS(lagrangevariabletype)%PTR, &
330  & err,error,*999)
331  rowidx=0
332  !Use Lagrange variable number of components here since we are only dealing with Lagrange variable scale factors
333  !\todo Currently Lagrange field variable component numbers must match each coupled dependent field variable component numbers. Generalise ordering
334  DO rowcomponentidx=1,lagrangevariable%NUMBER_OF_COMPONENTS
335  !Loop over element Lagrange variable rows
336  DO rowparameteridx=1,interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS
337  rowidx=rowidx+1
338  interfaceelementmatrix%MATRIX(rowidx,rowidx)=interfaceelementmatrix%MATRIX(rowidx,rowidx) * &
339  & interfaceinterpolation%DEPENDENT_INTERPOLATION(1)% &
340  & interpolation_parameters(lagrangevariabletype)%PTR%SCALE_FACTORS(rowparameteridx,rowcomponentidx)**2
341  ENDDO !rowParameterIdx
342  ENDDO !rowComponentIdx
343  ENDIF
344  ELSE
345  !Scale factor adjustment for the Lagrange Variable (columns)
346  IF(interfacedependentfield%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
347  CALL field_interpolationparametersscalefactorselementget(elementnumber, &
348  & interfaceinterpolation%DEPENDENT_INTERPOLATION(1)%INTERPOLATION_PARAMETERS(lagrangevariabletype)%PTR, &
349  & err,error,*999)
350  rowidx=0
351  !Use Lagrange variable number of components here since we are only dealing with Lagrange variable scale factors
352  !\todo Currently Lagrange field variable component numbers must match each coupled dependent field variable component numbers. Generalise ordering
353  DO rowcomponentidx=1,lagrangevariable%NUMBER_OF_COMPONENTS
354  rowmeshcomponentnumber=interfacematrixvariable%COMPONENTS(rowcomponentidx)%MESH_COMPONENT_NUMBER
355  coupledmeshbasis=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(rowmeshcomponentnumber)%PTR%TOPOLOGY% &
356  & elements%ELEMENTS(coupledmeshelementnumber)%BASIS
357  !Loop over element rows
358  DO rowparameteridx=1,coupledmeshbasis%NUMBER_OF_ELEMENT_PARAMETERS
359  rowidx=rowidx+1
360  colidx=0
361  !Loop over element columns
362  DO colcomponentidx=1,lagrangevariable%NUMBER_OF_COMPONENTS
363  DO colparameteridx=1,interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS
364  colidx=colidx+1
365  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx) * &
366  & interfaceinterpolation%DEPENDENT_INTERPOLATION(1)% &
367  & interpolation_parameters(lagrangevariabletype)%PTR%SCALE_FACTORS(colparameteridx,colcomponentidx)
368  ENDDO !colParameterIdx
369  ENDDO !colComponentIdx
370  ENDDO !rowParameterIdx
371  ENDDO !rowComponentIdx
372  ENDIF
373  !Scale factor adjustment for the row dependent variable
374  IF(coupledmeshdependentfield%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
375  CALL field_interpolationparametersscalefactorselementget(coupledmeshelementnumber, &
376  & interfaceequations%INTERPOLATION%VARIABLE_INTERPOLATION(coupledmeshidx)% &
377  & dependent_interpolation(1)%INTERPOLATION_PARAMETERS(coupledmeshvariabletype)%PTR,err,error,*999)
378  rowidx=0
379  DO rowcomponentidx=1,interfacematrixvariable%NUMBER_OF_COMPONENTS
380  !Loop over element rows
381  DO rowparameteridx=1,coupledmeshbasis%NUMBER_OF_ELEMENT_PARAMETERS
382  rowidx=rowidx+1
383  colidx=0
384  !Loop over element columns
385  DO colcomponentidx=1,lagrangevariable%NUMBER_OF_COMPONENTS
386  DO colparameteridx=1,interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS
387  colidx=colidx+1
388  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx)* &
389  & interfaceequations%INTERPOLATION%VARIABLE_INTERPOLATION(coupledmeshidx)% &
390  & dependent_interpolation(1)%INTERPOLATION_PARAMETERS(coupledmeshvariabletype)%PTR% &
391  & scale_factors(rowparameteridx,rowcomponentidx)
392  ENDDO !colParameterIdx
393  ENDDO !colComponentIdx
394  ENDDO !rowParameterIdx
395  ENDDO !rowComponentIdx
396  ENDIF
397  ENDIF
398  ENDIF
399  ENDDO ! coupledMeshIdx
400 
401  CASE(interface_condition_data_points_integration)
402 
403  matrixcoefficients(1)=1; !\todo: Change to interface mapping matrix coefficients
404  matrixcoefficients(2)=-1;
405  interfaceelementnumber = elementnumber!todo simplify
406  interface=>interfacecondition%INTERFACE
407  pointsconnectivity=>interface%pointsConnectivity
408  numberofinterfacemeshxi=pointsconnectivity%interfaceMesh%NUMBER_OF_DIMENSIONS
409  IF(ASSOCIATED(pointsconnectivity)) THEN
410  decompositionelementdata=>interfacecondition%LAGRANGE%LAGRANGE_FIELD%DECOMPOSITION%TOPOLOGY%dataPoints% &
411  & elementdatapoint(interfaceelementnumber)
412 
413  !Calculate PGSMI, update interface matrices with PGSMI, and update scale factors
414  DO coupledmeshidx=1,interface%NUMBER_OF_COUPLED_MESHES
415  IF(interfaceequations%INTERFACE_MATRICES%MATRICES(coupledmeshidx)%PTR%UPDATE_MATRIX) THEN
416  numberofmatrixcoupledelements=pointsconnectivity%coupledElements(interfaceelementnumber,coupledmeshidx)% &
417  & numberofcoupledelements
418  numberofcoupledmeshxi=interface%COUPLED_MESHES(coupledmeshidx)%PTR%NUMBER_OF_DIMENSIONS
419  coupledmeshgeometricfield=>interfacecondition%DEPENDENT%EQUATIONS_SETS(coupledmeshidx)%PTR% &
420  & geometry%GEOMETRIC_FIELD
421  coupledmeshdependentfield=>interfacecondition%DEPENDENT%EQUATIONS_SETS(coupledmeshidx)%PTR% &
422  & dependent%DEPENDENT_FIELD
423 
424  numberofcoupledmeshgeocomp=coupledmeshgeometricfield%VARIABLES(field_u_variable_type)%NUMBER_OF_COMPONENTS
425  interfaceelementmatrix=>interfaceequations%INTERFACE_MATRICES%MATRICES(coupledmeshidx)%PTR%ELEMENT_MATRIX
426  !mesh component number is the same for all geometric components in elasticity problems
427  meshcomponentnumber=coupledmeshdependentfield%VARIABLES(field_u_variable_type)%COMPONENTS(1)% &
428  & mesh_component_number
429  DO datapointidx=1,decompositionelementdata%numberOfProjectedData
430  localelementnumber=pointsconnectivity%pointsConnectivity(datapointidx,coupledmeshidx)% &
431  & coupledmeshelementnumber
432 
433  !Calculate the element index (non-conforming element) for this interface matrix
434  matrixelementidx=1
435  DO WHILE ((localelementnumber/=pointsconnectivity%coupledElements(interfaceelementnumber,coupledmeshidx)% &
436  & elementnumbers(matrixelementidx)).AND.(matrixelementidx/= &
437  & pointsconnectivity%coupledElements(interfaceelementnumber,coupledmeshidx)%numberOfCoupledElements))
438  matrixelementidx=matrixelementidx+1
439  ENDDO
440  xi(1:numberofcoupledmeshxi)=pointsconnectivity%pointsConnectivity(datapointidx,coupledmeshidx)% &
441  & xi(1:numberofcoupledmeshxi)
442  !Calculate PGSMI for each data point component
443  coupledmeshdependentbasis=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(meshcomponentnumber)%PTR% &
444  & topology%ELEMENTS%ELEMENTS(localelementnumber)%BASIS
445  DO rowcomponentidx=1,numberofcoupledmeshgeocomp
446  DO rowparameteridx=1,coupledmeshdependentbasis%NUMBER_OF_ELEMENT_PARAMETERS
447  pgmsi=basis_evaluate_xi(coupledmeshdependentbasis,rowparameteridx,no_part_deriv, &
448  & xi(1:numberofcoupledmeshxi),err,error)*matrixcoefficients(coupledmeshidx)
449  rowidx=rowparameteridx+coupledmeshdependentbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
450  colidx=datapointidx+decompositionelementdata%numberOfProjectedData*(rowcomponentidx-1)
451  interfaceelementmatrix%MATRIX(rowidx,colidx)=pgmsi !Update interface element matrix with contact point contribution
452  ENDDO !rowParameterIdx
453  ENDDO !rowComponentIdx
454  ENDDO !dataPointIdx
455 
456  !scale factor update
457  IF(coupledmeshdependentfield%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
458  DO datapointidx=1,decompositionelementdata%numberOfProjectedData
459  localelementnumber=pointsconnectivity%pointsConnectivity(datapointidx,coupledmeshidx)% &
460  & coupledmeshelementnumber
461  CALL field_interpolationparametersscalefactorselementget(localelementnumber,interfaceequations% &
462  & interpolation%VARIABLE_INTERPOLATION(coupledmeshidx)%DEPENDENT_INTERPOLATION(1)% &
463  & interpolation_parameters(field_u_variable_type)%PTR,err,error,*999)
464  !Calculate the element index (non-conforming element) for this interface matrix
465  matrixelementidx=1
466  DO WHILE ((localelementnumber/=pointsconnectivity%coupledElements(interfaceelementnumber,coupledmeshidx)% &
467  & elementnumbers(matrixelementidx)).AND.(matrixelementidx/= &
468  & pointsconnectivity%coupledElements(interfaceelementnumber,coupledmeshidx)%numberOfCoupledElements))
469  matrixelementidx=matrixelementidx+1
470  ENDDO
471  coupledmeshdependentbasis=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(meshcomponentnumber)%PTR% &
472  & topology%ELEMENTS%ELEMENTS(localelementnumber)%BASIS
473  DO rowcomponentidx=1,numberofcoupledmeshgeocomp
474  DO rowparameteridx=1,coupledmeshdependentbasis%NUMBER_OF_ELEMENT_PARAMETERS
475  rowidx=rowparameteridx+coupledmeshdependentbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
476  colidx=datapointidx+decompositionelementdata%numberOfProjectedData*(rowcomponentidx-1)
477  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx)* &
478  & interfaceequations%INTERPOLATION%VARIABLE_INTERPOLATION(coupledmeshidx)% &
479  & dependent_interpolation(1)%INTERPOLATION_PARAMETERS(field_u_variable_type)% &
480  & ptr%SCALE_FACTORS(rowparameteridx,rowcomponentidx)
481  ENDDO !rowParameterIdx
482  ENDDO !rowComponentIdx
483  ENDDO !dataPointIdx
484  ENDIF !.NOT. FIELD_NO_SCALING
485 
486  ENDIF !UPDATE_MATRIX
487  ENDDO !coupledMeshIdx
488  ELSE
489  CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
490  ENDIF
491 
492  CASE DEFAULT
493  localerror="Interface condition integration type "//trim(number_to_vstring(interfacecondition%integrationType, &
494  & "*",err,error))// " is not valid."
495  CALL flagerror(localerror,err,error,*999)
496  END SELECT !interfaceCondition%integrationType
497 
498  CASE(interface_condition_augmented_lagrange_method)
499  CALL flag_error("Not implemented.",err,error,*999)
500  CASE DEFAULT
501  localerror="Interface condition method "//trim(number_to_vstring(interfacecondition%METHOD,"*",err,error))// &
502  & " is not valid."
503  CALL flag_error(localerror,err,error,*999)
504  END SELECT
505 
506  exits("FieldContinuity_FiniteElementCalculate")
507  RETURN
508 999 errorsexits("FieldContinuity_FiniteElementCalculate",err,error)
509  RETURN 1
510 
511  END SUBROUTINE fieldcontinuity_finiteelementcalculate
512 
513  !
514  !================================================================================================================================
515  !
516 
518  SUBROUTINE frictionlesscontact_finiteelementcalculate(interfaceCondition,interfaceElementNumber,err,error,*)
519 
520  !Argument variables
521  TYPE(interface_condition_type), POINTER :: interfacecondition
522  INTEGER(INTG), INTENT(IN) :: interfaceelementnumber
523  INTEGER(INTG), INTENT(OUT) :: err
524  TYPE(varying_string), INTENT(OUT) :: error
525  !Local Variables
526  TYPE(interface_equations_type), POINTER :: interfaceequations
527  TYPE(interface_type), POINTER :: interface
528  TYPE(interfacepointsconnectivitytype), POINTER :: pointsconnectivity
529  TYPE(decompositionelementdatapointstype), POINTER :: decompositionelementdata
530  TYPE(field_type), POINTER :: coupledmeshdependentfield,penaltyfield
531  TYPE(interface_penalty_type), POINTER :: interfacepenalty
532  TYPE(field_interpolated_point_ptr_type), POINTER :: interpolatedpoints(:)
533  TYPE(field_interpolated_point_type), POINTER :: interpolatedpoint
534  TYPE(field_interpolation_parameters_ptr_type), POINTER :: interpolationparameters(:)
535  TYPE(field_interpolated_point_metrics_ptr_type), POINTER :: interpolatedpointsmetrics(:)
536  TYPE(basis_type), POINTER :: coupledmeshdependentbasis
537  TYPE(element_matrix_type), POINTER :: interfaceelementmatrix
538  TYPE(interface_matrix_type), POINTER :: penaltymatrix
539  INTEGER(INTG) :: meshcomponentnumber,numberofcoupledmeshgeocomp,numberofinterfacemeshxi,numberofcoupledmeshxi, &
540  & numberOfMatrixCoupledElements,localDof
541  INTEGER(INTG) :: datapointidx,coupledmeshidx,xiidx,localelementnumber,localfacelinenumber,matrixelementidx,rowcomponentidx, &
542  & rowParameterIdx,rowIdx,colIdx,componentIdx,globalDataPointNumber
543  INTEGER(INTG) :: matrixcoefficients(2)
544  REAL(DP) :: pgmsi,contactstiffness
545  REAL(DP) :: positionpoint(3),normalpoint(3),tangentspoint(3,3),xi(3)
546  LOGICAL :: reversenormal
547  REAL(DP), ALLOCATABLE :: gaps(:),gapscomponents(:,:),normals(:,:)
548  LOGICAL, ALLOCATABLE :: orthogonallyprojected(:)
549 
550 
551  TYPE(varying_string) :: localerror
552 
553  enters("FrictionlessContact_FiniteElementCalculate",err,error,*999)
554 
555  IF(ASSOCIATED(interfacecondition)) THEN
556  interfaceequations=>interfacecondition%INTERFACE_EQUATIONS
557  IF(ASSOCIATED(interfaceequations)) THEN
558  interface=>interfacecondition%INTERFACE
559  IF(ASSOCIATED(interface)) THEN
560  SELECT CASE(interfacecondition%METHOD)
561  CASE(interface_condition_point_to_point_method)
562  CALL flagerror("Not implemented.",err,error,*999)
563  CASE(interface_condition_lagrange_multipliers_method,interface_condition_penalty_method)
564  SELECT CASE(interfacecondition%integrationType)
565  CASE(interface_condition_gauss_integration)
566  CALL flagerror("Mesh connectivity is not implemented for frictionless contact.",err,error,*999)
567  CASE(interface_condition_data_points_integration)
568  matrixcoefficients(1)=1; !\todo: Change to interface mapping matrix coefficients
569  matrixcoefficients(2)=-1;
570  pointsconnectivity=>interface%pointsConnectivity
571  numberofinterfacemeshxi=pointsconnectivity%interfaceMesh%NUMBER_OF_DIMENSIONS
572  IF(ASSOCIATED(pointsconnectivity)) THEN
573  decompositionelementdata=>interfacecondition%LAGRANGE%LAGRANGE_FIELD%DECOMPOSITION%TOPOLOGY%dataPoints% &
574  & elementdatapoint(interfaceelementnumber)
575  !###################################################################################################################
576 
577  !Test if datapoints were orthogonally projected.
578  !\todo: Allow the user to choose to only include orthogonally projected points or not (check is commented when populating element matrix below).
579  ALLOCATE(orthogonallyprojected(decompositionelementdata%numberOfProjectedData),stat=err)
580  IF(err/=0) CALL flagerror("Could not allocate orthogonal projected logicals.",err,error,*999)
581  orthogonallyprojected=.true. !Initialise orthogonal projected logicals
582  DO coupledmeshidx=1,interface%NUMBER_OF_COUPLED_MESHES
583  coupledmeshdependentfield=>interfacecondition%DEPENDENT%EQUATIONS_SETS(coupledmeshidx)%PTR% &
584  & dependent%DEPENDENT_FIELD
585  !mesh component number is the same for all geometric components in elasticity problems
586  DO datapointidx=1,decompositionelementdata%numberOfProjectedData
587  globaldatapointnumber=decompositionelementdata%dataIndices(datapointidx)%globalNumber
588  DO xiidx=1,SIZE(pointsconnectivity%pointsConnectivity(globaldatapointnumber,coupledmeshidx)%reducedXi,1)
589  IF(abs(pointsconnectivity%pointsConnectivity(globaldatapointnumber,coupledmeshidx)%reducedXi(xiidx)) &
590  & < zero_tolerance) THEN
591  orthogonallyprojected(datapointidx)=.false.
592  ENDIF
593  ENDDO !xiIdx
594  ENDDO !dataPointIdx
595  ENDDO !coupledMeshIdx
596 
597  !###################################################################################################################
598 
599  !Allocate memory for local allocatable variables
600  ALLOCATE(gaps(decompositionelementdata%numberOfProjectedData),stat=err)
601  IF(err/=0) CALL flagerror("Could not allocate gaps.",err,error,*999)
602  gaps=0.0_dp !Initialise gap functions
603  ALLOCATE(gapscomponents(3,decompositionelementdata%numberOfProjectedData),stat=err)
604  IF(err/=0) CALL flagerror("Could not allocate component gaps.",err,error,*999)
605  gapscomponents=0.0_dp !Initialise gap functions
606  ALLOCATE(normals(3,decompositionelementdata%numberOfProjectedData),stat=err)
607  IF(err/=0) CALL flagerror("Could not allocate normals.",err,error,*999)
608  normals=0.0_dp !Initialise gap functions
609 
610  !Calculate Gap for each data point
611  !\todo: This is only required if only penetration is to penalized (ie seperation of meshes allowed.)
612  ! If a no seperation condition is also required then calculation of the gap is not required.
613  ! Need to allow user to choose which type of problem to solve.
614  DO coupledmeshidx=1,interface%NUMBER_OF_COUPLED_MESHES
615  coupledmeshdependentfield=>interfacecondition%DEPENDENT%FIELD_VARIABLES(coupledmeshidx)%PTR%FIELD
616  numberofcoupledmeshgeocomp=coupledmeshdependentfield%GEOMETRIC_FIELD%VARIABLE_TYPE_MAP(field_u_variable_type)% &
617  & ptr%NUMBER_OF_COMPONENTS
618  NULLIFY(interpolatedpoints)
619  NULLIFY(interpolationparameters)
620  CALL field_interpolation_parameters_initialise(coupledmeshdependentfield,interpolationparameters,err,error, &
621  & *999,field_geometric_components_type)
622  CALL field_interpolated_points_initialise(interpolationparameters,interpolatedpoints,err,error,*999, &
623  & field_geometric_components_type)
624  interpolatedpoint=>interpolatedpoints(field_u_variable_type)%PTR
625  DO datapointidx=1,decompositionelementdata%numberOfProjectedData
626  globaldatapointnumber=decompositionelementdata%dataIndices(datapointidx)%globalNumber
627  !Only interpolate if orthogonally projected
628  !\todo: Allow the user to choose to only include orthogonally projected points or not (currenlty commented out).
629  !IF(orthogonallyProjected(dataPointIdx)) THEN
630  localelementnumber=pointsconnectivity%pointsConnectivity(globaldatapointnumber,coupledmeshidx)% &
631  & coupledmeshelementnumber
632  localfacelinenumber=coupledmeshdependentfield%DECOMPOSITION%TOPOLOGY%ELEMENTS%ELEMENTS(localelementnumber)% &
633  & element_faces(pointsconnectivity%pointsConnectivity(globaldatapointnumber,coupledmeshidx)% &
634  & elementlinefacenumber)
635  SELECT CASE(numberofinterfacemeshxi) !Use face/line interpolation parameters for normal calculation
636  CASE(1)
637  CALL field_interpolation_parameters_line_get(field_values_set_type,localfacelinenumber, &
638  & interpolationparameters(field_u_variable_type)%PTR,err,error,*999,field_geometric_components_type)
639  CASE(2)
640  SELECT CASE(pointsconnectivity%pointsConnectivity(globaldatapointnumber,coupledmeshidx)% &
641  & elementlinefacenumber)
642  CASE(1,3,5)
643  reversenormal=.false.
644  CASE(2,4,6)
645  reversenormal=.true.
646  END SELECT
647  CALL field_interpolation_parameters_face_get(field_values_set_type,localfacelinenumber, &
648  & interpolationparameters(field_u_variable_type)%PTR,err,error,*999,field_geometric_components_type)
649  END SELECT
650  ! Determine the gap.
651  ! \todo: Note that FIELD_INTERPOLATE_XI(FIRST_PART_DERIV by default calculates NO_PART_DERIV too
652  ! and is used because the FIRST_PART_DERIV is also need for the surface normal calculation. However the
653  ! normal is only calculated for one of the coupled bodies so unnecessary computation. Need to generalize
654  ! FIELD_INTERPOLATE_XI to allow the user to specify which PART_DERIV to calculate.
655  CALL field_interpolate_xi(first_part_deriv,pointsconnectivity%pointsConnectivity(globaldatapointnumber, &
656  & coupledmeshidx)%reducedXi(:),interpolatedpoint,err,error,*999,field_geometric_components_type) !Interpolate contact data points on each surface
657  gapscomponents(1:numberofcoupledmeshgeocomp,datapointidx)=gapscomponents(1:numberofcoupledmeshgeocomp, &
658  & datapointidx)+interpolatedpoint%VALUES(1:numberofcoupledmeshgeocomp,no_part_deriv)* &
659  & matrixcoefficients(coupledmeshidx) !Calculate 3 components gap function for each contact point
660  !Calculate surface normal (use 2nd coupled mesh surface normal)
661  !\todo: Allow the user to choose which surface normal to calculate or alternatively allow for a weighted average of the two.
662  IF (coupledmeshidx==2) THEN
663  CALL field_interpolatedpointsmetricsinitialise(interpolatedpoints,interpolatedpointsmetrics, &
664  & err,error,*999)
665  CALL field_interpolated_point_metrics_calculate(numberofcoupledmeshgeocomp,interpolatedpointsmetrics &
666  & (field_u_variable_type)%PTR,err,error,*999)
667  CALL field_positionnormaltangentscalculateintptmetric(interpolatedpointsmetrics &
668  & (field_u_variable_type)%PTR,reversenormal,positionpoint,normalpoint,tangentspoint,err,error,*999)
669  normals(1:numberofcoupledmeshgeocomp,datapointidx)=normalpoint(1:numberofcoupledmeshgeocomp)
670  CALL field_interpolatedpointsmetricsfinalise(interpolatedpointsmetrics,err,error,*999)
671  ENDIF !coupledMeshIdx==1
672  !ENDIF !orthogonallyProjected(dataPointIdx)
673  ENDDO !dataPointIdx
674  CALL field_interpolation_parameters_finalise(interpolationparameters,err,error,*999)
675  CALL field_interpolated_points_finalise(interpolatedpoints,err,error,*999)
676  ENDDO !coupledMeshIdx
677 
678  !###################################################################################################################
679 
680  !Calcualte 1 component gap
681  DO datapointidx=1,decompositionelementdata%numberOfProjectedData
682  gaps(datapointidx)=dot_product(gapscomponents(1:numberofcoupledmeshgeocomp,datapointidx), &
683  & normals(1:numberofcoupledmeshgeocomp,datapointidx))
684  ENDDO !dataPointIdx
685 
686  !###################################################################################################################
687 
688  !Calculate PGSMI, update interface matrices with PGSMI, and update scale factors
689  DO coupledmeshidx=1,interface%NUMBER_OF_COUPLED_MESHES
690  IF(interfaceequations%INTERFACE_MATRICES%MATRICES(coupledmeshidx)%PTR%UPDATE_MATRIX) THEN
691  numberofmatrixcoupledelements=pointsconnectivity%coupledElements(interfaceelementnumber,coupledmeshidx)% &
692  & numberofcoupledelements
693  numberofcoupledmeshxi=interface%COUPLED_MESHES(coupledmeshidx)%PTR%NUMBER_OF_DIMENSIONS
694  numberofcoupledmeshgeocomp=interfacecondition%DEPENDENT%EQUATIONS_SETS(coupledmeshidx)%PTR%GEOMETRY% &
695  & geometric_field%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR%NUMBER_OF_COMPONENTS
696  coupledmeshdependentfield=>interfacecondition%DEPENDENT%EQUATIONS_SETS(coupledmeshidx)%PTR% &
697  & dependent%DEPENDENT_FIELD
698  interfaceelementmatrix=>interfaceequations%INTERFACE_MATRICES%MATRICES(coupledmeshidx)%PTR%ELEMENT_MATRIX
699  DO datapointidx=1,decompositionelementdata%numberOfProjectedData
700  globaldatapointnumber=decompositionelementdata%dataIndices(datapointidx)%globalNumber
701  !\todo: Allow the user to choose gap tolerance or default to zero tolerance (currently commented out).
702  !IF(gaps(dataPointIdx)>1.0E-10) THEN !Only add contact point contribution if the gap is a penetration
703  localelementnumber=pointsconnectivity%pointsConnectivity(globaldatapointnumber,coupledmeshidx)% &
704  & coupledmeshelementnumber
705  !Calculate the element index (non-conforming element) for this interface matrix
706  matrixelementidx=1
707  DO WHILE (localelementnumber/=pointsconnectivity%coupledElements(interfaceelementnumber,coupledmeshidx)% &
708  & elementnumbers(matrixelementidx))
709  matrixelementidx=matrixelementidx+1
710  ENDDO
711  CALL field_interpolationparametersscalefactorselementget(localelementnumber,interfaceequations% &
712  & interpolation%VARIABLE_INTERPOLATION(coupledmeshidx)%DEPENDENT_INTERPOLATION(1)% &
713  & interpolation_parameters(field_u_variable_type)%PTR,err,error,*999)
714  xi(1:numberofcoupledmeshxi)=pointsconnectivity%pointsConnectivity(globaldatapointnumber,coupledmeshidx)% &
715  & xi(1:numberofcoupledmeshxi)
716  DO rowcomponentidx=1,numberofcoupledmeshgeocomp
717  meshcomponentnumber=coupledmeshdependentfield%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR% &
718  & components(rowcomponentidx)%MESH_COMPONENT_NUMBER
719  !Calculate PGSMI for each data point component
720  coupledmeshdependentbasis=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(meshcomponentnumber)%PTR% &
721  & topology%ELEMENTS%ELEMENTS(localelementnumber)%BASIS
722  !\todo: Loop over the number of coupled mesh dependent basis element parameters on the contact face to save a bit of computation.
723  DO rowparameteridx=1,coupledmeshdependentbasis%NUMBER_OF_ELEMENT_PARAMETERS
724  pgmsi=basis_evaluate_xi(coupledmeshdependentbasis,rowparameteridx,no_part_deriv, &
725  & xi(1:numberofcoupledmeshxi),err,error)*normals(rowcomponentidx,datapointidx)* &
726  & matrixcoefficients(coupledmeshidx)
727  rowidx=coupledmeshdependentbasis%NUMBER_OF_ELEMENT_PARAMETERS*numberofmatrixcoupledelements* &
728  & (rowcomponentidx-1)+coupledmeshdependentbasis%NUMBER_OF_ELEMENT_PARAMETERS* &
729  & (matrixelementidx-1)+rowparameteridx
730  colidx=datapointidx
731  !Update interface element matrix with contact point contribution
732  !\todo: Seperate multiplication of scale factors if required.
733  interfaceelementmatrix%MATRIX(rowidx,colidx)=pgmsi*interfaceequations%INTERPOLATION% &
734  & variable_interpolation(coupledmeshidx)%DEPENDENT_INTERPOLATION(1)% &
735  & interpolation_parameters(field_u_variable_type)%PTR%SCALE_FACTORS(rowparameteridx,rowcomponentidx)
736  ENDDO !rowParameterIdx
737  ENDDO !rowComponentIdx
738  !ENDIF !gaps(dataPointIdx)>ZERO_TOLERANCE
739  ENDDO !dataPointIdx
740 
741  IF(interfaceequations%INTERFACE_MATRICES%MATRICES(coupledmeshidx)%PTR%FIRST_ASSEMBLY) &
742  & interfaceequations%INTERFACE_MATRICES%MATRICES(coupledmeshidx)%PTR%FIRST_ASSEMBLY=.false.
743  ENDIF !UPDATE_MATRIX
744  ENDDO !coupledMeshIdx
745 
746  !###################################################################################################################
747 
748  !Deallocate memory
749  IF(ALLOCATED(orthogonallyprojected)) DEALLOCATE(orthogonallyprojected)
750  IF(ALLOCATED(gapscomponents)) DEALLOCATE(gapscomponents)
751  IF(ALLOCATED(gaps)) DEALLOCATE(gaps)
752 
753  !###################################################################################################################
754 
755  !Calculate penalty matrix if required
756  IF(interfacecondition%METHOD==interface_condition_penalty_method) THEN
757  interfacepenalty=>interfacecondition%PENALTY
758  IF(ASSOCIATED(interfacepenalty)) THEN
759  penaltyfield=>interfacepenalty%PENALTY_FIELD
760  IF(ASSOCIATED(penaltyfield)) THEN
761  penaltymatrix=>interfaceequations%INTERFACE_MATRICES%MATRICES(interfaceequations% &
762  & interface_matrices%NUMBER_OF_INTERFACE_MATRICES)%PTR
763  IF(ASSOCIATED(penaltymatrix)) THEN
764  IF(penaltymatrix%FIRST_ASSEMBLY .AND. penaltymatrix%UPDATE_MATRIX) THEN
765  DO componentidx=1,penaltyfield%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR%NUMBER_OF_COMPONENTS
766  SELECT CASE(penaltyfield%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR% &
767  & components(componentidx)%INTERPOLATION_TYPE)
768  CASE(field_constant_interpolation)
769  CALL field_parameter_set_get_constant(penaltyfield,field_u_variable_type,field_values_set_type, &
770  & componentidx,contactstiffness,err,error,*999)
771  DO datapointidx=1,decompositionelementdata%numberOfProjectedData
772  penaltymatrix%ELEMENT_MATRIX%MATRIX(datapointidx,datapointidx)=-(1.0_dp/contactstiffness)
773  ENDDO !dataPointIdx
774  CASE(field_element_based_interpolation)
775  localdof=penaltyfield%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR%COMPONENTS(componentidx)% &
776  & param_to_dof_map%ELEMENT_PARAM2DOF_MAP%ELEMENTS(interfaceelementnumber)
777  CALL field_parameter_set_get_local_dof(penaltyfield,field_u_variable_type,field_values_set_type, &
778  & localdof,contactstiffness,err,error,*999)
779  DO datapointidx=1,decompositionelementdata%numberOfProjectedData
780  penaltymatrix%ELEMENT_MATRIX%MATRIX(datapointidx,datapointidx)=-(1.0_dp/contactstiffness)
781  ENDDO !dataPointIdx
782  CASE(field_data_point_based_interpolation)
783  DO datapointidx=1,decompositionelementdata%numberOfProjectedData
784  localdof=penaltyfield%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR%COMPONENTS(componentidx)% &
785  & param_to_dof_map%DATA_POINT_PARAM2DOF_MAP%DATA_POINTS(datapointidx)
786  CALL field_parameter_set_get_local_dof(penaltyfield,field_u_variable_type,field_values_set_type, &
787  & localdof,contactstiffness,err,error,*999)
788  penaltymatrix%ELEMENT_MATRIX%MATRIX(datapointidx,datapointidx)=-(1.0_dp/contactstiffness)
789  ENDDO !dataPointIdx
790  CASE DEFAULT
791  localerror="The interpolation type for component number "// &
792  & trim(number_to_vstring(componentidx,"*",err,error))// &
793  & " of variable type "//trim(number_to_vstring(field_u_variable_type,"*",err,error))// &
794  & " of field number "//trim(number_to_vstring(penaltyfield%USER_NUMBER,"*",err,error))//" is "// &
795  & trim(number_to_vstring(penaltyfield%VARIABLE_TYPE_MAP(field_u_variable_type)%PTR%COMPONENTS &
796  & (componentidx)%INTERPOLATION_TYPE,"*", err,error))// " which is invalid for penalty field."
797  CALL flagerror(localerror,err,error,*999)
798  END SELECT
799  ENDDO !componentIdx
800  ENDIF
801  ELSE
802  CALL flagerror("Interface penalty matrix is not associated.",err,error,*999)
803  ENDIF
804  ELSE
805  CALL flagerror("Interface penalty field is not associated.",err,error,*999)
806  ENDIF
807  ELSE
808  CALL flagerror("Interface penalty is not associated.",err,error,*999)
809  ENDIF
810  ENDIF
811  ELSE
812  CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
813  ENDIF
814  CASE DEFAULT
815  localerror="Interface condition integration type "//trim(number_to_vstring(interfacecondition%integrationType, &
816  & "*",err,error))// " is not valid."
817  CALL flagerror(localerror,err,error,*999)
818  END SELECT !interfaceCondition%integrationType
819  CASE(interface_condition_augmented_lagrange_method)
820  CALL flagerror("Not implemented.",err,error,*999)
821  CASE DEFAULT
822  localerror="Interface condition method "//trim(number_to_vstring(interfacecondition%METHOD,"*",err,error))// &
823  & " is not valid."
824  CALL flagerror(localerror,err,error,*999)
825  END SELECT !interfaceCondition%METHOD
826  ELSE
827  CALL flagerror("Interface is not associated.",err,error,*999)
828  ENDIF
829  ELSE
830  CALL flagerror("Interface equations is not associated.",err,error,*999)
831  ENDIF
832  ELSE
833  CALL flagerror("Interface condition is not associated.",err,error,*999)
834  ENDIF
835 
836  exits("FrictionlessContact_FiniteElementCalculate")
837  RETURN
838 999 errorsexits("FrictionlessContact_FiniteElementCalculate",err,error)
839  RETURN 1
840 
841  END SUBROUTINE frictionlesscontact_finiteelementcalculate
842 
843  !
844  !================================================================================================================================
845  !
846 
848  !Note First interface matrix must be fluid equations set's interface matrix
849  !SUBROUTINE FluidSolidOperator_FiniteElementCalculate(interfaceCondition,elementNumber,err,error,*)
850 
851  !
852  !================================================================================================================================
853  !
854 
856  !Note First interface matrix must be solid equations set's interface matrix
857  SUBROUTINE solidfluidoperator_finiteelementcalculate(interfaceCondition,elementNumber,err,error,*)
858 
859  !Argument variables
860  TYPE(interface_condition_type), POINTER :: interfacecondition
861  TYPE(interface_equations_type), POINTER :: interfaceequations
862  INTEGER(INTG), INTENT(IN) :: elementnumber
863  INTEGER(INTG), INTENT(OUT) :: err
864  TYPE(varying_string), INTENT(OUT) :: error
865  !Local Variables
866  INTEGER(INTG) :: gausspoint, rowcomponentidx, rowidx, rowparameteridx, colcomponentidx, colidx, colparameteridx
867  INTEGER(INTG) :: rowmeshcomponentnumber,derivativeidx,derivative,localelementnode,interfacenode,interfacederivative
868  INTEGER(INTG) :: coupledmeshelementnumber,coupledmeshidx,coupledmeshvariabletype,lagrangevariabletype
869  INTEGER(INTG) :: connectedline,decompositionlinenumber,locallinenodeidx,connectedface,decompositionfacenumber,localfacenodeidx
870  REAL(DP) :: xi(3),rwg,pgmsi,pgnsi,matrixcoefficient
871  TYPE(basis_type), POINTER :: interfacedependentbasis,coupledmeshbasis,interfacegeometricbasis, &
872  & interfaceConnectivityBasis
873  TYPE(quadrature_scheme_type), POINTER :: interfacequadraturescheme
874  TYPE(field_type), POINTER :: coupledmeshdependentfield,interfacedependentfield,interfacegeometricfield
875  TYPE(field_variable_type), POINTER :: interfacematrixvariable,lagrangevariable
876  TYPE(element_matrix_type), POINTER :: interfaceelementmatrix
877  TYPE(interface_equations_domain_interpolation_type), POINTER :: interfaceinterpolation
878  TYPE(interface_element_connectivity_type), POINTER :: elementconnectivity
879  TYPE(domain_line_type), POINTER :: coupledmeshdomainline
880  TYPE(domain_face_type), POINTER :: coupledmeshdomainface
881  TYPE(varying_string) :: localerror
882 
883  enters("SolidFluidOperator_FiniteElementCalculate",err,error,*999)
884 
885  IF(.NOT.ASSOCIATED(interfacecondition)) CALL flagerror("Interface condition is not associated.",err,error,*999)
886  IF(.NOT.ASSOCIATED(interfacecondition%INTERFACE_EQUATIONS)) CALL flagerror("Interface equations is not associated." &
887  & ,err,error,*999)
888  IF(.NOT.ASSOCIATED(interfacecondition%INTERFACE)) CALL flagerror("Interface is not associated.",err,error,*999)
889 
890  interfaceequations=>interfacecondition%INTERFACE_EQUATIONS
891 
892  !===============================================================================================================================
893  !Select Interface method
894  SELECT CASE(interfacecondition%METHOD)
895  CASE(interface_condition_point_to_point_method)
896  CALL flagerror("Not implemented.",err,error,*999)
897  CASE(interface_condition_lagrange_multipliers_method)
898  !=============================================================================================================================
899  !Select Integration type
900  SELECT CASE(interfacecondition%integrationType)
901  CASE(interface_condition_gauss_integration)
902  !Pointers to interface variables (columns of interface element matrix)
903  interfaceinterpolation=>interfaceequations%INTERPOLATION%INTERFACE_INTERPOLATION
904  interfacegeometricfield=>interfaceinterpolation%GEOMETRIC_FIELD
905  interfacedependentfield=>interfaceinterpolation%DEPENDENT_FIELD
906  interfacegeometricbasis=>interfacegeometricfield%DECOMPOSITION%DOMAIN(interfacegeometricfield% &
907  & decomposition%MESH_COMPONENT_NUMBER)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(elementnumber)%BASIS
908  interfacedependentbasis=>interfacedependentfield%DECOMPOSITION%DOMAIN(interfacedependentfield% &
909  & decomposition%MESH_COMPONENT_NUMBER)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(elementnumber)%BASIS
910  !Integrate using the interface quadrature scheme
911  interfacequadraturescheme=>interfacegeometricbasis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
912  lagrangevariable=>interfaceequations%INTERFACE_MAPPING%LAGRANGE_VARIABLE
913  !lagrangeVariableNumberOfComponents=>interfaceEquations%INTERFACE_MAPPING%LAGRANGE_VARIABLE%NUMBER_OF_COMPONENTS
914  lagrangevariabletype=lagrangevariable%VARIABLE_TYPE
915  !Get element interpolation parameters from the first geometric interpolation set (to get Jacobian for interface surface integral)
916  CALL field_interpolation_parameters_element_get(field_values_set_type,elementnumber,interfaceinterpolation% &
917  & geometric_interpolation(1)%INTERPOLATION_PARAMETERS(field_u_variable_type)%PTR,err,error,*999)
918  !\todo Use matrix coefficients in solver routines when assembling solver matrices instead of multiplying them here
919  matrixcoefficient=1.0_dp
920  !Loop over interface matrices (1st solid, 2nd fluid)
921  DO coupledmeshidx=1,interfaceequations%INTERFACE_MATRICES%NUMBER_OF_INTERFACE_MATRICES
922  IF(interfaceequations%INTERFACE_MATRICES%MATRICES(coupledmeshidx)%PTR%UPDATE_MATRIX) THEN
923  !\todo Use matrix coefficients in solver routines when assembling solver matrices instead of multiplying them here
924  IF(coupledmeshidx>1) THEN
925  matrixcoefficient=-1.0_dp
926  ENDIF
927  !Pointers to the coupledMeshIdx'th coupled mesh variables (rows of interface element matrix)
928  coupledmeshdependentfield=>interfaceequations%INTERPOLATION%VARIABLE_INTERPOLATION(coupledmeshidx)%DEPENDENT_FIELD
929  elementconnectivity=>interfacecondition%INTERFACE%MESH_CONNECTIVITY%ELEMENT_CONNECTIVITY(elementnumber,coupledmeshidx)
930  coupledmeshelementnumber=elementconnectivity%COUPLED_MESH_ELEMENT_NUMBER
931  interfacematrixvariable=> &
932  & interfaceequations%INTERFACE_MAPPING%INTERFACE_MATRIX_ROWS_TO_VAR_MAPS(coupledmeshidx)%VARIABLE
933  coupledmeshvariabletype=interfacematrixvariable%VARIABLE_TYPE
934  interfaceelementmatrix=>interfaceequations%INTERFACE_MATRICES%MATRICES(coupledmeshidx)%PTR%ELEMENT_MATRIX
935  interfaceconnectivitybasis=>interfacecondition%INTERFACE%MESH_CONNECTIVITY%BASIS
936 
937  !coupledMeshDependentInterpolation=>interfaceEquations%INTERPOLATION%VARIABLE_INTERPOLATION(coupledMeshIdx)% &
938  ! & DEPENDENT_INTERPOLATION
939 
940  !=======================================================================================================================
941  !Loop over gauss points
942  DO gausspoint=1,interfacequadraturescheme%NUMBER_OF_GAUSS
943  !CALL FIELD_INTERPOLATE_GAUSS(FIRST_PART_DERIV,BASIS_DEFAULT_QUADRATURE_SCHEME,GaussPoint, &
944  ! & coupledMeshDependentInterpolation%GEOMETRIC_INTERPOLATION(1)%INTERPOLATED_POINT(FIELD_U_VARIABLE_TYPE)%PTR, &
945  ! & err,error,*999)
946  !CALL FIELD_INTERPOLATED_POINT_METRICS_CALCULATE(interfaceGeometricBasis%NUMBER_OF_XI,interfaceInterpolation% &
947  ! & GEOMETRIC_INTERPOLATION(1)%INTERPOLATED_POINT_METRICS(FIELD_U_VARIABLE_TYPE)%PTR,err,error,*999)
948  !=====================================================================================================================
949  !Interpolates field at given gauss point, includes first partial derivatives
950  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,gausspoint,interfaceinterpolation% &
951  & geometric_interpolation(1)%INTERPOLATED_POINT(field_u_variable_type)%PTR,err,error,*999)
952  !Calculates the interpolated point metrics and the associated interpolated point
953  CALL field_interpolated_point_metrics_calculate(interfacegeometricbasis%NUMBER_OF_XI,interfaceinterpolation% &
954  & geometric_interpolation(1)%INTERPOLATED_POINT_METRICS(field_u_variable_type)%PTR,err,error,*999)
955  !=====================================================================================================================
956  ! R W G = GAUSSWEIGTHS * JACOBIAN
957  rwg=interfaceinterpolation%GEOMETRIC_INTERPOLATION(1)%INTERPOLATED_POINT_METRICS(field_u_variable_type)%PTR% &
958  & jacobian*interfacequadraturescheme%GAUSS_WEIGHTS(gausspoint)
959  IF(interfacecondition%METHOD==interface_condition_penalty_method .AND. &
960  & coupledmeshidx==interfaceequations%INTERFACE_MATRICES%NUMBER_OF_INTERFACE_MATRICES) THEN
961  CALL flagerror("Not implemented.",err,error,*999)
962  ELSE
963  !===================================================================================================================
964  !\todo defaults to first mesh component, generalise
965  !TODO Originally XI=...
966  xi(1:SIZE(elementconnectivity%XI,1))=interfaceoperators_interftocoupledmeshgausstransform( &
967  & elementconnectivity,interfaceconnectivitybasis,gausspoint,err,error)
968  ! Loop over number of Lagrange variable components as not all components in the dependent field variable may be coupled
969  !\todo Currently Lagrange field variable component numbers must match each coupled dependent field variable component numbers. Generalise ordering
970  DO rowcomponentidx=1,lagrangevariable%NUMBER_OF_COMPONENTS
971  rowmeshcomponentnumber=interfacematrixvariable%COMPONENTS(rowcomponentidx)%MESH_COMPONENT_NUMBER
972  coupledmeshbasis=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(rowmeshcomponentnumber)%PTR%TOPOLOGY% &
973  & elements%ELEMENTS(coupledmeshelementnumber)%BASIS
974 
975  SELECT CASE(interfacedependentbasis%NUMBER_OF_XI)
976 
977  CASE(1) !1D interface (line)
978  connectedline=elementconnectivity%CONNECTED_LINE
979  decompositionlinenumber=coupledmeshdependentfield%DECOMPOSITION%TOPOLOGY% &
980  & elements%ELEMENTS(coupledmeshelementnumber)%ELEMENT_LINES(connectedline)
981  coupledmeshdomainline=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(rowmeshcomponentnumber)%PTR%TOPOLOGY% &
982  & lines%LINES(decompositionlinenumber)
983  DO locallinenodeidx=1,coupledmeshbasis%NUMBER_OF_NODES_IN_LOCAL_LINE(connectedline)
984  localelementnode=coupledmeshbasis%NODE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx,connectedline)
985  DO derivativeidx=1,coupledmeshdomainline%BASIS%NUMBER_OF_DERIVATIVES(locallinenodeidx)
986  !???????????????????
987  derivative=coupledmeshbasis%DERIVATIVE_NUMBERS_IN_LOCAL_LINE(locallinenodeidx,connectedline)
988  derivative=coupledmeshdomainline%DERIVATIVES_IN_LINE(1,derivativeidx,locallinenodeidx)
989  !???????????????????
990  rowparameteridx=coupledmeshbasis%ELEMENT_PARAMETER_INDEX(derivative,localelementnode)
991  !===========================================================================================================
992  ! P G M S I - this represents the D E P E N D E N T _ F I E L D S (solid, fluid)
993  !Evaluates the appropriate partial derivative index at position XI for the basis
994  pgmsi=basis_evaluate_xi(coupledmeshbasis,rowparameteridx,no_part_deriv,xi,err,error)
995  rowidx=rowparameteridx+coupledmeshbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
996  DO interfacenode=1,interfacedependentbasis%NUMBER_OF_NODES
997  DO interfacederivative=1,interfacedependentbasis%NUMBER_OF_DERIVATIVES(interfacenode)
998  !\todo requires equal number of nodes between interface mesh and coupled mesh. Generalize
999  colparameteridx=interfacedependentbasis%ELEMENT_PARAMETER_INDEX(interfacederivative,interfacenode)
1000  !=======================================================================================================
1001  ! P G N S I - this represents the L A M B D A
1002  pgnsi=interfacequadraturescheme%GAUSS_BASIS_FNS(colparameteridx,no_part_deriv,gausspoint)
1003  colidx=colparameteridx+interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
1004  !\todo Use matrix coefficients in solver routines when assembling solver matrices instead of multiplying them here
1005  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx)+ &
1006  & pgnsi*pgmsi*rwg*matrixcoefficient
1007  ENDDO !interfaceDerivative
1008  ENDDO !interfaceNode
1009  ENDDO !derivativeIdx
1010  ENDDO !localLineNodeIdx
1011 
1012  CASE(2) !2D interface (face)
1013 
1014  SELECT CASE(coupledmeshbasis%NUMBER_OF_XI)
1015 
1016  CASE(2) !Coupled Mesh has 2 xi directions
1017  DO localelementnode=1,coupledmeshbasis%NUMBER_OF_NODES
1018  DO derivative=1,coupledmeshbasis%NUMBER_OF_DERIVATIVES(localelementnode)
1019  rowparameteridx=coupledmeshbasis%ELEMENT_PARAMETER_INDEX(derivative,localelementnode)
1020  !=========================================================================================================
1021  ! P G M S I
1022  pgmsi=basis_evaluate_xi(coupledmeshbasis,rowparameteridx,no_part_deriv, &
1023  & xi(1:coupledmeshbasis%NUMBER_OF_XI),err,error)
1024  rowidx=rowparameteridx+coupledmeshbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
1025  DO interfacenode=1,interfacedependentbasis%NUMBER_OF_NODES
1026  DO interfacederivative=1,interfacedependentbasis%NUMBER_OF_DERIVATIVES(interfacenode)
1027  !\todo requires equal number of nodes between interface mesh and coupled mesh. Generalize
1028  colparameteridx=interfacedependentbasis%ELEMENT_PARAMETER_INDEX(interfacederivative,interfacenode)
1029  !=====================================================================================================
1030  ! P G N S I
1031  pgnsi=interfacequadraturescheme%GAUSS_BASIS_FNS(colparameteridx,no_part_deriv,gausspoint)
1032  colidx=colparameteridx+interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
1033  !\todo Use matrix coefficients in solver routines when assembling solver matrices instead of multiplying them here
1034  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx)+ &
1035  & pgnsi*pgmsi*rwg*matrixcoefficient
1036  ENDDO !interfaceDerivative
1037  ENDDO !interfaceNode
1038  ENDDO !derivative
1039  ENDDO !localElementNode
1040 
1041  CASE(3) !Coupled Mesh has 3 xi directions
1042  connectedface = elementconnectivity%CONNECTED_FACE
1043  decompositionfacenumber=coupledmeshdependentfield%DECOMPOSITION%TOPOLOGY% &
1044  & elements%ELEMENTS(coupledmeshelementnumber)%ELEMENT_FACES(connectedface)
1045  coupledmeshdomainface=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(rowmeshcomponentnumber)%PTR%TOPOLOGY% &
1046  & faces%FACES(decompositionfacenumber)
1047  DO localfacenodeidx=1,coupledmeshbasis%NUMBER_OF_NODES_IN_LOCAL_FACE(connectedface)
1048  localelementnode=coupledmeshbasis%NODE_NUMBERS_IN_LOCAL_FACE(localfacenodeidx,connectedface)
1049  DO derivativeidx=1,coupledmeshdomainface%BASIS%NUMBER_OF_DERIVATIVES(localfacenodeidx)
1050  derivative=coupledmeshbasis% &
1051  & derivative_numbers_in_local_face(derivativeidx,localfacenodeidx,connectedface)
1052  rowparameteridx=coupledmeshbasis%ELEMENT_PARAMETER_INDEX(derivative,localelementnode)
1053  !=========================================================================================================
1054  ! P G M S I
1055  pgmsi=basis_evaluate_xi(coupledmeshbasis,rowparameteridx,no_part_deriv, &
1056  & xi(1:coupledmeshbasis%NUMBER_OF_XI),err,error)
1057  rowidx=rowparameteridx+coupledmeshbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
1058  DO interfacenode=1,interfacedependentbasis%NUMBER_OF_NODES
1059  DO interfacederivative=1,interfacedependentbasis%NUMBER_OF_DERIVATIVES(interfacenode)
1060  !\todo requires equal number of nodes between interface mesh and coupled mesh. Generalize
1061  colparameteridx=interfacedependentbasis%ELEMENT_PARAMETER_INDEX(interfacederivative,interfacenode)
1062  !=====================================================================================================
1063  ! P G N S I
1064  pgnsi=interfacequadraturescheme%GAUSS_BASIS_FNS(colparameteridx,no_part_deriv,gausspoint)
1065  colidx=colparameteridx+interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS*(rowcomponentidx-1)
1066  !\todo Use matrix coefficients in solver routines when assembling solver matrices instead of multiplying them here
1067  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx)+ &
1068  & pgnsi*pgmsi*rwg*matrixcoefficient
1069  ENDDO !interfaceDerivative
1070  ENDDO !interfaceNode
1071  ENDDO !derivativeIdx
1072  ENDDO !FaceNodeIdx
1073 
1074  END SELECT !coupledMeshBasis%NUMBER_OF_XI
1075 
1076  END SELECT !interfaceDependentBasis%NUMBER_OF_XI
1077 
1078  ENDDO !rowComponentIdx
1079  ENDIF
1080  ENDDO !GaussPoint
1081 
1082  !Scale factor adjustment
1083  !\todo check if scale factor adjustments are already made elsewhere eg when calculating the interface matrix contribution to the residual for non-linear problems
1084  !\todo update looping of variables/components for non-zero matrix elements as done above
1085  IF(interfacecondition%METHOD==interface_condition_penalty_method .AND. &
1086  & coupledmeshidx==interfaceequations%INTERFACE_MATRICES%NUMBER_OF_INTERFACE_MATRICES) THEN
1087  CALL flagerror("Not implemented.",err,error,*999)
1088  ELSE
1089  !Scale factor adjustment for the Lagrange Variable (columns)
1090  IF(interfacedependentfield%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
1091  CALL field_interpolationparametersscalefactorselementget(elementnumber, &
1092  & interfaceinterpolation%DEPENDENT_INTERPOLATION(1)%INTERPOLATION_PARAMETERS(lagrangevariabletype)%PTR, &
1093  & err,error,*999)
1094  rowidx=0
1095  !Use Lagrange variable number of components here since we are only dealing with Lagrange variable scale factors
1096  !\todo Currently Lagrange field variable component numbers must match each coupled dependent field variable component numbers. Generalise ordering
1097  DO rowcomponentidx=1,lagrangevariable%NUMBER_OF_COMPONENTS
1098  rowmeshcomponentnumber=interfacematrixvariable%COMPONENTS(rowcomponentidx)%MESH_COMPONENT_NUMBER
1099  coupledmeshbasis=>coupledmeshdependentfield%DECOMPOSITION%DOMAIN(rowmeshcomponentnumber)%PTR%TOPOLOGY% &
1100  & elements%ELEMENTS(coupledmeshelementnumber)%BASIS
1101  !Loop over element rows
1102  DO rowparameteridx=1,coupledmeshbasis%NUMBER_OF_ELEMENT_PARAMETERS
1103  rowidx=rowidx+1
1104  colidx=0
1105  !Loop over element columns
1106  DO colcomponentidx=1,lagrangevariable%NUMBER_OF_COMPONENTS
1107  DO colparameteridx=1,interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS
1108  colidx=colidx+1
1109  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx) * &
1110  & interfaceinterpolation%DEPENDENT_INTERPOLATION(1)% &
1111  & interpolation_parameters(lagrangevariabletype)%PTR%SCALE_FACTORS(colparameteridx,colcomponentidx)
1112  ENDDO !colParameterIdx
1113  ENDDO !colComponentIdx
1114  ENDDO !rowParameterIdx
1115  ENDDO !rowComponentIdx
1116  ENDIF
1117  !Scale factor adjustment for the row dependent variable
1118  IF(coupledmeshdependentfield%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
1119  CALL field_interpolationparametersscalefactorselementget(coupledmeshelementnumber, &
1120  & interfaceequations%INTERPOLATION%VARIABLE_INTERPOLATION(coupledmeshidx)% &
1121  & dependent_interpolation(1)%INTERPOLATION_PARAMETERS(coupledmeshvariabletype)%PTR,err,error,*999)
1122  rowidx=0
1123  DO rowcomponentidx=1,interfacematrixvariable%NUMBER_OF_COMPONENTS
1124  !Loop over element rows
1125  DO rowparameteridx=1,coupledmeshbasis%NUMBER_OF_ELEMENT_PARAMETERS
1126  rowidx=rowidx+1
1127  colidx=0
1128  !Loop over element columns
1129  DO colcomponentidx=1,lagrangevariable%NUMBER_OF_COMPONENTS
1130  DO colparameteridx=1,interfacedependentbasis%NUMBER_OF_ELEMENT_PARAMETERS
1131  colidx=colidx+1
1132  interfaceelementmatrix%MATRIX(rowidx,colidx)=interfaceelementmatrix%MATRIX(rowidx,colidx)* &
1133  & interfaceequations%INTERPOLATION%VARIABLE_INTERPOLATION(coupledmeshidx)% &
1134  & dependent_interpolation(1)%INTERPOLATION_PARAMETERS(coupledmeshvariabletype)%PTR% &
1135  & scale_factors(rowparameteridx,rowcomponentidx)
1136  ENDDO !colParameterIdx
1137  ENDDO !colComponentIdx
1138  ENDDO !rowParameterIdx
1139  ENDDO !rowComponentIdx
1140  ENDIF
1141  ENDIF
1142  ENDIF
1143  ENDDO ! coupledMeshIdx
1144 
1145  CASE(interface_condition_data_points_integration)
1146  CALL flagerror("Not implemented.",err,error,*999)
1147  CASE DEFAULT
1148  localerror="Interface condition integration type "//trim(number_to_vstring(interfacecondition%integrationType, &
1149  & "*",err,error))// " is not valid."
1150  CALL flagerror(localerror,err,error,*999)
1151  END SELECT !interfaceCondition%integrationType
1152 
1153  CASE(interface_condition_augmented_lagrange_method)
1154  CALL flagerror("Not implemented.",err,error,*999)
1155  CASE DEFAULT
1156  localerror="Interface condition method "//trim(number_to_vstring(interfacecondition%METHOD,"*",err,error))// &
1157  & " is not valid."
1158  CALL flagerror(localerror,err,error,*999)
1159  END SELECT
1160 
1161  exits("SolidFluidOperator_FiniteElementCalculate")
1162  RETURN
1163 999 errorsexits("SolidFluidOperator_FiniteElementCalculate",err,error)
1164  RETURN 1
1165 
1166  END SUBROUTINE solidfluidoperator_finiteelementcalculate
1167 
1168  !
1169  !================================================================================================================================
1170  !
1171 
1172  FUNCTION interfaceoperators_interftocoupledmeshgausstransform(elementConnectivity,interfaceConnectivityBasis,GaussPoint,err,error)
1173 
1174  !Argument variables
1175  TYPE(interface_element_connectivity_type), POINTER :: elementconnectivity
1176  TYPE(basis_type), POINTER :: interfaceconnectivitybasis
1177  INTEGER(INTG), INTENT(IN) :: gausspoint
1178  INTEGER(INTG), INTENT(OUT) :: err
1179  TYPE(varying_string), INTENT(OUT) :: error
1180  !Function variable
1181  REAL(DP) :: interfaceoperators_interftocoupledmeshgausstransform(size(elementconnectivity%xi,1))
1182  !Local Variables
1183  INTEGER(INTG) :: rowparameteridx
1184 
1185  enters("InterfaceOperators_InterfToCoupledMeshGaussTransform",err,error,*999)
1186 
1187  interfaceoperators_interftocoupledmeshgausstransform=0.0_dp
1188  DO rowparameteridx = 1,interfaceconnectivitybasis%NUMBER_OF_ELEMENT_PARAMETERS
1189  interfaceoperators_interftocoupledmeshgausstransform(:)= interfaceoperators_interftocoupledmeshgausstransform(:) + &
1190  & interfaceconnectivitybasis%QUADRATURE%QUADRATURE_SCHEME_MAP &
1191  & (basis_default_quadrature_scheme)%PTR%GAUSS_BASIS_FNS(rowparameteridx,no_part_deriv,gausspoint) * &
1192  & elementconnectivity%XI(:,1,rowparameteridx)
1193  ENDDO
1194 
1195  exits("InterfaceOperators_InterfToCoupledMeshGaussTransform")
1196  RETURN
1197 999 errors("InterfaceOperators_InterfToCoupledMeshGaussTransform",err,error)
1198  exits("InterfaceOperators_InterfToCoupledMeshGaussTransform")
1199  RETURN
1200 
1201  END FUNCTION interfaceoperators_interftocoupledmeshgausstransform
1202 
1203 END MODULE interface_operators_routines
basis_routines
This module contains all basis function routines.
Definition: basis_routines.f90:45
base_routines::enters
subroutine, public enters(NAME, ERR, ERROR,)
Records the entry into the named procedure and initialises the error code.
Definition: base_routines.f90:363
kinds::ptr
integer, parameter ptr
Pointer integer kind.
Definition: kinds.f90:58
types::domain_face_type
Contains the information for a face in a domain.
Definition: types.f90:644
strings::number_to_vstring
Converts a number to its equivalent varying string representation.
Definition: strings.f90:161
basis_routines::basis_evaluate_xi
Evaluates the appropriate partial derivative index for the specificied basis function at a Xi locatio...
Definition: basis_routines.f90:168
types::interface_matrix_type
Contains information about an interface matrix.
Definition: types.f90:1978
constants::no_part_deriv
integer(intg), parameter no_part_deriv
No partial derivative i.e., u.
Definition: constants.f90:177
interface_conditions_constants::interface_condition_lagrange_multipliers_method
integer(intg), parameter interface_condition_lagrange_multipliers_method
Lagrange multipliers interface condition method.
Definition: interface_conditions_constants.f90:57
types::interface_penalty_type
Contains information about the penalty field information for an interface condition.
Definition: types.f90:2129
strings
This module contains all string manipulation and transformation routines.
Definition: strings.f90:45
types::interface_condition_type
Contains information for the interface condition data.
Definition: types.f90:2155
constants::first_part_deriv
integer(intg), parameter first_part_deriv
First partial derivative i.e., du/ds.
Definition: constants.f90:178
timer
This module contains routines for timing the program.
Definition: timer_f.f90:45
types::decompositionelementdatapointstype
Contains information on the projected data points on an element, for decomposition since data points ...
Definition: types.f90:1034
types::field_type
Contains information for a field defined on a region.
Definition: types.f90:1346
iso_varying_string
This module provides an iso_varying_string module, conformant to the API specified in ISO/IEC 1539-2:...
Definition: iso_varying_string.f90:62
interface_conditions_constants::interface_condition_augmented_lagrange_method
integer(intg), parameter interface_condition_augmented_lagrange_method
Augmented Lagrange multiplers interface condition method.
Definition: interface_conditions_constants.f90:58
basis_routines::basis_default_quadrature_scheme
integer(intg), parameter, public basis_default_quadrature_scheme
Identifier for the default quadrature scheme.
Definition: basis_routines.f90:125
iso_varying_string::varying_string
Definition: iso_varying_string.f90:75
interface_conditions_constants::interface_condition_gauss_integration
integer(intg), parameter interface_condition_gauss_integration
Gauss points integration type, i.e. Loop over element Gauss points and sum up their contribution...
Definition: interface_conditions_constants.f90:101
interface_mapping_routines
This module contains all interface mapping routines.
Definition: interface_mapping_routines.f90:45
constants
This module contains all program wide constants.
Definition: constants.f90:45
types::domain_line_type
Contains the information for a line in a domain.
Definition: types.f90:622
base_routines::exits
subroutine, public exits(NAME)
Records the exit out of the named procedure.
Definition: base_routines.f90:504
types::field_interpolated_point_metrics_ptr_type
Definition: types.f90:1124
types
This module contains all type definitions in order to avoid cyclic module references.
Definition: types.f90:70
interface_conditions_constants::interface_condition_data_points_integration
integer(intg), parameter interface_condition_data_points_integration
Data points integration type i.e. Loop over data points and sum up their contribution.
Definition: interface_conditions_constants.f90:102
base_routines
This module contains all the low-level base routines e.g., all debug, control, and low-level communic...
Definition: base_routines.f90:45
iso_varying_string::trim
Definition: iso_varying_string.f90:204
interface_conditions_constants::interface_condition_penalty_method
integer(intg), parameter interface_condition_penalty_method
Penalty interface condition method.
Definition: interface_conditions_constants.f90:59
types::element_matrix_type
Contains information for an element matrix.
Definition: types.f90:1387
types::field_interpolated_point_ptr_type
Definition: types.f90:1136
types::interface_equations_domain_interpolation_type
Contains information about the interpolation for a domain (interface or coupled mesh) in the interfac...
Definition: types.f90:2089
interface_conditions_constants
This module defines all constants shared across interface condition routines.
Definition: interface_conditions_constants.f90:45
interface_conditions_constants::interface_condition_point_to_point_method
integer(intg), parameter interface_condition_point_to_point_method
Point to point interface condition method.
Definition: interface_conditions_constants.f90:60
types::field_interpolated_point_type
Contains the interpolated value (and the derivatives wrt xi) of a field at a point. Old CMISS name XG.
Definition: types.f90:1129
types::quadrature_scheme_type
Contains information for a particular quadrature scheme.
Definition: types.f90:141
interface_equations_routines
This module handles all interface equations routines.
Definition: interface_equations_routines.f90:45
matrix_vector
This module contains all routines dealing with (non-distributed) matrix and vectors types...
Definition: matrix_vector.f90:129
types::field_variable_type
Contains information for a field variable defined on a field.
Definition: types.f90:1289
types::interfacepointsconnectivitytype
Contains information on the data point coupling/points connectivity between meshes in the an interfac...
Definition: types.f90:2218
types::interface_element_connectivity_type
Contains information on the mesh connectivity for a given coupled mesh element.
Definition: types.f90:2185
base_routines::errors
subroutine, public errors(NAME, ERR, ERROR)
Records the exiting error of the subroutine.
Definition: base_routines.f90:481
types::interface_type
Contains information for the interface data.
Definition: types.f90:2228
types::basis_type
Contains all information about a basis .
Definition: types.f90:184
base_routines::flag_error
Flags an error condition.
Definition: base_routines.f90:226
base_routines::flagerror
Flags an error condition.
Definition: base_routines.f90:240
types::field_interpolation_parameters_ptr_type
Definition: types.f90:1151
constants::zero_tolerance
real(dp), parameter zero_tolerance
Definition: constants.f90:70
kinds
This module contains all kind definitions.
Definition: kinds.f90:45
types::interface_equations_type
Contains information about the interface equations for an interface condition.
Definition: types.f90:2110
input_output
This module handles all formating and input and output.
Definition: input_output.f90:45