OpenCMISS-Iron Internal API Documentation
bioelectric_finite_elasticity_routines.f90
Go to the documentation of this file.
1 
43 
45 
46 
47 MODULE bioelectric_finite_elasticity_routines
48 
49  USE base_routines
50  USE basis_routines
51  USE bioelectric_routines
52  USE biodomain_equation_routines
53  USE constants
54  USE control_loop_routines
55  USE equations_routines
56  USE equations_mapping_routines
57  USE equations_matrices_routines
58  USE equations_set_constants
59  USE field_io_routines
60  USE field_routines
61  USE finite_elasticity_routines
62  USE input_output
63  USE iso_varying_string
64  USE kinds
65  USE maths
66  USE problem_constants
67  USE strings
68  USE solver_routines
69  USE types
70 
71 #include "macros.h"
72 
73  IMPLICIT NONE
74 
75  PUBLIC bioelectricfiniteelasticity_equationssetsetup
76  PUBLIC bioelectricfiniteelasticity_equationssetsolutionmethodset
77 
78  PUBLIC bioelectric_finite_elasticity_problem_setup
79  PUBLIC bioelectricfiniteelasticity_problemspecificationset
80 
81  PUBLIC bioelectricfiniteelasticity_finiteelementcalculate
82 
83  PUBLIC bioelectric_finite_elasticity_pre_solve
84  PUBLIC bioelectric_finite_elasticity_post_solve
85 
86  PUBLIC bioelectricfiniteelasticity_controllooppreloop
87  PUBLIC bioelectricfiniteelasticity_controllooppostloop
88 
89  PUBLIC bioelectricfiniteelasticity_updategeometricfield
90 
91 CONTAINS
92 
93  !
94  !================================================================================================================================
95  !
96 
98  SUBROUTINE bioelectricfiniteelasticity_equationssetsolutionmethodset(EQUATIONS_SET,SOLUTION_METHOD,ERR,ERROR,*)
99 
100  !Argument variables
101  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
102  INTEGER(INTG), INTENT(IN) :: SOLUTION_METHOD
103  INTEGER(INTG), INTENT(OUT) :: ERR
104  TYPE(varying_string), INTENT(OUT) :: ERROR
105  !Local Variables
106  TYPE(varying_string) :: LOCAL_ERROR
107 
108  enters("BioelectricFiniteElasticity_EquationsSetSolutionMethodSet",err,error,*999)
109 
110  IF(ASSOCIATED(equations_set)) THEN
111  IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
112  CALL flagerror("Equations set specification is not allocated.",err,error,*999)
113  ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
114  CALL flagerror("Equations set specification must have three entries for a "// &
115  & "Bioelectric-finite elasticity type equations set.",err,error,*999)
116  END IF
117  SELECT CASE(equations_set%SPECIFICATION(3))
118  CASE(equations_set_standard_monodomain_elasticity_subtype, &
119  & equations_set_1d3d_monodomain_elasticity_subtype, &
120  & equations_set_monodomain_elasticity_w_titin_subtype, &
121  & equations_set_monodomain_elasticity_velocity_subtype, &
122  & equations_set_1d3d_monodomain_active_strain_subtype)
123  SELECT CASE(solution_method)
124  CASE(equations_set_fem_solution_method)
125  equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
126  CASE(equations_set_bem_solution_method)
127  CALL flagerror("Not implemented.",err,error,*999)
128  CASE(equations_set_fd_solution_method)
129  CALL flagerror("Not implemented.",err,error,*999)
130  CASE(equations_set_fv_solution_method)
131  CALL flagerror("Not implemented.",err,error,*999)
132  CASE(equations_set_gfem_solution_method)
133  CALL flagerror("Not implemented.",err,error,*999)
134  CASE(equations_set_gfv_solution_method)
135  CALL flagerror("Not implemented.",err,error,*999)
136  CASE DEFAULT
137  local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
138  CALL flagerror(local_error,err,error,*999)
139  END SELECT
140  CASE DEFAULT
141  local_error="Equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
142  & " is not valid for a bioelectrics finite elasticity equation type of a multi physics equations set class."
143  CALL flagerror(local_error,err,error,*999)
144  END SELECT
145  ELSE
146  CALL flagerror("Equations set is not associated.",err,error,*999)
147  ENDIF
148 
149  exits("BioelectricFiniteElasticity_EquationsSetSolutionMethodSet")
150  RETURN
151 999 errors("BioelectricFiniteElasticity_EquationsSetSolutionMethodSet",err,error)
152  exits("BioelectricFiniteElasticity_EquationsSetSolutionMethodSet")
153  RETURN 1
154 
155  END SUBROUTINE bioelectricfiniteelasticity_equationssetsolutionmethodset
156 
157  !
158  !================================================================================================================================
159  !
160 
162  SUBROUTINE bioelectricfiniteelasticity_equationssetsetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)
163 
164  !Argument variables
165  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
166  TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
167  INTEGER(INTG), INTENT(OUT) :: ERR
168  TYPE(varying_string), INTENT(OUT) :: ERROR
169 
170  enters("BioelectricFiniteElasticity_EquationsSetSetup",err,error,*999)
171 
172  CALL flagerror("BioelectricFiniteElasticity_EquationsSetSetup is not implemented.",err,error,*999)
173 
174  exits("BioelectricFiniteElasticity_EquationsSetSetup")
175  RETURN
176 999 errors("BioelectricFiniteElasticity_EquationsSetSetup",err,error)
177  exits("BioelectricFiniteElasticity_EquationsSetSetup")
178  RETURN 1
179 
180  END SUBROUTINE bioelectricfiniteelasticity_equationssetsetup
181 
182  !
183  !================================================================================================================================
184  !
185 
187  SUBROUTINE bioelectricfiniteelasticity_finiteelementcalculate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)
188 
189  !Argument variables
190  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
191  INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
192  INTEGER(INTG), INTENT(OUT) :: ERR
193  TYPE(varying_string), INTENT(OUT) :: ERROR
194 
195  enters("BioelectricFiniteElasticity_FiniteElementCalculate",err,error,*999)
196 
197  CALL flagerror("BioelectricFiniteElasticity_FiniteElementCalculate is not implemented.",err,error,*999)
198 
199  exits("BioelectricFiniteElasticity_FiniteElementCalculate")
200  RETURN
201 999 errors("BioelectricFiniteElasticity_FiniteElementCalculate",err,error)
202  exits("BioelectricFiniteElasticity_FiniteElementCalculate")
203  RETURN 1
204 
205  END SUBROUTINE bioelectricfiniteelasticity_finiteelementcalculate
206 
207  !
208  !================================================================================================================================
209  !
210 
212  SUBROUTINE bioelectricfiniteelasticity_problemspecificationset(problem,problemSpecification,err,error,*)
213 
214  !Argument variables
215  TYPE(problem_type), POINTER :: problem
216  INTEGER(INTG), INTENT(IN) :: problemSpecification(:)
217  INTEGER(INTG), INTENT(OUT) :: err
218  TYPE(varying_string), INTENT(OUT) :: error
219  !Local Variables
220  TYPE(varying_string) :: localError
221  INTEGER(INTG) :: problemSubtype
222 
223  enters("BioelectricFiniteElasticity_ProblemSpecificationSet",err,error,*999)
224 
225  IF(ASSOCIATED(problem)) THEN
226  IF(SIZE(problemspecification,1)==3) THEN
227  problemsubtype=problemspecification(3)
228  SELECT CASE(problemsubtype)
229  CASE(problem_gudunov_monodomain_simple_elasticity_subtype, &
230  & problem_gudunov_monodomain_1d3d_elasticity_subtype, &
231  & problem_monodomain_elasticity_w_titin_subtype, &
232  & problem_monodomain_elasticity_velocity_subtype, &
233  & problem_monodomain_1d3d_active_strain_subtype)
234  !ok
235  CASE DEFAULT
236  localerror="The third problem specification of "//trim(numbertovstring(problemsubtype,"*",err,error))// &
237  & " is not valid for a bioelectric finite elasticity type of a multi physics problem."
238  CALL flagerror(localerror,err,error,*999)
239  END SELECT
240  IF(ALLOCATED(problem%specification)) THEN
241  CALL flagerror("Problem specification is already allocated.",err,error,*999)
242  ELSE
243  ALLOCATE(problem%specification(3),stat=err)
244  IF(err/=0) CALL flagerror("Could not allocate problem specification.",err,error,*999)
245  END IF
246  problem%specification(1:3)=[problem_multi_physics_class,problem_bioelectric_finite_elasticity_type, &
247  & problemsubtype]
248  ELSE
249  CALL flagerror("Bioelectric finite elasticity problem specification must have three entries.",err,error,*999)
250  END IF
251  ELSE
252  CALL flagerror("Problem is not associated.",err,error,*999)
253  END IF
254 
255  exits("BioelectricFiniteElasticity_ProblemSpecificationSet")
256  RETURN
257 999 errors("BioelectricFiniteElasticity_ProblemSpecificationSet",err,error)
258  exits("BioelectricFiniteElasticity_ProblemSpecificationSet")
259  RETURN 1
260 
261  END SUBROUTINE bioelectricfiniteelasticity_problemspecificationset
262 
263  !
264  !================================================================================================================================
265  !
266 
268  SUBROUTINE bioelectric_finite_elasticity_problem_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)
269 
270  !Argument variables
271  TYPE(problem_type), POINTER :: PROBLEM
272  TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
273  INTEGER(INTG), INTENT(OUT) :: ERR
274  TYPE(varying_string), INTENT(OUT) :: ERROR
275  !Local Variables
276  TYPE(control_loop_type), POINTER :: CONTROL_LOOP,CONTROL_LOOP_ROOT
277  TYPE(control_loop_type), POINTER :: MONODOMAIN_SUB_LOOP,ELASTICITY_SUB_LOOP
278  TYPE(solver_type), POINTER :: SOLVER
279  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
280  TYPE(cellml_equations_type), POINTER :: CELLML_EQUATIONS
281  TYPE(solvers_type), POINTER :: SOLVERS,MONODOMAIN_SOLVERS,ELASTICITY_SOLVERS
282  TYPE(varying_string) :: LOCAL_ERROR
283 
284  enters("BIOELECTRIC_FINITE_ELASTICITY_PROBLEM_SETUP",err,error,*999)
285 
286  NULLIFY(control_loop)
287  NULLIFY(monodomain_sub_loop)
288  NULLIFY(elasticity_sub_loop)
289  NULLIFY(solver)
290  NULLIFY(solvers)
291  NULLIFY(monodomain_solvers)
292  NULLIFY(elasticity_solvers)
293  NULLIFY(solver_equations)
294  NULLIFY(cellml_equations)
295 
296  IF(ASSOCIATED(problem)) THEN
297  IF(.NOT.ALLOCATED(problem%specification)) THEN
298  CALL flagerror("Problem specification is not allocated.",err,error,*999)
299  ELSE IF(SIZE(problem%specification,1)<3) THEN
300  CALL flagerror("Problem specification must have three entries for a bioelectric-finite elasticity problem.",err,error,*999)
301  END IF
302  SELECT CASE(problem%SPECIFICATION(3))
303 
304  !--------------------------------------------------------------------
305  ! Transient Gudunov monodomain, simple finite elasticity
306  !--------------------------------------------------------------------
307  CASE(problem_gudunov_monodomain_simple_elasticity_subtype,problem_gudunov_monodomain_1d3d_elasticity_subtype, &
308  & problem_monodomain_elasticity_w_titin_subtype,problem_monodomain_elasticity_velocity_subtype, &
309  & problem_monodomain_1d3d_active_strain_subtype)
310  SELECT CASE(problem_setup%SETUP_TYPE)
311  CASE(problem_setup_initial_type)
312  SELECT CASE(problem_setup%ACTION_TYPE)
313  CASE(problem_setup_start_action)
314  !Do nothing
315  CASE(problem_setup_finish_action)
316  !Do nothing
317  CASE DEFAULT
318  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
319  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
320  & " is invalid for a bioelectrics finite elasticity equation."
321  CALL flagerror(local_error,err,error,*999)
322  END SELECT
323  CASE(problem_setup_control_type)
324  SELECT CASE(problem_setup%ACTION_TYPE)
325  CASE(problem_setup_start_action)
326  !Set up a time control loop
327  CALL control_loop_create_start(problem,control_loop,err,error,*999)
328  CALL control_loop_type_set(control_loop,problem_control_time_loop_type,err,error,*999)
329  CALL control_loop_number_of_sub_loops_set(control_loop,2,err,error,*999)
330  CALL control_loop_output_type_set(control_loop,control_loop_progress_output,err,error,*999)
331 
332  !Set up the control sub loop for monodomain
333  CALL control_loop_sub_loop_get(control_loop,1,monodomain_sub_loop,err,error,*999)
334  CALL control_loop_label_set(monodomain_sub_loop,'MONODOMAIN_TIME_LOOP',err,error,*999)
335  CALL control_loop_type_set(monodomain_sub_loop,problem_control_time_loop_type,err,error,*999)
336  CALL control_loop_output_type_set(monodomain_sub_loop,control_loop_progress_output,err,error,*999)
337 
338  IF(problem%specification(3)==problem_monodomain_elasticity_velocity_subtype) THEN
339  !Set up the control sub loop for finite elasicity
340  CALL control_loop_sub_loop_get(control_loop,2,elasticity_sub_loop,err,error,*999)
341  CALL control_loop_label_set(elasticity_sub_loop,'ELASTICITY_WHILE_LOOP',err,error,*999)
342  CALL control_loop_type_set(elasticity_sub_loop,problem_control_while_loop_type,err,error,*999)
343  CALL control_loop_output_type_set(elasticity_sub_loop,control_loop_progress_output,err,error,*999)
344  ELSE
345  !Set up the control sub loop for finite elasicity
346  CALL control_loop_sub_loop_get(control_loop,2,elasticity_sub_loop,err,error,*999)
347  CALL control_loop_label_set(elasticity_sub_loop,'ELASTICITY_LOAD_INCREMENT_LOOP',err,error,*999)
348  CALL control_loop_type_set(elasticity_sub_loop,problem_control_load_increment_loop_type,err,error,*999)
349  CALL control_loop_output_type_set(elasticity_sub_loop,control_loop_progress_output,err,error,*999)
350  ENDIF
351  CASE(problem_setup_finish_action)
352  !Finish the control loops
353  control_loop_root=>problem%CONTROL_LOOP
354  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
355  CALL control_loop_create_finish(control_loop,err,error,*999)
356  !Sub-loops are finished when parent is finished
357  CASE DEFAULT
358  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
359  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
360  & " is invalid for a bioelectrics finite elasticity equation."
361  CALL flagerror(local_error,err,error,*999)
362  END SELECT
363  CASE(problem_setup_solvers_type)
364  !Get the control loop
365  control_loop_root=>problem%CONTROL_LOOP
366  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
367  SELECT CASE(problem_setup%ACTION_TYPE)
368  CASE(problem_setup_start_action)
369  !Get the monodomain sub loop
370  CALL control_loop_sub_loop_get(control_loop,1,monodomain_sub_loop,err,error,*999)
371  !Start the solvers creation
372  CALL solvers_create_start(monodomain_sub_loop,monodomain_solvers,err,error,*999)
373  CALL solvers_number_set(monodomain_solvers,2,err,error,*999)
374  !Set the first solver to be a differential-algebraic equations solver
375  NULLIFY(solver)
376  CALL solvers_solver_get(monodomain_solvers,1,solver,err,error,*999)
377  CALL solver_type_set(solver,solver_dae_type,err,error,*999)
378  CALL solver_label_set(solver,"ODE Solver",err,error,*999)
379  CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
380  !Set the second solver to be a dynamic solver
381  NULLIFY(solver)
382  CALL solvers_solver_get(monodomain_solvers,2,solver,err,error,*999)
383  CALL solver_type_set(solver,solver_dynamic_type,err,error,*999)
384  CALL solver_dynamic_order_set(solver,solver_dynamic_first_order,err,error,*999)
385  CALL solver_label_set(solver,"Parabolic solver",err,error,*999)
386  CALL solver_dynamic_degree_set(solver,solver_dynamic_first_degree,err,error,*999)
387  CALL solver_dynamic_scheme_set(solver,solver_dynamic_crank_nicolson_scheme,err,error,*999)
388  CALL solver_dynamic_restart_set(solver,.true.,err,error,*999)
389  CALL solver_library_type_set(solver,solver_cmiss_library,err,error,*999)
390 
391  !Get the finite elasticity sub loop
392  CALL control_loop_sub_loop_get(control_loop,2,elasticity_sub_loop,err,error,*999)
393  !Start the solvers creation
394  CALL solvers_create_start(elasticity_sub_loop,elasticity_solvers,err,error,*999)
395  CALL solvers_number_set(elasticity_solvers,1,err,error,*999)
396  !Set the finite elasticity solver to be a nonlinear solver
397  NULLIFY(solver)
398  CALL solvers_solver_get(elasticity_solvers,1,solver,err,error,*999)
399  CALL solver_type_set(solver,solver_nonlinear_type,err,error,*999)
400  CALL solver_library_type_set(solver,solver_petsc_library,err,error,*999)
401  CASE(problem_setup_finish_action)
402  !Get the monodomain solvers
403  CALL control_loop_sub_loop_get(control_loop,1,monodomain_sub_loop,err,error,*999)
404  CALL control_loop_solvers_get(monodomain_sub_loop,monodomain_solvers,err,error,*999)
405  !Finish the solvers creation
406  CALL solvers_create_finish(monodomain_solvers,err,error,*999)
407 
408  !Get the finite elasticity solvers
409  CALL control_loop_sub_loop_get(control_loop,2,elasticity_sub_loop,err,error,*999)
410  CALL control_loop_solvers_get(elasticity_sub_loop,elasticity_solvers,err,error,*999)
411  !Finish the solvers creation
412  CALL solvers_create_finish(elasticity_solvers,err,error,*999)
413  CASE DEFAULT
414  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
415  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
416  & " is invalid for a bioelectrics finite elasticity equation."
417  CALL flagerror(local_error,err,error,*999)
418  END SELECT
419  CASE(problem_setup_solver_equations_type)
420  SELECT CASE(problem_setup%ACTION_TYPE)
421  CASE(problem_setup_start_action)
422  !Get the control loop and solvers
423  control_loop_root=>problem%CONTROL_LOOP
424  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
425 
426  !Get the monodomain sub loop and solvers
427  CALL control_loop_sub_loop_get(control_loop,1,monodomain_sub_loop,err,error,*999)
428  CALL control_loop_solvers_get(monodomain_sub_loop,monodomain_solvers,err,error,*999)
429  !Create the solver equations for the second (parabolic) solver
430  NULLIFY(solver)
431  CALL solvers_solver_get(monodomain_solvers,2,solver,err,error,*999)
432  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
433  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
434  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_first_order_dynamic,err,error,*999)
435  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
436 
437  !Get the finite elasticity sub loop and solvers
438  CALL control_loop_sub_loop_get(control_loop,2,elasticity_sub_loop,err,error,*999)
439  CALL control_loop_solvers_get(elasticity_sub_loop,elasticity_solvers,err,error,*999)
440  !Get the finite elasticity solver and create the finite elasticity solver equations
441  NULLIFY(solver)
442  NULLIFY(solver_equations)
443  CALL solvers_solver_get(elasticity_solvers,1,solver,err,error,*999)
444  CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
445  CALL solver_equations_linearity_type_set(solver_equations,solver_equations_nonlinear,err,error,*999)
446  CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_static,err,error,*999)
447  CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
448  CASE(problem_setup_finish_action)
449  !Get the control loop
450  control_loop_root=>problem%CONTROL_LOOP
451  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
452 
453  !Get the monodomain sub loop and solvers
454  CALL control_loop_sub_loop_get(control_loop,1,monodomain_sub_loop,err,error,*999)
455  CALL control_loop_solvers_get(monodomain_sub_loop,monodomain_solvers,err,error,*999)
456  !Get the solver equations for the second (parabolic) solver
457  NULLIFY(solver)
458  NULLIFY(solver_equations)
459  CALL solvers_solver_get(monodomain_solvers,2,solver,err,error,*999)
460  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
461  !Finish the solver equations creation
462  CALL solver_equations_create_finish(solver_equations,err,error,*999)
463 
464  !Get the finite elasticity sub loop and solvers
465  CALL control_loop_sub_loop_get(control_loop,2,elasticity_sub_loop,err,error,*999)
466  CALL control_loop_solvers_get(elasticity_sub_loop,elasticity_solvers,err,error,*999)
467  !Finish the creation of the finite elasticity solver equations
468  NULLIFY(solver)
469  NULLIFY(solver_equations)
470  CALL solvers_solver_get(elasticity_solvers,1,solver,err,error,*999)
471  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
472  CALL solver_equations_create_finish(solver_equations,err,error,*999)
473  CASE DEFAULT
474  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
475  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
476  & " is invalid for a bioelectrics finite elasticity equation."
477  CALL flagerror(local_error,err,error,*999)
478  END SELECT
479  CASE(problem_setup_cellml_equations_type)
480  SELECT CASE(problem_setup%ACTION_TYPE)
481  CASE(problem_setup_start_action)
482  !Get the control loop
483  control_loop_root=>problem%CONTROL_LOOP
484  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
485  CALL control_loop_sub_loop_get(control_loop,1,monodomain_sub_loop,err,error,*999)
486  !Get the solvers
487  CALL control_loop_solvers_get(monodomain_sub_loop,monodomain_solvers,err,error,*999)
488  !Create the CellML equations for the first DAE solver
489  CALL solvers_solver_get(monodomain_solvers,1,solver,err,error,*999)
490  CALL cellml_equations_create_start(solver,cellml_equations,err,error,*999)
491  CASE(problem_setup_finish_action)
492  !Get the control loop
493  control_loop_root=>problem%CONTROL_LOOP
494  CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
495  CALL control_loop_sub_loop_get(control_loop,1,monodomain_sub_loop,err,error,*999)
496  !Get the solvers
497  CALL control_loop_solvers_get(monodomain_sub_loop,monodomain_solvers,err,error,*999)
498  !Get the CellML equations for the first DAE solver
499  CALL solvers_solver_get(monodomain_solvers,1,solver,err,error,*999)
500  CALL solver_cellml_equations_get(solver,cellml_equations,err,error,*999)
501  !Finish the CellML equations creation
502  CALL cellml_equations_create_finish(cellml_equations,err,error,*999)
503  CASE DEFAULT
504  local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
505  & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
506  & " is invalid for a bioelectrics finite elasticity equation."
507  CALL flagerror(local_error,err,error,*999)
508  END SELECT
509  CASE DEFAULT
510  local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
511  & " is invalid for a bioelectrics finite elasticity equation."
512  CALL flagerror(local_error,err,error,*999)
513  END SELECT
514  CASE DEFAULT
515  local_error="The problem subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
516  & " does not equal a transient monodomain quasistatic finite elasticity equation subtype."
517  CALL flagerror(local_error,err,error,*999)
518  END SELECT
519  ELSE
520  CALL flagerror("Problem is not associated.",err,error,*999)
521  ENDIF
522 
523  exits("BIOELECTRIC_FINITE_ELASTICITY_PROBLEM_SETUP")
524  RETURN
525 999 errorsexits("BIOELECTRIC_FINITE_ELASTICITY_PROBLEM_SETUP",err,error)
526  RETURN 1
527 
528  END SUBROUTINE bioelectric_finite_elasticity_problem_setup
529 
530  !
531  !================================================================================================================================
532  !
533 
535  SUBROUTINE bioelectric_finite_elasticity_pre_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
536 
537  !Argument variables
538  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
539  TYPE(solver_type), POINTER :: SOLVER
540  INTEGER(INTG), INTENT(OUT) :: ERR
541  TYPE(varying_string), INTENT(OUT) :: ERROR
542 
543  !Local Variables
544  TYPE(varying_string) :: LOCAL_ERROR
545 
546  enters("BIOELECTRIC_FINITE_ELASTICITY_PRE_SOLVE",err,error,*999)
547 
548  IF(ASSOCIATED(control_loop)) THEN
549  IF(ASSOCIATED(solver)) THEN
550  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
551  IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
552  CALL flagerror("Problem specification is not allocated.",err,error,*999)
553  ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
554  CALL flagerror("Problem specification must have three entries for a bioelectric-finite elasticity problem.", &
555  & err,error,*999)
556  END IF
557  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
558  CASE(problem_gudunov_monodomain_simple_elasticity_subtype,problem_gudunov_monodomain_1d3d_elasticity_subtype, &
559  & problem_monodomain_elasticity_w_titin_subtype,problem_monodomain_1d3d_active_strain_subtype)
560  SELECT CASE(control_loop%LOOP_TYPE)
561  CASE(problem_control_time_loop_type)
562  CALL biodomain_pre_solve(solver,err,error,*999)
563  CASE(problem_control_load_increment_loop_type)
564  CALL finite_elasticity_pre_solve(control_loop,solver,err,error,*999)
565  CASE DEFAULT
566  local_error="Control loop loop type "//trim(number_to_vstring(control_loop%LOOP_TYPE,"*",err,error))// &
567  & " is not valid for a bioelectrics finite elasticity type of a multi physics problem class."
568  CALL flagerror(local_error,err,error,*999)
569  END SELECT
570  CASE(problem_monodomain_elasticity_velocity_subtype)
571  SELECT CASE(control_loop%LOOP_TYPE)
572  CASE(problem_control_time_loop_type)
573  CALL biodomain_pre_solve(solver,err,error,*999)
574  CASE(problem_control_while_loop_type)
575  CALL finite_elasticity_pre_solve(control_loop,solver,err,error,*999)
576  CASE DEFAULT
577  local_error="Control loop loop type "//trim(number_to_vstring(control_loop%LOOP_TYPE,"*",err,error))// &
578  & " is not valid for a bioelectrics finite elasticity type of a multi physics problem class."
579  CALL flagerror(local_error,err,error,*999)
580  END SELECT
581  CASE DEFAULT
582  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
583  & " is not valid for a bioelectrics finite elasticity type of a multi physics problem class."
584  CALL flagerror(local_error,err,error,*999)
585  END SELECT
586  ELSE
587  CALL flagerror("Problem is not associated.",err,error,*999)
588  ENDIF
589  ELSE
590  CALL flagerror("Solver is not associated.",err,error,*999)
591  ENDIF
592  ELSE
593  CALL flagerror("Control loop is not associated.",err,error,*999)
594  ENDIF
595 
596  exits("BIOELECTRIC_FINITE_ELASTICITY_PRE_SOLVE")
597  RETURN
598 999 errorsexits("BIOELECTRIC_FINITE_ELASTICITY_PRE_SOLVE",err,error)
599  RETURN 1
600 
601  END SUBROUTINE bioelectric_finite_elasticity_pre_solve
602 
603  !
604  !================================================================================================================================
605  !
606 
608  SUBROUTINE bioelectric_finite_elasticity_post_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)
609 
610  !Argument variables
611  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
612  TYPE(solver_type), POINTER :: SOLVER
613  INTEGER(INTG), INTENT(OUT) :: ERR
614  TYPE(varying_string), INTENT(OUT) :: ERROR
615 
616  !Local Variables
617  TYPE(varying_string) :: LOCAL_ERROR
618 
619  enters("BIOELECTRIC_FINITE_ELASTICITY_POST_SOLVE",err,error,*999)
620 
621  IF(ASSOCIATED(control_loop)) THEN
622  IF(ASSOCIATED(solver)) THEN
623  IF(ASSOCIATED(control_loop%PROBLEM)) THEN
624  IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
625  CALL flagerror("Problem specification is not allocated.",err,error,*999)
626  ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
627  CALL flagerror("Problem specification must have three entries for a bioelectric-finite elasticity problem.", &
628  & err,error,*999)
629  END IF
630  SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
631  CASE(problem_gudunov_monodomain_simple_elasticity_subtype,problem_gudunov_monodomain_1d3d_elasticity_subtype, &
632  & problem_monodomain_elasticity_w_titin_subtype,problem_monodomain_elasticity_velocity_subtype, &
633  & problem_monodomain_1d3d_active_strain_subtype)
634  SELECT CASE(solver%SOLVE_TYPE)
635  CASE(solver_dae_type)
636  CALL bioelectric_post_solve(solver,err,error,*999)
637  CASE(solver_dynamic_type)
638  CALL bioelectric_post_solve(solver,err,error,*999)
639  CASE(solver_nonlinear_type)
640  CALL finite_elasticity_post_solve(control_loop,solver,err,error,*999)
641  CASE DEFAULT
642  local_error="Solver solve type "//trim(number_to_vstring(solver%SOLVE_TYPE,"*",err,error))// &
643  & " is not valid for a bioelectrics finite elasticity type of a multi physics problem class."
644  CALL flagerror(local_error,err,error,*999)
645  END SELECT
646  CASE DEFAULT
647  local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
648  & " is not valid for a bioelectrics finite elasticity type of a multi physics problem class."
649  CALL flagerror(local_error,err,error,*999)
650  END SELECT
651  ELSE
652  CALL flagerror("Problem is not associated.",err,error,*999)
653  ENDIF
654  ELSE
655  CALL flagerror("Solver is not associated.",err,error,*999)
656  ENDIF
657  ELSE
658  CALL flagerror("Control loop is not associated.",err,error,*999)
659  ENDIF
660 
661  exits("BIOELECTRIC_FINITE_ELASTICITY_POST_SOLVE")
662  RETURN
663 999 errorsexits("BIOELECTRIC_FINITE_ELASTICITY_POST_SOLVE",err,error)
664  RETURN 1
665 
666  END SUBROUTINE bioelectric_finite_elasticity_post_solve
667 
668  !
669  !================================================================================================================================
670  !
671 
673  SUBROUTINE bioelectricfiniteelasticity_controllooppreloop(CONTROL_LOOP,ERR,ERROR,*)
674 
675  !Argument variables
676  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
677  INTEGER(INTG), INTENT(OUT) :: ERR
678  TYPE(varying_string), INTENT(OUT) :: ERROR
679  !Local Variables
680  TYPE(problem_type), POINTER :: PROBLEM
681  TYPE(varying_string) :: LOCAL_ERROR
682 
683  enters("BioelectricFiniteElasticity_ControlLoopPreLoop",err,error,*999)
684 
685  IF(ASSOCIATED(control_loop)) THEN
686  IF(control_loop%NUMBER_OF_SUB_LOOPS==0) THEN
687  problem=>control_loop%PROBLEM
688  IF(ASSOCIATED(problem)) THEN
689  IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
690  CALL flagerror("Problem specification is not allocated.",err,error,*999)
691  ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
692  CALL flagerror("Problem specification must have three entries for a bioelectric-finite elasticity problem.", &
693  & err,error,*999)
694  END IF
695  SELECT CASE(problem%SPECIFICATION(2))
696  CASE(problem_bioelectric_finite_elasticity_type)
697  SELECT CASE(control_loop%LOOP_TYPE)
698  CASE(problem_control_time_loop_type)
699  !do nothing ???
700  CASE(problem_control_load_increment_loop_type)
701  SELECT CASE(problem%SPECIFICATION(3))
702  CASE(problem_gudunov_monodomain_1d3d_elasticity_subtype,problem_monodomain_1d3d_active_strain_subtype)
703  CALL bioelectricfiniteelasticity_independentfieldinterpolate(control_loop,err,error,*999)
704  CASE(problem_monodomain_elasticity_w_titin_subtype)
705  CALL bioelectric_finite_elasticity_compute_titin(control_loop,err,error,*999)
706  CALL bioelectricfiniteelasticity_independentfieldinterpolate(control_loop,err,error,*999)
707  CASE DEFAULT
708  local_error="The third problem specification of "// &
709  & trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
710  & " is not valid for bioelectric finite elasticity problem."
711  CALL flagerror(local_error,err,error,*999)
712  END SELECT
713  CALL finiteelasticity_controltimelooppreloop(control_loop,err,error,*999)
714  CASE(problem_control_while_loop_type)
715  SELECT CASE(problem%SPECIFICATION(3))
716  CASE(problem_monodomain_elasticity_velocity_subtype)
717  IF(control_loop%WHILE_LOOP%ITERATION_NUMBER==1) THEN
718  CALL bioelectricfiniteelasticity_independentfieldinterpolate(control_loop,err,error,*999)
719  ENDIF
720  CALL bioelectricfiniteelasticity_computefibrestretch(control_loop,err,error,*999)
721  CALL bioelectricfiniteelasticity_forcelengthvelocityrelation(control_loop,err,error,*999)
722  CASE DEFAULT
723  local_error="The third problem specification of "// &
724  & trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
725  & " is not valid for bioelectric finite elasticity problem."
726  CALL flagerror(local_error,err,error,*999)
727  END SELECT
728  CALL finiteelasticity_controltimelooppreloop(control_loop,err,error,*999)
729  CASE DEFAULT
730  local_error="Control loop loop type "//trim(number_to_vstring(control_loop%LOOP_TYPE,"*",err,error))// &
731  & " is not valid for bioelectric finite elasticity problem type."
732  CALL flagerror(local_error,err,error,*999)
733  END SELECT
734  CASE DEFAULT
735  local_error="The second problem specification of "// &
736  & trim(number_to_vstring(problem%SPECIFICATION(2),"*",err,error))// &
737  & " is not valid for a multi physics problem."
738  CALL flagerror(local_error,err,error,*999)
739  END SELECT
740  ELSE
741  CALL flagerror("Control loop problem is not associated.",err,error,*999)
742  ENDIF
743  ELSE
744  !the main time loop - do nothing!
745  ENDIF
746  ELSE
747  CALL flagerror("Control loop is not associated.",err,error,*999)
748  ENDIF
749 
750  exits("BioelectricFiniteElasticity_ControlLoopPreLoop")
751  RETURN
752 999 errors("BioelectricFiniteElasticity_ControlLoopPreLoop",err,error)
753  exits("BioelectricFiniteElasticity_ControlLoopPreLoop")
754  RETURN 1
755 
756  END SUBROUTINE bioelectricfiniteelasticity_controllooppreloop
757 
758  !
759  !================================================================================================================================
760  !
761 
763  SUBROUTINE bioelectricfiniteelasticity_computefibrestretch(CONTROL_LOOP,ERR,ERROR,*)
764 
765  !Argument variables
766  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
767  INTEGER(INTG), INTENT(OUT) :: ERR
768  TYPE(varying_string), INTENT(OUT) :: ERROR
769  !Local Variables
770  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
771  TYPE(field_type), POINTER :: DEPENDENT_FIELD,FIBRE_FIELD,GEOMETRIC_FIELD,INDEPENDENT_FIELD
772  TYPE(solver_type), POINTER :: SOLVER
773  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
774  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
775  TYPE(solvers_type), POINTER :: SOLVERS
776  TYPE(varying_string) :: LOCAL_ERROR
777  TYPE(basis_type), POINTER :: DEPENDENT_BASIS
778  TYPE(equations_type), POINTER :: EQUATIONS
779  TYPE(quadrature_scheme_type), POINTER :: DEPENDENT_QUADRATURE_SCHEME
780  TYPE(field_interpolation_parameters_type), POINTER :: GEOMETRIC_INTERPOLATION_PARAMETERS, &
781  & FIBRE_INTERPOLATION_PARAMETERS,DEPENDENT_INTERPOLATION_PARAMETERS
782  TYPE(field_interpolated_point_type), POINTER :: GEOMETRIC_INTERPOLATED_POINT,FIBRE_INTERPOLATED_POINT, &
783  & DEPENDENT_INTERPOLATED_POINT
784  TYPE(field_interpolated_point_metrics_type), POINTER :: GEOMETRIC_INTERPOLATED_POINT_METRICS, &
785  & DEPENDENT_INTERPOLATED_POINT_METRICS
786  TYPE(basis_type), POINTER :: GEOMETRIC_BASIS
787  TYPE(decomposition_type), POINTER :: DECOMPOSITION
788  TYPE(field_variable_type), POINTER :: FIELD_VAR_U1
789  INTEGER(INTG) :: DEPENDENT_NUMBER_OF_GAUSS_POINTS
790  INTEGER(INTG) :: MESH_COMPONENT_NUMBER,NUMBER_OF_ELEMENTS,FIELD_VAR_TYPE
791  INTEGER(INTG) :: equations_set_idx,gauss_idx,dof_idx,element_idx,idx
792  REAL(DP) :: DZDNU(3,3),DZDNUT(3,3),AZL(3,3)
793 
794  enters("BioelectricFiniteElasticity_ComputeFibreStretch",err,error,*999)
795 
796  NULLIFY(solvers)
797  NULLIFY(solver)
798  NULLIFY(solver_equations)
799  NULLIFY(field_var_u1)
800  NULLIFY(dependent_basis)
801  NULLIFY(equations)
802  NULLIFY(dependent_field,fibre_field,geometric_field,independent_field)
803  NULLIFY(dependent_quadrature_scheme)
804  NULLIFY(geometric_interpolation_parameters,fibre_interpolation_parameters,dependent_interpolation_parameters)
805  NULLIFY(geometric_interpolated_point,fibre_interpolated_point,dependent_interpolated_point)
806  NULLIFY(geometric_interpolated_point_metrics,dependent_interpolated_point_metrics)
807 
808  IF(ASSOCIATED(control_loop)) THEN
809  IF(control_loop%NUMBER_OF_SUB_LOOPS==0) THEN
810  SELECT CASE(control_loop%LOOP_TYPE)
811  CASE(problem_control_time_loop_type)
812  !do nothing
813  CASE(problem_control_load_increment_loop_type)
814  !do nothing
815  CASE(problem_control_while_loop_type)
816  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
817  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
818  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
819  !Loop over the equations sets associated with the solver
820  IF(ASSOCIATED(solver_equations)) THEN
821  solver_mapping=>solver_equations%SOLVER_MAPPING
822  IF(ASSOCIATED(solver_mapping)) THEN
823  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
824  equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%PTR
825  IF(ASSOCIATED(equations_set)) THEN
826  geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
827  IF(.NOT.ASSOCIATED(geometric_field)) THEN
828  local_error="Geometric field is not associated for equations set index "// &
829  & trim(number_to_vstring(equations_set_idx,"*",err,error))//" in the solver mapping."
830  CALL flagerror(local_error,err,error,*999)
831  ENDIF
832  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
833  IF(.NOT.ASSOCIATED(dependent_field)) THEN
834  local_error="Dependent field is not associated for equations set index "// &
835  & trim(number_to_vstring(equations_set_idx,"*",err,error))//" in the solver mapping."
836  CALL flagerror(local_error,err,error,*999)
837  ENDIF
838  independent_field=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
839  IF(.NOT.ASSOCIATED(independent_field)) THEN
840  local_error="Independent field is not associated for equations set index "// &
841  & trim(number_to_vstring(equations_set_idx,"*",err,error))//" in the solver mapping."
842  CALL flagerror(local_error,err,error,*999)
843  ENDIF
844  equations=>equations_set%EQUATIONS
845  IF(.NOT.ASSOCIATED(equations)) THEN
846  local_error="Equations is not associated for equations set index "// &
847  & trim(number_to_vstring(equations_set_idx,"*",err,error))//" in the solver mapping."
848  CALL flagerror(local_error,err,error,*999)
849  ENDIF
850  fibre_field=>equations%INTERPOLATION%FIBRE_FIELD
851  IF(.NOT.ASSOCIATED(fibre_field)) THEN
852  local_error="Fibre field is not associated for equations set index "// &
853  & trim(number_to_vstring(equations_set_idx,"*",err,error))//" in the solver mapping."
854  CALL flagerror(local_error,err,error,*999)
855  ENDIF
856  ELSE
857  local_error="Equations set is not associated for equations set index "// &
858  & trim(number_to_vstring(equations_set_idx,"*",err,error))//" in the solver mapping."
859  CALL flagerror(local_error,err,error,*999)
860  ENDIF
861  ENDDO !equations_set_idx
862  ELSE
863  CALL flagerror("Solver equations solver mapping is not associated.",err,error,*999)
864  ENDIF
865  ELSE
866  CALL flagerror("Solver solver equations are not associated.",err,error,*999)
867  ENDIF
868 
869  CALL field_variable_get(independent_field,field_u1_variable_type,field_var_u1,err,error,*999)
870 
871  decomposition=>dependent_field%DECOMPOSITION
872  mesh_component_number=decomposition%MESH_COMPONENT_NUMBER
873 
874  IF(control_loop%WHILE_LOOP%ITERATION_NUMBER==1) THEN
875  !copy the old fibre stretch to the previous values parameter set
876  CALL field_parameter_sets_copy(independent_field,field_u1_variable_type, &
877  & field_values_set_type,field_previous_values_set_type,1.0_dp,err,error,*999)
878  ENDIF
879 
880  number_of_elements=dependent_field%GEOMETRIC_FIELD%DECOMPOSITION%TOPOLOGY%ELEMENTS%NUMBER_OF_ELEMENTS
881 
882  !loop over the elements of the finite elasticity mesh (internal and boundary elements)
883  DO element_idx=1,number_of_elements
884 
885  dependent_basis=>decomposition%DOMAIN(mesh_component_number)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%BASIS
886  dependent_quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
887  dependent_number_of_gauss_points=dependent_quadrature_scheme%NUMBER_OF_GAUSS
888  geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
889  & topology%ELEMENTS%ELEMENTS(element_idx)%BASIS
890 
891  !Initialise tensors and matrices
892  dzdnu=0.0_dp
893  DO idx=1,3
894  dzdnu(idx,idx)=1.0_dp
895  ENDDO
896 
897  !Grab interpolation parameters
898  field_var_type=equations_set%EQUATIONS%EQUATIONS_MAPPING%NONLINEAR_MAPPING%RESIDUAL_VARIABLES(1)%PTR%VARIABLE_TYPE
899  dependent_interpolation_parameters=>equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR
900  geometric_interpolation_parameters=>equations%INTERPOLATION%GEOMETRIC_INTERP_PARAMETERS(field_u_variable_type)%PTR
901  fibre_interpolation_parameters=>equations%INTERPOLATION%FIBRE_INTERP_PARAMETERS(field_u_variable_type)%PTR
902 
903  CALL field_interpolation_parameters_element_get(field_values_set_type,element_idx, &
904  & geometric_interpolation_parameters,err,error,*999)
905  CALL field_interpolation_parameters_element_get(field_values_set_type,element_idx, &
906  & fibre_interpolation_parameters,err,error,*999)
907  CALL field_interpolation_parameters_element_get(field_values_set_type,element_idx, &
908  & dependent_interpolation_parameters,err,error,*999)
909 
910  !Point interpolation pointer
911  geometric_interpolated_point=>equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR
912  geometric_interpolated_point_metrics=>equations%INTERPOLATION% &
913  & geometric_interp_point_metrics(field_u_variable_type)%PTR
914  fibre_interpolated_point=>equations%INTERPOLATION%FIBRE_INTERP_POINT(field_u_variable_type)%PTR
915  dependent_interpolated_point=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_type)%PTR
916  dependent_interpolated_point_metrics=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT_METRICS(field_var_type)%PTR
917 
918  !Loop over gauss points
919  DO gauss_idx=1,dependent_number_of_gauss_points
920 
921  !Interpolate dependent, geometric, fibre and materials fields
922  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,gauss_idx, &
923  & dependent_interpolated_point,err,error,*999)
924  CALL field_interpolated_point_metrics_calculate(dependent_basis%NUMBER_OF_XI, &
925  & dependent_interpolated_point_metrics,err,error,*999)
926  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,gauss_idx, &
927  & geometric_interpolated_point,err,error,*999)
928  CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI, &
929  & geometric_interpolated_point_metrics,err,error,*999)
930  CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,gauss_idx, &
931  & fibre_interpolated_point,err,error,*999)
932 
933  !Calculate F=dZ/dNU, the deformation gradient tensor at the gauss point
934  CALL finiteelasticity_gaussdeformationgradienttensor(dependent_interpolated_point_metrics, &
935  & geometric_interpolated_point_metrics,fibre_interpolated_point,dzdnu,err,error,*999)
936 
937  !compute C=F^T F
938  CALL matrix_transpose(dzdnu,dzdnut,err,error,*999)
939  CALL matrix_product(dzdnut,dzdnu,azl,err,error,*999)
940 
941  !store the fibre stretch lambda_f, i.e., sqrt(C_11) or AZL(1,1)
942  dof_idx=field_var_u1%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,element_idx)
943  CALL field_parameter_set_update_local_dof(independent_field,field_u1_variable_type,field_values_set_type, &
944  & dof_idx,sqrt(azl(1,1)),err,error,*999)
945 
946  ENDDO !gauss_idx
947  ENDDO !element_idx
948 
949  !now the ghost elements -- get the relevant info from the other computational nodes
950  CALL field_parameter_set_update_start(independent_field,field_u1_variable_type,field_values_set_type,err,error,*999)
951  CALL field_parameter_set_update_finish(independent_field,field_u1_variable_type,field_values_set_type,err,error,*999)
952 
953  CASE DEFAULT
954  local_error="Control loop type "//trim(number_to_vstring(control_loop%LOOP_TYPE,"*",err,error))// &
955  & " is not valid for a bioelectrics finite elasticity type of a multi physics problem class."
956  CALL flagerror(local_error,err,error,*999)
957  END SELECT
958  ELSE !the control loop contains subloops
959  !do nothing
960  ENDIF
961  ELSE
962  CALL flagerror("Control loop is not associated.",err,error,*999)
963  ENDIF
964 
965  exits("BioelectricFiniteElasticity_ComputeFibreStretch")
966  RETURN
967 999 errors("BioelectricFiniteElasticity_ComputeFibreStretch",err,error)
968  exits("BioelectricFiniteElasticity_ComputeFibreStretch")
969  RETURN 1
970 
971  END SUBROUTINE bioelectricfiniteelasticity_computefibrestretch
972 
973  !
974  !================================================================================================================================
975  !
976 
978  SUBROUTINE bioelectricfiniteelasticity_forcelengthvelocityrelation(CONTROL_LOOP,ERR,ERROR,*)
979 
980  !Argument variables
981  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
982  INTEGER(INTG), INTENT(OUT) :: ERR
983  TYPE(varying_string), INTENT(OUT) :: ERROR
984  !Local Variables
985  TYPE(control_loop_type), POINTER :: CONTROL_LOOP_PARENT
986  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
987  TYPE(field_type), POINTER :: DEPENDENT_FIELD,FIBRE_FIELD,GEOMETRIC_FIELD,INDEPENDENT_FIELD
988  TYPE(solver_type), POINTER :: SOLVER
989  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
990  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
991  TYPE(solvers_type), POINTER :: SOLVERS
992  TYPE(varying_string) :: LOCAL_ERROR
993  TYPE(basis_type), POINTER :: DEPENDENT_BASIS
994  TYPE(equations_type), POINTER :: EQUATIONS
995  TYPE(quadrature_scheme_type), POINTER :: DEPENDENT_QUADRATURE_SCHEME
996  TYPE(field_interpolation_parameters_type), POINTER :: GEOMETRIC_INTERPOLATION_PARAMETERS, &
997  & FIBRE_INTERPOLATION_PARAMETERS,DEPENDENT_INTERPOLATION_PARAMETERS
998  TYPE(field_interpolated_point_type), POINTER :: GEOMETRIC_INTERPOLATED_POINT,FIBRE_INTERPOLATED_POINT, &
999  & DEPENDENT_INTERPOLATED_POINT
1000  TYPE(field_interpolated_point_metrics_type), POINTER :: GEOMETRIC_INTERPOLATED_POINT_METRICS, &
1001  & DEPENDENT_INTERPOLATED_POINT_METRICS
1002  TYPE(basis_type), POINTER :: GEOMETRIC_BASIS
1003  TYPE(decomposition_type), POINTER :: DECOMPOSITION
1004  TYPE(field_variable_type), POINTER :: FIELD_VAR_U,FIELD_VAR_U1
1005  INTEGER(INTG) :: DEPENDENT_NUMBER_OF_GAUSS_POINTS
1006  INTEGER(INTG) :: MESH_COMPONENT_NUMBER,NUMBER_OF_ELEMENTS
1007  INTEGER(INTG) :: equations_set_idx,gauss_idx,dof_idx,element_idx
1008  INTEGER(INTG) :: ITERATION_NUMBER,MAXIMUM_NUMBER_OF_ITERATIONS
1009  REAL(DP) :: LENGTH_HS,LENGTH_HS_0,ACTIVE_STRESS,FIBRE_STRETCH,FIBRE_STRETCH_OLD
1010  REAL(DP) :: FACTOR_LENGTH,FACTOR_VELO,SARCO_LENGTH,VELOCITY,VELOCITY_MAX,TIME_STEP,kappa,A,S,d,c
1011 
1012 !tomo
1013  REAL(DP) :: VELOCITY_AVERAGE,STRETCH_AVERAGE,OLD_STRETCH_AVERAGE
1014  INTEGER(INTG) :: counter
1015 
1016  enters("BioelectricFiniteElasticity_ForceLengthVelocityRelation",err,error,*999)
1017 
1018  NULLIFY(control_loop_parent)
1019  NULLIFY(equations_set)
1020  NULLIFY(solver_mapping)
1021  NULLIFY(geometric_basis)
1022  NULLIFY(decomposition)
1023  NULLIFY(solvers)
1024  NULLIFY(solver)
1025  NULLIFY(solver_equations)
1026  NULLIFY(field_var_u,field_var_u1)
1027  NULLIFY(dependent_basis)
1028  NULLIFY(equations)
1029  NULLIFY(dependent_field,fibre_field,geometric_field,independent_field)
1030  NULLIFY(dependent_quadrature_scheme)
1031  NULLIFY(geometric_interpolation_parameters,fibre_interpolation_parameters,dependent_interpolation_parameters)
1032  NULLIFY(geometric_interpolated_point,fibre_interpolated_point,dependent_interpolated_point)
1033  NULLIFY(geometric_interpolated_point_metrics,dependent_interpolated_point_metrics)
1034 
1035  IF(ASSOCIATED(control_loop)) THEN
1036  IF(control_loop%NUMBER_OF_SUB_LOOPS==0) THEN
1037  SELECT CASE(control_loop%LOOP_TYPE)
1038  CASE(problem_control_time_loop_type)
1039  !do nothing
1040  CASE(problem_control_load_increment_loop_type)
1041  !do nothing
1042  CASE(problem_control_while_loop_type)
1043  CALL control_loop_get(control_loop%PROBLEM%CONTROL_LOOP,control_loop_node,control_loop_parent,err,error,*999)
1044  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1045  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1046  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
1047  !Loop over the equations sets associated with the solver
1048  IF(ASSOCIATED(solver_equations)) THEN
1049  solver_mapping=>solver_equations%SOLVER_MAPPING
1050  IF(ASSOCIATED(solver_mapping)) THEN
1051  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
1052  equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%PTR
1053  IF(ASSOCIATED(equations_set)) THEN
1054  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
1055  IF(.NOT.ASSOCIATED(dependent_field)) THEN
1056  local_error="Dependent field is not associated for equations set index "// &
1057  & trim(number_to_vstring(equations_set_idx,"*",err,error))//" in the solver mapping."
1058  CALL flagerror(local_error,err,error,*999)
1059  ENDIF
1060  independent_field=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
1061  IF(.NOT.ASSOCIATED(independent_field)) THEN
1062  local_error="Independent field is not associated for equations set index "// &
1063  & trim(number_to_vstring(equations_set_idx,"*",err,error))//" in the solver mapping."
1064  CALL flagerror(local_error,err,error,*999)
1065  ENDIF
1066  ELSE
1067  local_error="Equations set is not associated for equations set index "// &
1068  & trim(number_to_vstring(equations_set_idx,"*",err,error))//" in the solver mapping."
1069  CALL flagerror(local_error,err,error,*999)
1070  ENDIF
1071  ENDDO !equations_set_idx
1072  ELSE
1073  CALL flagerror("Solver equations solver mapping is not associated.",err,error,*999)
1074  ENDIF
1075  ELSE
1076  CALL flagerror("Solver solver equations are not associated.",err,error,*999)
1077  ENDIF
1078 
1079  CALL field_variable_get(independent_field,field_u_variable_type,field_var_u,err,error,*999)
1080  CALL field_variable_get(independent_field,field_u1_variable_type,field_var_u1,err,error,*999)
1081 
1082  !get the inital half-sarcomere length
1083  dof_idx=field_var_u1%COMPONENTS(2)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
1084  CALL field_parameter_set_get_local_dof(independent_field,field_u1_variable_type, &
1085  & field_values_set_type,dof_idx,length_hs_0,err,error,*999)
1086 
1087  !get the maximum contraction velocity
1088  dof_idx=field_var_u1%COMPONENTS(3)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
1089  CALL field_parameter_set_get_local_dof(independent_field,field_u1_variable_type, &
1090  & field_values_set_type,dof_idx,velocity_max,err,error,*999)
1091 
1092  iteration_number=control_loop%WHILE_LOOP%ITERATION_NUMBER
1093  maximum_number_of_iterations=control_loop%WHILE_LOOP%MAXIMUM_NUMBER_OF_ITERATIONS
1094  !in the first iteration store the unaltered homogenized active stress field
1095  IF(iteration_number==1) THEN
1096  CALL field_parameter_sets_copy(independent_field,field_u_variable_type, &
1097  & field_values_set_type,field_previous_values_set_type,1.0_dp,err,error,*999)
1098  ELSE
1099  !restore the solution of the previous time step
1100  CALL field_parameter_sets_copy(dependent_field,field_u_variable_type, &
1101  & field_previous_values_set_type,field_values_set_type,1.0_dp,err,error,*999)
1102  ENDIF
1103 
1104  time_step=control_loop_parent%TIME_LOOP%TIME_INCREMENT
1105 
1106  decomposition=>dependent_field%DECOMPOSITION
1107  mesh_component_number=decomposition%MESH_COMPONENT_NUMBER
1108 
1109 !tomo start
1110  velocity_average=0.0_dp
1111  stretch_average=0.0_dp
1112  old_stretch_average=0.0_dp
1113  counter=0
1114 !tomo end
1115 
1116  number_of_elements=dependent_field%GEOMETRIC_FIELD%DECOMPOSITION%TOPOLOGY%ELEMENTS%NUMBER_OF_ELEMENTS
1117 
1118  !loop over the elements of the finite elasticity mesh (internal and boundary elements)
1119  DO element_idx=1,number_of_elements
1120 
1121  dependent_basis=>decomposition%DOMAIN(mesh_component_number)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%BASIS
1122  dependent_quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
1123  dependent_number_of_gauss_points=dependent_quadrature_scheme%NUMBER_OF_GAUSS
1124 
1125  !Loop over gauss points
1126  DO gauss_idx=1,dependent_number_of_gauss_points
1127 
1128  !get the unaltered ACTIVE_STRESS at the GP
1129  dof_idx=field_var_u%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,element_idx)
1130  CALL field_parameter_set_get_local_dof(independent_field,field_u_variable_type,field_previous_values_set_type, &
1131  & dof_idx,active_stress,err,error,*999)
1132 
1133  ! FORCE-LENGTH RELATION -------------------------------------------------------------------------------
1134 
1135  !get the current fibre stretch at the GP
1136  dof_idx=field_var_u1%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,element_idx)
1137  CALL field_parameter_set_get_local_dof(independent_field,field_u1_variable_type, &
1138  & field_values_set_type,dof_idx,fibre_stretch,err,error,*999)
1139 
1140  !compute the current half-sarcomere length at the GP: l_hs = lambda_f * l_hs_0
1141  length_hs=length_hs_0*fibre_stretch
1142 
1143  !compute the scale factor (0,1) due to sarcomere F-l relation of Gordon, A. M., A.F. Huxley, and F.J. Julian.
1144  !The variation in isometric tension with sarcomere length in vertebrate muscle fibres.
1145  !The Journal of Physiology 184.1 (1966): 170-192.
1146  sarco_length=2.0_dp*length_hs
1147  IF(sarco_length.LE.1.27_dp) THEN
1148  factor_length=0.0_dp
1149  ELSEIF(sarco_length.LE.1.7_dp) THEN
1150  factor_length=1.6047_dp*sarco_length-2.0379_dp
1151  ELSEIF(sarco_length.LE.2.34_dp) THEN
1152  factor_length=0.4844_dp*sarco_length-0.1334_dp
1153  ELSEIF(sarco_length.LE.2.51_dp) THEN
1154  factor_length=1.0_dp
1155  ELSEIF(sarco_length.LE.3.94_dp) THEN
1156  factor_length=-0.6993_dp*sarco_length+2.7552_dp
1157  ELSE
1158  factor_length=0.0_dp
1159  ENDIF
1160 
1161  !multiply the ACTIVE_STRESS with the scale factor
1162  active_stress=active_stress*factor_length
1163 
1164  !update the ACTIVE_STRESS at GP
1165  dof_idx=field_var_u%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,element_idx)
1166  CALL field_parameter_set_update_local_dof(independent_field,field_u_variable_type,field_values_set_type, &
1167  & dof_idx,active_stress,err,error,*999)
1168 
1169 
1170 
1171 
1172 
1173  ! FORCE-VELOCITY RELATION -------------------------------------------------------------------------------
1174 
1175  !get fibre stretch at the GP of the previous time step
1176  dof_idx=field_var_u1%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,element_idx)
1177  CALL field_parameter_set_get_local_dof(independent_field,field_u1_variable_type, &
1178  & field_previous_values_set_type,dof_idx,fibre_stretch_old,err,error,*999)
1179 
1180  !compute the contraction velocity
1181  velocity=(fibre_stretch-fibre_stretch_old)/time_step
1182 
1183 !!!! original
1184  !NOTE: VELOCITY_MAX is the max shortening velocity, and hence negative
1185  IF(velocity<velocity_max) THEN
1186  CALL flag_warning('Exceeded maximum contraction velocity (shortening).',err,error,*999)
1187 ! VELOCITY=VELOCITY_MAX
1188  !damping
1189  IF(iteration_number<(maximum_number_of_iterations/2)) THEN
1190  velocity=velocity*1.0_dp/dble((maximum_number_of_iterations/2)-iteration_number)
1191  ENDIF
1192  ELSEIF(velocity>(abs(velocity_max))) THEN
1193  CALL flag_warning('Exceeded maximum contraction velocity (lengthening).',err,error,*999)
1194 !!! VELOCITY=-VELOCITY_MAX
1195 ! !damping
1196 ! IF(ITERATION_NUMBER<(MAXIMUM_NUMBER_OF_ITERATIONS/2)) THEN
1197 ! VELOCITY=VELOCITY*1.0_DP/DBLE((MAXIMUM_NUMBER_OF_ITERATIONS/2)-ITERATION_NUMBER)
1198 ! ENDIF
1199  ENDIF
1200 
1201  !compute scale factor (FACTOR_VELO)
1202  kappa=0.24_dp !make mat param TODO
1203  a=1.6_dp !make mat param TODO
1204  s=2.5_dp !make mat param TODO
1205  IF(velocity<0.0_dp) THEN
1206  !shortening contraction
1207  factor_velo=(1-velocity/velocity_max)/(1+velocity/velocity_max/kappa)
1208  ELSE
1209  !lengthening contraction
1210  d=kappa*(1-a)
1211  c=velocity/velocity_max*s*(kappa+1)
1212  factor_velo=1+c*(a-1)/(d+c)
1213  ENDIF
1214 
1215  !multiply the ACTIVE_STRESS with the scale factor
1216  active_stress=active_stress*factor_velo
1217 
1218  !update the ACTIVE_STRESS at GP
1219  dof_idx=field_var_u%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,element_idx)
1220  CALL field_parameter_set_update_local_dof(independent_field,field_u_variable_type,field_values_set_type, &
1221  & dof_idx,active_stress,err,error,*999)
1222 
1223  ENDDO !gauss_idx
1224  ENDDO !element_idx
1225 !!!!end original
1226 
1227 
1228 
1229 
1230 !!!!tomo --- try averaging
1231 !!! VELOCITY_AVERAGE=VELOCITY_AVERAGE+VELOCITY
1232 !!! STRETCH_AVERAGE=STRETCH_AVERAGE+FIBRE_STRETCH
1233 !!! OLD_STRETCH_AVERAGE=OLD_STRETCH_AVERAGE+FIBRE_STRETCH_OLD
1234 !!!
1235 !!! counter=counter+1
1236 
1237 !!! ENDDO !gauss_idx
1238 !!! ENDDO !element_idx
1239 
1240 !!! VELOCITY=VELOCITY_AVERAGE/REAL(counter)
1241 !!! FIBRE_STRETCH=STRETCH_AVERAGE/REAL(counter)
1242 !!! FIBRE_STRETCH_OLD=OLD_STRETCH_AVERAGE/REAL(counter)
1243 
1244 !!! VELOCITY_AVERAGE=(FIBRE_STRETCH-FIBRE_STRETCH_OLD)/TIME_STEP
1245 
1246 !!!! VELOCITY=0.0_DP
1247 
1248 !!!! !damping
1249 !!!! IF(ITERATION_NUMBER<(MAXIMUM_NUMBER_OF_ITERATIONS/2)) THEN
1250 !!!! VELOCITY=VELOCITY*1.0_DP/DBLE((MAXIMUM_NUMBER_OF_ITERATIONS/2)-ITERATION_NUMBER)
1251 !!!! ENDIF
1252 
1253 !!! LOCAL_ERROR="######### VELOCITY: "//TRIM(NUMBER_TO_VSTRING(VELOCITY,"*",ERR,ERROR))//" #########"
1254 !!! CALL FLAG_WARNING(LOCAL_ERROR,ERR,ERROR,*999)
1255 !!! LOCAL_ERROR="######### VELOCITY AVERAGE: "//TRIM(NUMBER_TO_VSTRING(VELOCITY_AVERAGE,"*",ERR,ERROR))//" #########"
1256 !!! CALL FLAG_WARNING(LOCAL_ERROR,ERR,ERROR,*999)
1257 !!! LOCAL_ERROR="######### STRETCH: "//TRIM(NUMBER_TO_VSTRING(FIBRE_STRETCH,"*",ERR,ERROR))//" #########"
1258 !!! CALL FLAG_WARNING(LOCAL_ERROR,ERR,ERROR,*999)
1259 !!! LOCAL_ERROR="######### OLD STRETCH: "//TRIM(NUMBER_TO_VSTRING(FIBRE_STRETCH_OLD,"*",ERR,ERROR))//" #########"
1260 !!! CALL FLAG_WARNING(LOCAL_ERROR,ERR,ERROR,*999)
1261 
1262 
1263 !!!
1264 !!!
1265 !!! VELOCITY=VELOCITY_AVERAGE
1266 !!!
1267 !!!
1268 !!!
1269 
1270 !!! !compute scale factor (FACTOR_VELO)
1271 !!! kappa=0.24_DP !make mat param TODO
1272 !!! A=1.6_DP !make mat param TODO
1273 !!! S=2.5_DP !make mat param TODO
1274 !!! IF(VELOCITY<0.0_DP) THEN
1275 !!! !shortening contraction
1276 !!! FACTOR_VELO=(1-VELOCITY/VELOCITY_MAX)/(1+VELOCITY/VELOCITY_MAX/kappa)
1277 !!! IF(VELOCITY<VELOCITY_MAX) THEN
1278 !!! CALL FLAG_WARNING('Exceeded maximum contraction velocity (shortening).',ERR,ERROR,*999)
1279 !!! ENDIF
1280 !!! ELSE
1281 !!! !lengthening contraction
1282 !!! d=kappa*(1-A)
1283 !!! c=VELOCITY/VELOCITY_MAX*S*(kappa+1)
1284 !!! FACTOR_VELO=1+c*(A-1)/(d+c)
1285 !!! IF(VELOCITY>ABS(VELOCITY_MAX)) THEN
1286 !!! CALL FLAG_WARNING('Exceeded maximum contraction velocity (lengthening).',ERR,ERROR,*999)
1287 !!! ENDIF
1288 !!! ENDIF
1289 
1290 !!! DO element_idx=1,NUMBER_OF_ELEMENTS
1291 
1292 !!! DEPENDENT_BASIS=>DECOMPOSITION%DOMAIN(MESH_COMPONENT_NUMBER)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(element_idx)%BASIS
1293 !!! DEPENDENT_QUADRATURE_SCHEME=>DEPENDENT_BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP(BASIS_DEFAULT_QUADRATURE_SCHEME)%PTR
1294 !!! DEPENDENT_NUMBER_OF_GAUSS_POINTS=DEPENDENT_QUADRATURE_SCHEME%NUMBER_OF_GAUSS
1295 
1296 !!! !Loop over gauss points
1297 !!! DO gauss_idx=1,DEPENDENT_NUMBER_OF_GAUSS_POINTS
1298 
1299 !!! !get the ACTIVE_STRESS at the GP
1300 !!! dof_idx=FIELD_VAR_U%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,element_idx)
1301 !!! CALL FIELD_PARAMETER_SET_GET_LOCAL_DOF(INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE, &
1302 !!! & dof_idx,ACTIVE_STRESS,ERR,ERROR,*999)
1303 
1304 !!! !multiply the ACTIVE_STRESS with the scale factor
1305 !!! ACTIVE_STRESS=ACTIVE_STRESS*FACTOR_VELO
1306 
1307 !!! !update the ACTIVE_STRESS at GP
1308 !!! dof_idx=FIELD_VAR_U%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,element_idx)
1309 !!! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_DOF(INDEPENDENT_FIELD,FIELD_U_VARIABLE_TYPE,FIELD_VALUES_SET_TYPE, &
1310 !!! & dof_idx,ACTIVE_STRESS,ERR,ERROR,*999)
1311 !!!
1312 !!! ENDDO !gauss_idx
1313 !!! ENDDO !element_idx
1314 !tomo end
1315 
1316  !now the ghost elements -- get the relevant info from the other computational nodes
1317  CALL field_parameter_set_update_start(independent_field,field_u_variable_type,field_values_set_type,err,error,*999)
1318  CALL field_parameter_set_update_finish(independent_field,field_u_variable_type,field_values_set_type,err,error,*999)
1319 
1320  CASE DEFAULT
1321  local_error="Control loop type "//trim(number_to_vstring(control_loop%LOOP_TYPE,"*",err,error))// &
1322  & " is not valid for a bioelectrics finite elasticity type of a multi physics problem class."
1323  CALL flagerror(local_error,err,error,*999)
1324  END SELECT
1325  ELSE !the control loop contains subloops
1326  !do nothing
1327  ENDIF
1328  ELSE
1329  CALL flagerror("Control loop is not associated.",err,error,*999)
1330  ENDIF
1331 
1332 
1333  exits("BioelectricFiniteElasticity_ForceLengthVelocityRelation")
1334  RETURN
1335 999 errors("BioelectricFiniteElasticity_ForceLengthVelocityRelation",err,error)
1336  exits("BioelectricFiniteElasticity_ForceLengthVelocityRelation")
1337  RETURN 1
1338 
1339  END SUBROUTINE bioelectricfiniteelasticity_forcelengthvelocityrelation
1340 
1341  !
1342  !================================================================================================================================
1343  !
1344 
1346  SUBROUTINE bioelectricfiniteelasticity_controllooppostloop(CONTROL_LOOP,ERR,ERROR,*)
1347 
1348  !Argument variables
1349  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
1350  INTEGER(INTG), INTENT(OUT) :: ERR
1351  TYPE(varying_string), INTENT(OUT) :: ERROR
1352  !Local Variables
1353  TYPE(problem_type), POINTER :: PROBLEM
1354  INTEGER(INTG) :: equations_set_idx
1355  TYPE(control_loop_time_type), POINTER :: TIME_LOOP
1356  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1357  TYPE(field_type), POINTER :: DEPENDENT_FIELD
1358  TYPE(region_type), POINTER :: DEPENDENT_REGION
1359  TYPE(solver_type), POINTER :: SOLVER
1360  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
1361  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
1362  TYPE(solvers_type), POINTER :: SOLVERS
1363  TYPE(varying_string) :: FILENAME,LOCAL_ERROR,METHOD
1364  TYPE(control_loop_type), POINTER :: ELASTICITY_SUB_LOOP,BIOELECTRIC_SUB_LOOP
1365 
1366  enters("BioelectricFiniteElasticity_ControlLoopPostLoop",err,error,*999)
1367 
1368  IF(ASSOCIATED(control_loop)) THEN
1369  IF(control_loop%NUMBER_OF_SUB_LOOPS==0) THEN
1370  problem=>control_loop%PROBLEM
1371  IF(ASSOCIATED(problem)) THEN
1372  IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
1373  CALL flagerror("Problem specification is not allocated.",err,error,*999)
1374  ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
1375  CALL flagerror("Problem specification must have three entries for a bioelectric-finite elasticity problem.", &
1376  & err,error,*999)
1377  END IF
1378  SELECT CASE(control_loop%LOOP_TYPE)
1379  CASE(problem_control_time_loop_type)
1380  SELECT CASE(problem%SPECIFICATION(2))
1381  CASE(problem_bioelectric_finite_elasticity_type)
1382  !the monodomain time loop - output of the monodomain fields
1383  CALL biodomain_control_loop_post_loop(control_loop,err,error,*999)
1384  CASE DEFAULT
1385  local_error="Problem type "//trim(number_to_vstring(problem%SPECIFICATION(2),"*",err,error))// &
1386  & " is not valid for a multi physics problem class."
1387  CALL flagerror(local_error,err,error,*999)
1388  END SELECT
1389  CASE(problem_control_load_increment_loop_type)
1390  CALL bioelectricfiniteelasticity_updategeometricfield(control_loop,.false.,err,error,*999)
1391  CASE(problem_control_while_loop_type)
1392  CALL bioelectricfiniteelasticity_convergencecheck(control_loop,err,error,*999)
1393  CASE DEFAULT
1394  !do nothing
1395  END SELECT
1396  ELSE
1397  CALL flagerror("Control loop problem is not associated.",err,error,*999)
1398  ENDIF
1399  ELSE
1400  !the main time loop - output the finite elasticity fields
1401  IF(control_loop%OUTPUT_TYPE>=control_loop_progress_output) THEN
1402  !Export the dependent field for this time step
1403  time_loop=>control_loop%TIME_LOOP
1404  IF(ASSOCIATED(time_loop)) THEN
1405  problem=>control_loop%PROBLEM
1406  IF(ASSOCIATED(problem)) THEN
1407  NULLIFY(solvers)
1408  NULLIFY(solver)
1409  NULLIFY(solver_equations)
1410  NULLIFY(elasticity_sub_loop)
1411  !Get the solver. The first solver of the second sub loop will contain the finite elasticity dependent field equation set
1412  CALL control_loop_sub_loop_get(control_loop,2,elasticity_sub_loop,err,error,*999)
1413  CALL control_loop_solvers_get(elasticity_sub_loop,solvers,err,error,*999)
1414  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1415  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
1416  !Loop over the equations sets associated with the solver
1417  IF(ASSOCIATED(solver_equations)) THEN
1418  solver_mapping=>solver_equations%SOLVER_MAPPING
1419  IF(ASSOCIATED(solver_mapping)) THEN
1420  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
1421  equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%PTR
1422  IF(ASSOCIATED(equations_set)) THEN
1423  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
1424  NULLIFY(dependent_region)
1425  CALL field_region_get(dependent_field,dependent_region,err,error,*999)
1426  filename="MainTime_"//trim(number_to_vstring(dependent_region%USER_NUMBER,"*",err,error))// &
1427  & "_"//trim(number_to_vstring(time_loop%GLOBAL_ITERATION_NUMBER,"*",err,error))
1428  method="FORTRAN"
1429  CALL field_io_nodes_export(dependent_region%FIELDS,filename,method,err,error,*999)
1430  ELSE
1431  local_error="Equations set is not associated for equations set index "// &
1432  & trim(number_to_vstring(equations_set_idx,"*",err,error))// &
1433  & " in the solver mapping."
1434  CALL flagerror(local_error,err,error,*999)
1435  ENDIF
1436  ENDDO !equations_set_idx
1437  ELSE
1438  CALL flagerror("Solver equations solver mapping is not associated.",err,error,*999)
1439  ENDIF
1440  ELSE
1441  CALL flagerror("Solver solver equations are not associated.",err,error,*999)
1442  ENDIF
1443  IF((problem%SPECIFICATION(3)==problem_gudunov_monodomain_1d3d_elasticity_subtype).OR. &
1444  & (problem%SPECIFICATION(3)==problem_monodomain_elasticity_w_titin_subtype).OR. &
1445  & (problem%SPECIFICATION(3)==problem_monodomain_elasticity_velocity_subtype).OR. &
1446  & (problem%SPECIFICATION(3)==problem_monodomain_1d3d_active_strain_subtype)) THEN
1447  NULLIFY(solvers)
1448  NULLIFY(solver)
1449  NULLIFY(solver_equations)
1450  NULLIFY(bioelectric_sub_loop)
1451  !Get the solver. The second solver of the first sub loop will contain the bioelectrics equation set
1452  CALL control_loop_sub_loop_get(control_loop,1,bioelectric_sub_loop,err,error,*999)
1453  CALL control_loop_solvers_get(bioelectric_sub_loop,solvers,err,error,*999)
1454  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
1455  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
1456  !Loop over the equations sets associated with the solver
1457  IF(ASSOCIATED(solver_equations)) THEN
1458  solver_mapping=>solver_equations%SOLVER_MAPPING
1459  IF(ASSOCIATED(solver_mapping)) THEN
1460  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
1461  equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%PTR
1462  IF(ASSOCIATED(equations_set)) THEN
1463  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
1464  NULLIFY(dependent_region)
1465  CALL field_region_get(dependent_field,dependent_region,err,error,*999)
1466  filename="MainTime_M_"//trim(number_to_vstring(dependent_region%USER_NUMBER,"*",err,error))// &
1467  & "_"//trim(number_to_vstring(time_loop%GLOBAL_ITERATION_NUMBER,"*",err,error))
1468  method="FORTRAN"
1469  CALL field_io_nodes_export(dependent_region%FIELDS,filename,method,err,error,*999)
1470 
1471  WRITE(*,*) time_loop%ITERATION_NUMBER
1472 
1473  ELSE
1474  local_error="Equations set is not associated for equations set index "// &
1475  & trim(number_to_vstring(equations_set_idx,"*",err,error))// &
1476  & " in the solver mapping."
1477  CALL flagerror(local_error,err,error,*999)
1478  ENDIF
1479  ENDDO !equations_set_idx
1480  ELSE
1481  CALL flagerror("Solver equations solver mapping is not associated.",err,error,*999)
1482  ENDIF
1483  ELSE
1484  CALL flagerror("Solver solver equations are not associated.",err,error,*999)
1485  ENDIF
1486  ENDIF !PROBLEM_GUDUNOV_MONODOMAIN_1D3D_ELASTICITY_SUBTYPE,PROBLEM_MONODOMAIN_ELASTICITY_W_TITIN_SUBTYPE
1487  ELSE
1488  CALL flagerror("Control loop problem is not associated.",err,error,*999)
1489  ENDIF
1490  ELSE
1491  CALL flagerror("Time loop is not associated.",err,error,*999)
1492  ENDIF
1493  ENDIF
1494  ENDIF
1495  ELSE
1496  CALL flagerror("Control loop is not associated.",err,error,*999)
1497  ENDIF
1498 
1499  exits("BioelectricFiniteElasticity_ControlLoopPostLoop")
1500  RETURN
1501 999 errors("BioelectricFiniteElasticity_ControlLoopPostLoop",err,error)
1502  exits("BioelectricFiniteElasticity_ControlLoopPostLoop")
1503  RETURN 1
1504 
1505  END SUBROUTINE bioelectricfiniteelasticity_controllooppostloop
1506 
1507 
1508  !
1509  !================================================================================================================================
1510  !
1511 
1513  SUBROUTINE bioelectricfiniteelasticity_convergencecheck(CONTROL_LOOP,ERR,ERROR,*)
1514 
1515  !Argument variables
1516  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
1517  INTEGER(INTG), INTENT(OUT) :: ERR
1518  TYPE(varying_string), INTENT(OUT) :: ERROR
1519  !Local Variables
1520  TYPE(problem_type), POINTER :: PROBLEM
1521  INTEGER(INTG) :: equations_set_idx,NUMBER_OF_NODES,dof_idx,node_idx
1522  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1523  TYPE(field_type), POINTER :: DEPENDENT_FIELD
1524  TYPE(solvers_type), POINTER :: SOLVERS
1525  TYPE(solver_type), POINTER :: SOLVER
1526  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
1527  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
1528  TYPE(varying_string) :: LOCAL_ERROR
1529  TYPE(field_variable_type), POINTER :: FIELD_VAR
1530  REAL(DP) :: x1,x2,x3,y1,y2,y3,my_sum
1531 
1532  enters("BioelectricFiniteElasticity_ConvergenceCheck",err,error,*999)
1533 
1534  NULLIFY(solvers)
1535  NULLIFY(solver)
1536  NULLIFY(solver_equations)
1537  NULLIFY(field_var)
1538 
1539  IF(ASSOCIATED(control_loop)) THEN
1540  IF(control_loop%NUMBER_OF_SUB_LOOPS==0) THEN
1541  problem=>control_loop%PROBLEM
1542  IF(ASSOCIATED(problem)) THEN
1543  SELECT CASE(control_loop%LOOP_TYPE)
1544  CASE(problem_control_time_loop_type)
1545  !do nothing
1546  CASE(problem_control_load_increment_loop_type)
1547  !do nothing
1548  CASE(problem_control_while_loop_type)
1549  CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
1550  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1551  CALL solver_solution_update(solver,err,error,*999) !tomo added this
1552  CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
1553  !Loop over the equations sets associated with the solver
1554  IF(ASSOCIATED(solver_equations)) THEN
1555  solver_mapping=>solver_equations%SOLVER_MAPPING
1556  IF(ASSOCIATED(solver_mapping)) THEN
1557  DO equations_set_idx=1,solver_mapping%NUMBER_OF_EQUATIONS_SETS
1558  equations_set=>solver_mapping%EQUATIONS_SETS(equations_set_idx)%PTR
1559  IF(ASSOCIATED(equations_set)) THEN
1560  dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
1561  ELSE
1562  local_error="Equations set is not associated for equations set index "// &
1563  & trim(number_to_vstring(equations_set_idx,"*",err,error))// &
1564  & " in the solver mapping."
1565  CALL flagerror(local_error,err,error,*999)
1566  ENDIF
1567  ENDDO !equations_set_idx
1568  ELSE
1569  CALL flagerror("Solver equations solver mapping is not associated.",err,error,*999)
1570  ENDIF
1571  ELSE
1572  CALL flagerror("Solver solver equations are not associated.",err,error,*999)
1573  ENDIF
1574 
1575  IF(control_loop%WHILE_LOOP%ITERATION_NUMBER==1) THEN
1576  !
1577  ELSE
1578  CALL field_variable_get(dependent_field,field_u_variable_type,field_var,err,error,*999)
1579 
1580  number_of_nodes=dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION% &
1581  & mesh_component_number)%PTR%MAPPINGS%NODES%NUMBER_OF_LOCAL
1582 
1583  my_sum=0.0_dp
1584  DO node_idx=1,number_of_nodes
1585 
1586  !get the current node position (x) and the node position of the last iteration (y)
1587  dof_idx=field_var%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)% &
1588  & derivatives(1)%VERSIONS(1)
1589  CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type, &
1590  & field_values_set_type,dof_idx,x1,err,error,*999)
1591  CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type, &
1592  & field_previous_iteration_values_set_type,dof_idx,y1,err,error,*999)
1593  dof_idx=field_var%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)% &
1594  & derivatives(1)%VERSIONS(1)
1595  CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type, &
1596  & field_values_set_type,dof_idx,x2,err,error,*999)
1597  CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type, &
1598  & field_previous_iteration_values_set_type,dof_idx,y2,err,error,*999)
1599  dof_idx=field_var%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)% &
1600  & derivatives(1)%VERSIONS(1)
1601  CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type, &
1602  & field_values_set_type,dof_idx,x3,err,error,*999)
1603  CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type, &
1604  & field_previous_iteration_values_set_type,dof_idx,y3,err,error,*999)
1605 
1606  !sum up
1607  my_sum=my_sum+sqrt((x1-y1)**2+(x2-y2)**2+(x3-y3)**2)
1608  ENDDO
1609 
1610  IF(my_sum<1.0e-06_dp) THEN !if converged then:
1611  control_loop%WHILE_LOOP%CONTINUE_LOOP=.false.
1612  CALL bioelectricfiniteelasticity_updategeometricfield(control_loop,.false.,err,error,*999)
1613  !copy the current solution to the previous solution
1614  CALL field_parameter_sets_copy(dependent_field,field_u_variable_type, &
1615  & field_values_set_type,field_previous_values_set_type,1.0_dp,err,error,*999)
1616  ELSEIF(control_loop%WHILE_LOOP%ITERATION_NUMBER==control_loop%WHILE_LOOP%MAXIMUM_NUMBER_OF_ITERATIONS) THEN
1617  CALL bioelectricfiniteelasticity_updategeometricfield(control_loop,.false.,err,error,*999)
1618  CALL flag_warning('----------- Maximum number of iterations in while loop reached. -----------',err,error,*999)
1619  !copy the current solution to the previous solution
1620  CALL field_parameter_sets_copy(dependent_field,field_u_variable_type, &
1621  & field_values_set_type,field_previous_values_set_type,1.0_dp,err,error,*999)
1622  ENDIF
1623 
1624  ENDIF
1625 
1626  !copy the current solution to the previous iteration solution
1627  CALL field_parameter_sets_copy(dependent_field,field_u_variable_type, &
1628  & field_values_set_type,field_previous_iteration_values_set_type,1.0_dp,err,error,*999)
1629 
1630  CASE DEFAULT
1631  !do nothing
1632  END SELECT
1633  ELSE
1634  CALL flagerror("Control loop problem is not associated.",err,error,*999)
1635  ENDIF
1636  ELSE !the main time loop
1637  !do nothing
1638  ENDIF
1639  ELSE
1640  CALL flagerror("Control loop is not associated.",err,error,*999)
1641  ENDIF
1642 
1643  exits("BioelectricFiniteElasticity_ConvergenceCheck")
1644  RETURN
1645 999 errors("BioelectricFiniteElasticity_ConvergenceCheck",err,error)
1646  exits("BioelectricFiniteElasticity_ConvergenceCheck")
1647  RETURN 1
1648 
1649  END SUBROUTINE bioelectricfiniteelasticity_convergencecheck
1650 
1651  !
1652  !================================================================================================================================
1653  !
1654 
1657  SUBROUTINE bioelectricfiniteelasticity_updategeometricfield(CONTROL_LOOP,CALC_CLOSEST_GAUSS_POINT,ERR,ERROR,*)
1658 
1659  !Argument variables
1660  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
1661  LOGICAL, INTENT(IN) :: CALC_CLOSEST_GAUSS_POINT
1662  INTEGER(INTG), INTENT(OUT) :: ERR
1663  TYPE(varying_string), INTENT(OUT) :: ERROR
1664  !Local Variables
1665  TYPE(control_loop_type), POINTER :: CONTROL_LOOP_ROOT,CONTROL_LOOP_PARENT,CONTROL_LOOP_ELASTICITY,CONTROL_LOOP_MONODOMAIN
1666  TYPE(problem_type), POINTER :: PROBLEM
1667  TYPE(solvers_type), POINTER :: SOLVERS
1668  TYPE(solver_type), POINTER :: SOLVER
1669  TYPE(field_type), POINTER :: INDEPENDENT_FIELD_ELASTICITY,GEOMETRIC_FIELD_MONODOMAIN,GEOMETRIC_FIELD_ELASTICITY
1670  TYPE(field_type), POINTER :: DEPENDENT_FIELD_MONODOMAIN,INDEPENDENT_FIELD_MONODOMAIN,DEPENDENT_FIELD_ELASTICITY
1671  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
1672  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
1673  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
1674  TYPE(varying_string) :: LOCAL_ERROR
1675  TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
1676  TYPE(field_interpolation_parameters_type), POINTER :: INTERPOLATION_PARAMETERS
1677  TYPE(domain_mapping_type), POINTER :: NODES_MAPPING
1678  TYPE(decomposition_elements_type), POINTER :: ELEMENTS_TOPOLOGY
1679  TYPE(field_variable_type), POINTER :: FIELD_VAR_DEP_M,FIELD_VAR_GEO_M,FIELD_VAR_IND_FE,FIELD_VAR_IND_M,FIELD_VAR_IND_M_2
1680  INTEGER(INTG) :: component_idx,element_idx,ne,start_elem,START_ELEMENT,start_element_idx
1681  INTEGER(INTG) :: DEPENDENT_FIELD_INTERPOLATION,GEOMETRIC_FIELD_INTERPOLATION
1682  INTEGER(INTG) :: node_idx,node_idx_2,GAUSS_POINT,gauss_idx,fibre_idx
1683  INTEGER(INTG) :: nodes_in_Xi_1,nodes_in_Xi_2,nodes_in_Xi_3,n3,n2,n1,dof_idx,dof_idx2,idx,my_element_idx
1684  INTEGER(INTG) :: offset,n4
1685  REAL(DP) :: XVALUE_M,XVALUE_FE,DIST_LEFT,DIST_RIGHT,VALUE,VALUE_LEFT,VALUE_RIGHT,DISTANCE,VELOCITY,VELOCITY_MAX,OLD_DIST
1686  REAL(DP) :: OLD_DIST_2,OLD_DIST_3,OLD_DIST_4
1687  REAL(DP) :: XI(3),PREVIOUS_NODE(3),DIST_INIT,SARCO_LENGTH_INIT,TIME_STEP,DIST
1688  LOGICAL :: OUTSIDE_NODE
1689  REAL(DP), POINTER :: GAUSS_POSITIONS(:,:)
1690 
1691  enters("BioelectricFiniteElasticity_UpdateGeometricField",err,error,*999)
1692 
1693  NULLIFY(control_loop_root)
1694  NULLIFY(control_loop_parent)
1695  NULLIFY(control_loop_elasticity)
1696  NULLIFY(control_loop_monodomain)
1697  NULLIFY(problem)
1698  NULLIFY(solvers)
1699  NULLIFY(solver)
1700  NULLIFY(independent_field_elasticity)
1701  NULLIFY(dependent_field_monodomain)
1702  NULLIFY(independent_field_monodomain)
1703  NULLIFY(dependent_field_elasticity)
1704  NULLIFY(solver_equations)
1705  NULLIFY(solver_mapping)
1706  NULLIFY(equations_set)
1707  NULLIFY(geometric_field_monodomain)
1708  NULLIFY(geometric_field_elasticity)
1709  NULLIFY(elements_topology)
1710  NULLIFY(interpolated_point)
1711  NULLIFY(interpolation_parameters)
1712  NULLIFY(field_var_dep_m)
1713  NULLIFY(field_var_geo_m)
1714  NULLIFY(field_var_ind_fe)
1715  NULLIFY(field_var_ind_m)
1716  NULLIFY(field_var_ind_m_2)
1717  NULLIFY(gauss_positions)
1718 
1719  IF(ASSOCIATED(control_loop)) THEN
1720  IF(control_loop%NUMBER_OF_SUB_LOOPS==0) THEN
1721  problem=>control_loop%PROBLEM
1722  IF(ASSOCIATED(problem)) THEN
1723  IF(.NOT.ALLOCATED(problem%specification)) THEN
1724  CALL flagerror("Problem specification is not allocated.",err,error,*999)
1725  ELSE IF(SIZE(problem%specification,1)<3) THEN
1726  CALL flagerror("Problem specification must have three entries for a bioelectric-finite elasticity problem.", &
1727  & err,error,*999)
1728  END IF
1729  SELECT CASE(problem%SPECIFICATION(2))
1730  CASE(problem_bioelectric_finite_elasticity_type)
1731  SELECT CASE(problem%SPECIFICATION(3))
1732 
1733  CASE(problem_gudunov_monodomain_simple_elasticity_subtype)
1734 
1735  control_loop_root=>problem%CONTROL_LOOP
1736  CALL control_loop_get(control_loop_root,control_loop_node,control_loop_parent,err,error,*999)
1737  !get the monodomain sub loop, solvers, solver, and finally geometric and field
1738  CALL control_loop_sub_loop_get(control_loop_parent,1,control_loop_monodomain,err,error,*999)
1739  CALL control_loop_solvers_get(control_loop_monodomain,solvers,err,error,*999)
1740  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
1741  solver_equations=>solver%SOLVER_EQUATIONS
1742  IF(ASSOCIATED(solver_equations)) THEN
1743  solver_mapping=>solver_equations%SOLVER_MAPPING
1744  IF(ASSOCIATED(solver_mapping)) THEN
1745  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
1746  IF(ASSOCIATED(equations_set)) THEN
1747  geometric_field_monodomain=>equations_set%GEOMETRY%GEOMETRIC_FIELD
1748  IF(.NOT.ASSOCIATED(geometric_field_monodomain)) THEN
1749  CALL flagerror("Geometric field is not associated.",err,error,*999)
1750  ENDIF
1751  ELSE
1752  CALL flagerror("Equations set is not associated.",err,error,*999)
1753  ENDIF
1754  ELSE
1755  CALL flagerror("Solver mapping is not associated.",err,error,*999)
1756  ENDIF
1757  ELSE
1758  CALL flagerror("Solver equations is not associated.",err,error,*999)
1759  ENDIF
1760  NULLIFY(solvers)
1761  NULLIFY(solver)
1762  NULLIFY(solver_mapping)
1763  NULLIFY(equations_set)
1764  NULLIFY(solver_equations)
1765  !get the finite elasticity sub loop, solvers, solver, and finally the dependent field
1766  CALL control_loop_sub_loop_get(control_loop_parent,2,control_loop_elasticity,err,error,*999)
1767  CALL control_loop_solvers_get(control_loop_elasticity,solvers,err,error,*999)
1768  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1769  solver_equations=>solver%SOLVER_EQUATIONS
1770  IF(ASSOCIATED(solver_equations)) THEN
1771  solver_mapping=>solver_equations%SOLVER_MAPPING
1772  IF(ASSOCIATED(solver_mapping)) THEN
1773  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
1774  IF(ASSOCIATED(equations_set)) THEN
1775  dependent_field_elasticity=>equations_set%DEPENDENT%DEPENDENT_FIELD
1776  IF(.NOT.ASSOCIATED(dependent_field_elasticity)) THEN
1777  CALL flagerror("Dependent field is not associated.",err,error,*999)
1778  ENDIF
1779  ELSE
1780  CALL flagerror("Equations set is not associated.",err,error,*999)
1781  ENDIF
1782  ELSE
1783  CALL flagerror("Solver mapping is not associated.",err,error,*999)
1784  ENDIF
1785  ELSE
1786  CALL flagerror("Solver equations is not associated.",err,error,*999)
1787  ENDIF
1788  DO component_idx=1,geometric_field_monodomain%VARIABLES(1)%NUMBER_OF_COMPONENTS
1789  !check for identical interpolation of the fields
1790  geometric_field_interpolation=geometric_field_monodomain%VARIABLES(1)%COMPONENTS(component_idx)%INTERPOLATION_TYPE
1791  dependent_field_interpolation=dependent_field_elasticity%VARIABLES(1)%COMPONENTS(component_idx)%INTERPOLATION_TYPE
1792  IF(geometric_field_interpolation==dependent_field_interpolation) THEN
1793  !copy the dependent field components to the geometric field
1794  CALL field_parameterstofieldparameterscopy(dependent_field_elasticity,field_u_variable_type, &
1795  & field_values_set_type,component_idx,geometric_field_monodomain,field_u_variable_type,field_values_set_type, &
1796  & component_idx,err,error,*999)
1797  ELSE
1798  local_error="The interpolation type of component number "//trim(number_to_vstring(component_idx,"*",err, &
1799  & error))//" of field number "//trim(number_to_vstring(geometric_field_monodomain%USER_NUMBER,"*",err, &
1800  & error))//" does not coincide with the interpolation type of field number " &
1801  & //trim(number_to_vstring(dependent_field_elasticity%USER_NUMBER,"*",err,error))//"."
1802  CALL flagerror(local_error,err,error,*999)
1803  ENDIF
1804  ENDDO
1805 
1806  CASE(problem_gudunov_monodomain_1d3d_elasticity_subtype,problem_monodomain_elasticity_w_titin_subtype, &
1807  & problem_monodomain_1d3d_active_strain_subtype)
1808 
1809  control_loop_root=>problem%CONTROL_LOOP
1810  CALL control_loop_get(control_loop_root,control_loop_node,control_loop_parent,err,error,*999)
1811  !get the monodomain sub loop, solvers, solver, and finally geometric field and dependent field
1812  CALL control_loop_sub_loop_get(control_loop_parent,1,control_loop_monodomain,err,error,*999)
1813  CALL control_loop_solvers_get(control_loop_monodomain,solvers,err,error,*999)
1814  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
1815  solver_equations=>solver%SOLVER_EQUATIONS
1816  IF(ASSOCIATED(solver_equations)) THEN
1817  solver_mapping=>solver_equations%SOLVER_MAPPING
1818  IF(ASSOCIATED(solver_mapping)) THEN
1819  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
1820  IF(ASSOCIATED(equations_set)) THEN
1821  geometric_field_monodomain=>equations_set%GEOMETRY%GEOMETRIC_FIELD
1822  IF(.NOT.ASSOCIATED(geometric_field_monodomain)) THEN
1823  CALL flagerror("Geometric field is not associated.",err,error,*999)
1824  ENDIF
1825  ! the Field_V_Variable_Type contains the 3D nodal positions
1826  dependent_field_monodomain=>equations_set%DEPENDENT%DEPENDENT_FIELD
1827  IF(.NOT.ASSOCIATED(dependent_field_monodomain)) THEN
1828  CALL flagerror("Dependent field is not associated.",err,error,*999)
1829  ENDIF
1830  independent_field_monodomain=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
1831  IF(.NOT.ASSOCIATED(independent_field_monodomain)) THEN
1832  CALL flagerror("Independent field is not associated.",err,error,*999)
1833  ENDIF
1834  ELSE
1835  CALL flagerror("Equations set is not associated.",err,error,*999)
1836  ENDIF
1837  ELSE
1838  CALL flagerror("Solver mapping is not associated.",err,error,*999)
1839  ENDIF
1840  ELSE
1841  CALL flagerror("Solver equations is not associated.",err,error,*999)
1842  ENDIF
1843  NULLIFY(solvers)
1844  NULLIFY(solver)
1845  NULLIFY(solver_mapping)
1846  NULLIFY(equations_set)
1847  NULLIFY(solver_equations)
1848  !get the finite elasticity sub loop, solvers, solver, and finally the dependent and independent fields
1849  CALL control_loop_sub_loop_get(control_loop_parent,2,control_loop_elasticity,err,error,*999)
1850  CALL control_loop_solvers_get(control_loop_elasticity,solvers,err,error,*999)
1851  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
1852  solver_equations=>solver%SOLVER_EQUATIONS
1853  IF(ASSOCIATED(solver_equations)) THEN
1854  solver_mapping=>solver_equations%SOLVER_MAPPING
1855  IF(ASSOCIATED(solver_mapping)) THEN
1856  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
1857  IF(ASSOCIATED(equations_set)) THEN
1858  independent_field_elasticity=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
1859  IF(.NOT.ASSOCIATED(independent_field_elasticity)) THEN
1860  CALL flagerror("Independent field is not associated.",err,error,*999)
1861  ENDIF
1862  geometric_field_elasticity=>equations_set%GEOMETRY%GEOMETRIC_FIELD
1863  IF(.NOT.ASSOCIATED(geometric_field_elasticity)) THEN
1864  CALL flagerror("Dependent field is not associated.",err,error,*999)
1865  ENDIF
1866  ELSE
1867  CALL flagerror("Equations set is not associated.",err,error,*999)
1868  ENDIF
1869  ELSE
1870  CALL flagerror("Solver mapping is not associated.",err,error,*999)
1871  ENDIF
1872  ELSE
1873  CALL flagerror("Solver equations is not associated.",err,error,*999)
1874  ENDIF
1875 
1876 
1877  node_idx=0
1878  node_idx_2=0
1879  fibre_idx=0
1880  CALL field_variable_get(dependent_field_monodomain,field_v_variable_type,field_var_dep_m,err,error,*999)
1881  CALL field_variable_get(geometric_field_monodomain,field_u_variable_type,field_var_geo_m,err,error,*999)
1882  CALL field_variable_get(independent_field_elasticity,field_v_variable_type,field_var_ind_fe,err,error,*999)
1883  CALL field_variable_get(independent_field_monodomain,field_u1_variable_type,field_var_ind_m,err,error,*999)
1884  CALL field_variable_get(independent_field_monodomain,field_u2_variable_type,field_var_ind_m_2,err,error,*999)
1885 
1886  nodes_mapping=>geometric_field_monodomain%DECOMPOSITION%DOMAIN(geometric_field_monodomain%DECOMPOSITION% &
1887  & mesh_component_number)%PTR%MAPPINGS%NODES
1888 
1889  elements_topology=>geometric_field_elasticity%DECOMPOSITION%TOPOLOGY%ELEMENTS
1890 
1891 
1892  !get the maximum contraction velocity
1893  dof_idx=field_var_ind_m_2%COMPONENTS(2)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
1894  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u2_variable_type, &
1895  & field_values_set_type,dof_idx,velocity_max,err,error,*999)
1896  !NOTE: VELOCITY_MAX is the max shortening velocity, and hence negative!!!
1897  !The max lengthening velocity is assumed to be abs(VELOCITY_MAX)/2.0
1898 
1899  !get the time step of the elasticity problem
1900  time_step=control_loop_parent%TIME_LOOP%TIME_INCREMENT
1901 
1902 
1903  !loop over the elements of the finite elasticity mesh (internal and boundary elements)
1904  !no need to consider ghost elements here since only bioelectrical fields are changed
1905  DO element_idx=1,elements_topology%NUMBER_OF_ELEMENTS
1906  ne=elements_topology%ELEMENTS(element_idx)%LOCAL_NUMBER
1907  my_element_idx=element_idx
1908 
1909  !the Field_V_Variable_Type of the FE independent field contains the number of nodes in each Xi-direction of the bioelectrics grid
1910  dof_idx=field_var_ind_fe%COMPONENTS(1)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
1911  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type,field_values_set_type, &
1912  & dof_idx,nodes_in_xi_1,err,error,*999)
1913  dof_idx=field_var_ind_fe%COMPONENTS(2)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
1914  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type,field_values_set_type, &
1915  & dof_idx,nodes_in_xi_2,err,error,*999)
1916  dof_idx=field_var_ind_fe%COMPONENTS(3)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
1917  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type,field_values_set_type, &
1918  & dof_idx,nodes_in_xi_3,err,error,*999)
1919  !beginning of a fibre in this element: 1=yes, 0=no
1920  dof_idx=field_var_ind_fe%COMPONENTS(4)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
1921  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type,field_values_set_type, &
1922  & dof_idx,start_elem,err,error,*999)
1923 
1924  !if there is no bioelectrics grid in this finite elasticity element, or the fibres don't begin in this element, jump to the next element
1925  IF((nodes_in_xi_1==0).OR.(nodes_in_xi_2==0).OR.(nodes_in_xi_3==0).OR.(start_elem==0)) cycle
1926 
1927  start_element=ne
1928  start_element_idx=my_element_idx
1929 
1930  !assume Xi(1) to be normal to the seed surface, i.e. the seed points have Xi(1)=0
1931  xi=[0.0_dp,1.0_dp/(REAL(2*nodes_in_xi_2)),1.0_dp/(REAL(2*nodes_in_xi_3))]
1932 
1933  !assume that the bioelectrics node numbers are increased in order Xi(1), Xi(2), Xi(3)
1934  DO n3=1,nodes_in_xi_3
1935  DO n2=1,nodes_in_xi_2
1936  fibre_idx=fibre_idx+1
1937 
1938  !loop over the FE elements that contain nodes of the very same fibres
1939  DO
1940  !get the finite elasticity dependent field interpolation parameters of this element
1941  interpolation_parameters=>equations_set%EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS &
1942  & (field_u_variable_type)%PTR
1943  CALL field_interpolation_parameters_element_get(field_values_set_type,ne,interpolation_parameters, &
1944  & err,error,*999)
1945  interpolated_point=>equations_set%EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_POINT(field_u_variable_type)%PTR
1946 
1947  !get the positions of the Gauss points of the Finite Elasticity element, GAUSS_POSITIONS(components,number_of_Gauss_points)
1948  gauss_positions=>geometric_field_elasticity%DECOMPOSITION%DOMAIN(geometric_field_elasticity%DECOMPOSITION% &
1949  & mesh_component_number)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(ne)%BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP( &
1950  & basis_default_quadrature_scheme)%PTR%GAUSS_POSITIONS
1951 
1952  DO n1=1,nodes_in_xi_1
1953  node_idx=node_idx+1
1954 
1955  !store the fibre number this bioelectrics node belongs to.
1956  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_v_variable_type, &
1957  & field_values_set_type,1,1,node_idx,3,fibre_idx,err,error,*999)
1958 
1959  !find the interpolated position of the bioelectric grid node from the FE dependent field
1960  CALL field_interpolate_xi(no_part_deriv,xi,interpolated_point,err,error,*999)
1961  !update the bioelectrics dependent field Field_V_Variable_Type
1962  !the Field_V_Variable_Type of the monodomain dependent field contains the nodal positions in 3D
1963  CALL field_parameter_set_update_local_node(dependent_field_monodomain,field_v_variable_type, &
1964  & field_values_set_type,1,1,node_idx,1,interpolated_point%VALUES(1,1),err,error,*999)
1965  CALL field_parameter_set_update_local_node(dependent_field_monodomain,field_v_variable_type, &
1966  & field_values_set_type,1,1,node_idx,2,interpolated_point%VALUES(2,1),err,error,*999)
1967  CALL field_parameter_set_update_local_node(dependent_field_monodomain,field_v_variable_type, &
1968  & field_values_set_type,1,1,node_idx,3,interpolated_point%VALUES(3,1),err,error,*999)
1969 
1970  IF((n1==1).AND.(ne==start_element)) THEN
1971  !a new line of bioelectrics grid nodes begins
1972  CALL field_parameter_set_update_local_node(geometric_field_monodomain,field_u_variable_type, &
1973  & field_values_set_type,1,1,node_idx,1,0.0_dp,err,error,*999)
1974  ELSE
1975  !get the position in 3D of the previous node
1976  dof_idx2=field_var_dep_m%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx-1)% &
1977  & derivatives(1)%VERSIONS(1)
1978  CALL field_parameter_set_get_local_dof(dependent_field_monodomain,field_v_variable_type, &
1979  & field_values_set_type,dof_idx2,previous_node(1),err,error,*999)
1980  dof_idx2=field_var_dep_m%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx-1)% &
1981  & derivatives(1)%VERSIONS(1)
1982  CALL field_parameter_set_get_local_dof(dependent_field_monodomain,field_v_variable_type, &
1983  & field_values_set_type,dof_idx2,previous_node(2),err,error,*999)
1984  dof_idx2=field_var_dep_m%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx-1)% &
1985  & derivatives(1)%VERSIONS(1)
1986  CALL field_parameter_set_get_local_dof(dependent_field_monodomain,field_v_variable_type, &
1987  & field_values_set_type,dof_idx2,previous_node(3),err,error,*999)
1988 
1989  !compute the distance between the previous node and the actual node
1990  dist=sqrt( &
1991  & (interpolated_point%VALUES(1,1)-previous_node(1))*(interpolated_point%VALUES(1,1)-previous_node(1))+ &
1992  & (interpolated_point%VALUES(2,1)-previous_node(2))*(interpolated_point%VALUES(2,1)-previous_node(2))+ &
1993  & (interpolated_point%VALUES(3,1)-previous_node(3))*(interpolated_point%VALUES(3,1)-previous_node(3)))
1994 
1995 
1996  !CONTRACTION VELOCITY CALCULATION
1997 
1998  !get the distance between the 2 nodes in the previous time step
1999  dof_idx=field_var_ind_m_2%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)% &
2000  & derivatives(1)%VERSIONS(1)
2001  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u2_variable_type, &
2002  & field_values_set_type,dof_idx,old_dist,err,error,*999)
2003 
2004  !get the distance between the 2 nodes before two time step
2005  dof_idx=field_var_ind_m_2%COMPONENTS(4)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)% &
2006  & derivatives(1)%VERSIONS(1)
2007  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u2_variable_type, &
2008  & field_values_set_type,dof_idx,old_dist_2,err,error,*999)
2009 
2010  !get the distance between the 2 nodes before three time step
2011  dof_idx=field_var_ind_m_2%COMPONENTS(5)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)% &
2012  & derivatives(1)%VERSIONS(1)
2013  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u2_variable_type, &
2014  & field_values_set_type,dof_idx,old_dist_3,err,error,*999)
2015 
2016  !get the distance between the 2 nodes before four time step
2017  dof_idx=field_var_ind_m_2%COMPONENTS(6)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)% &
2018  & derivatives(1)%VERSIONS(1)
2019  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u2_variable_type, &
2020  & field_values_set_type,dof_idx,old_dist_4,err,error,*999)
2021 
2022  !compute the new contraction velocity
2023  !VELOCITY=(DIST-OLD_DIST)/TIME_STEP
2024  velocity=0.25_dp*((dist-old_dist)/time_step+(dist-old_dist_2)/(2.0_dp*time_step)+ &
2025  & (dist-old_dist_3)/(3.0_dp*time_step)+(dist-old_dist_4)/(4.0_dp*time_step))
2026  IF(.NOT. calc_closest_gauss_point) THEN
2027  !NOTE: VELOCITY_MAX is the max shortening velocity, and hence negative!!!
2028  IF(velocity<velocity_max) THEN
2029  CALL flag_warning('Exceeded maximum contraction velocity (shortening).',err,error,*999)
2030  velocity=velocity_max
2031  !The max lengthening velocity is assumed to be VELOCITY_MAX/2.0
2032  ELSEIF(velocity>(abs(velocity_max)/2.0_dp)) THEN
2033  CALL flag_warning('Exceeded maximum contraction velocity (lengthening).',err,error,*999)
2034  velocity=-velocity_max/2.0_dp
2035  ENDIF
2036  ENDIF
2037 
2038  !store the relative contraction velocity in component 3 of the U2 variable of the monodomain independent field
2039  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_u2_variable_type, &
2040  & field_values_set_type,1,1,node_idx,3,velocity/abs(velocity_max),err,error,*999)
2041 
2042  !store the node distance for contraction velocity calculation
2043  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_u2_variable_type, &
2044  & field_values_set_type,1,1,node_idx,1,dist,err,error,*999)
2045 
2046  !store the node distance for contraction velocity calculation
2047  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_u2_variable_type, &
2048  & field_values_set_type,1,1,node_idx,4,old_dist,err,error,*999)
2049 
2050  !store the node distance for contraction velocity calculation
2051  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_u2_variable_type, &
2052  & field_values_set_type,1,1,node_idx,5,old_dist_2,err,error,*999)
2053 
2054  !store the node distance for contraction velocity calculation
2055  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_u2_variable_type, &
2056  & field_values_set_type,1,1,node_idx,6,old_dist_3,err,error,*999)
2057 
2058 
2059 
2060  !get the position in 1D of the previous node
2061  dof_idx2=field_var_geo_m%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx-1)% &
2062  & derivatives(1)%VERSIONS(1)
2063  CALL field_parameter_set_get_local_dof(geometric_field_monodomain,field_u_variable_type, &
2064  & field_values_set_type,dof_idx2,value_left,err,error,*999)
2065  !update the current 1D node position
2066  CALL field_parameter_set_update_local_node(geometric_field_monodomain,field_u_variable_type, &
2067  & field_values_set_type,1,1,node_idx,1,value_left+dist,err,error,*999)
2068 
2069  !get the initial sarcomere half length and initial node distance
2070  dof_idx=field_var_ind_m%COMPONENTS(2)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
2071  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u1_variable_type, &
2072  & field_values_set_type,dof_idx,sarco_length_init,err,error,*999)
2073  dof_idx=field_var_ind_m%COMPONENTS(3)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
2074  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u1_variable_type, &
2075  & field_values_set_type,dof_idx,dist_init,err,error,*999)
2076  !update the current sarcomere half length
2077  VALUE=sarco_length_init*dist/dist_init
2078  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_u1_variable_type, &
2079  & field_values_set_type,1,1,node_idx,1,VALUE,err,error,*999)
2080 
2081  !update the first node to the same value as the second node (no better info available)
2082  IF((n1==2).AND.(ne==start_element)) THEN
2083  !current sarcomere half length
2084  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_u1_variable_type, &
2085  & field_values_set_type,1,1,node_idx-1,1,VALUE,err,error,*999)
2086  !old node distance
2087  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_u2_variable_type, &
2088  & field_values_set_type,1,1,node_idx-1,1,dist,err,error,*999)
2089  !relative contraction velocity
2090  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_u2_variable_type, &
2091  & field_values_set_type,1,1,node_idx-1,3,velocity/abs(velocity_max),err,error,*999)
2092  ENDIF
2093 
2094  ENDIF !((n1==1).AND.(ne==START_ELEMENT))
2095 
2096  IF(calc_closest_gauss_point) THEN
2097  !calculate the closest finite elasticity Gauss point of each bioelectrics node
2098  distance=1000000.0_dp
2099  gauss_point=0
2100  DO gauss_idx=1,SIZE(gauss_positions,2)
2101  !compute the distance between the bioelectrics node and the Gauss point
2102  VALUE=sqrt( &
2103  & (xi(1)-gauss_positions(1,gauss_idx))*(xi(1)-gauss_positions(1,gauss_idx))+ &
2104  & (xi(2)-gauss_positions(2,gauss_idx))*(xi(2)-gauss_positions(2,gauss_idx))+ &
2105  & (xi(3)-gauss_positions(3,gauss_idx))*(xi(3)-gauss_positions(3,gauss_idx)))
2106  IF(VALUE<distance) THEN
2107  distance=VALUE
2108  gauss_point=gauss_idx
2109  ENDIF
2110  ENDDO !gauss_idx
2111  IF(gauss_point==0) CALL flag_warning("Closest Gauss Point not found",err,error,*999)
2112  !store the nearest Gauss Point info and the inElement info (local element number!!!)
2113  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_v_variable_type, &
2114  & field_values_set_type,1,1,node_idx,4,gauss_point,err,error,*999)
2115  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_v_variable_type, &
2116  & field_values_set_type,1,1,node_idx,5,ne,err,error,*999)
2117  ENDIF !CALC_CLOSEST_GAUSS_POINT
2118 
2119  IF(start_elem==1) THEN
2120  !fibres start in this element
2121  xi(1)=xi(1)+1.0_dp/(REAL(nodes_in_xi_1-1))
2122  ELSEIF(start_elem==0) THEN
2123  !fibres don't start in this element
2124  xi(1)=xi(1)+1.0_dp/(REAL(nodes_in_xi_1))
2125  ELSE
2126  local_error="The start element index is incorrect. The index is "// &
2127  & trim(number_to_vstring(start_elem,"*",err,error))// &
2128  & " and should be zero or one."
2129  CALL flagerror(local_error,err,error,*999)
2130  ENDIF
2131 
2132  ENDDO !n1
2133 
2134 
2135 
2136 
2137 !tomo new
2138  !smooth of the velocity field
2139 
2140  !arithmetic mean of all rel_velo values within one FE element
2141 ! VELOCITY=0.0_DP
2142 ! DO n1=1,nodes_in_Xi_1
2143 ! node_idx_2=node_idx_2+1
2144 ! dof_idx=FIELD_VAR_IND_M_2%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx_2)% &
2145 ! & DERIVATIVES(1)%VERSIONS(1)
2146 ! CALL FIELD_PARAMETER_SET_GET_LOCAL_DOF(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2147 ! & FIELD_VALUES_SET_TYPE,dof_idx,VALUE,ERR,ERROR,*999)
2148 ! VELOCITY=VELOCITY+VALUE
2149 ! ENDDO
2150 ! VELOCITY=VELOCITY/nodes_in_Xi_1
2151 
2152 ! node_idx_2=node_idx_2-nodes_in_Xi_1
2153 ! DO n1=1,nodes_in_Xi_1
2154 ! node_idx_2=node_idx_2+1
2155 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_NODE(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2156 ! & FIELD_VALUES_SET_TYPE,1,1,node_idx_2,3,VELOCITY,ERR,ERROR,*999)
2157 ! ENDDO
2158 
2159 !--------------------------------------------------------------------------
2160 
2161 ! !moving average
2162 ! offset=3
2163 !
2164 ! !do the first three nodes of a fibre manually - arithmetic mean
2165 ! VELOCITY=0.0_DP
2166 
2167 ! node_idx_2=node_idx_2+1
2168 !
2169 ! dof_idx=FIELD_VAR_IND_M_2%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx_2)% &
2170 ! & DERIVATIVES(1)%VERSIONS(1)
2171 ! CALL FIELD_PARAMETER_SET_GET_LOCAL_DOF(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2172 ! & FIELD_VALUES_SET_TYPE,dof_idx,VALUE,ERR,ERROR,*999)
2173 ! VELOCITY=VELOCITY+VALUE
2174 
2175 ! node_idx_2=node_idx_2+1
2176 ! dof_idx=FIELD_VAR_IND_M_2%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx_2)% &
2177 ! & DERIVATIVES(1)%VERSIONS(1)
2178 ! CALL FIELD_PARAMETER_SET_GET_LOCAL_DOF(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2179 ! & FIELD_VALUES_SET_TYPE,dof_idx,VALUE,ERR,ERROR,*999)
2180 ! VELOCITY=VELOCITY+VALUE
2181 
2182 ! node_idx_2=node_idx_2+1
2183 ! dof_idx=FIELD_VAR_IND_M_2%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx_2)% &
2184 ! & DERIVATIVES(1)%VERSIONS(1)
2185 ! CALL FIELD_PARAMETER_SET_GET_LOCAL_DOF(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2186 ! & FIELD_VALUES_SET_TYPE,dof_idx,VALUE,ERR,ERROR,*999)
2187 ! VELOCITY=VELOCITY+VALUE
2188 !
2189 ! VELOCITY=VELOCITY/offset
2190 
2191 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_NODE(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2192 ! & FIELD_VALUES_SET_TYPE,1,1,node_idx_2-2,3,VELOCITY,ERR,ERROR,*999)
2193 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_NODE(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2194 ! & FIELD_VALUES_SET_TYPE,1,1,node_idx_2-1,3,VELOCITY,ERR,ERROR,*999)
2195 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_NODE(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2196 ! & FIELD_VALUES_SET_TYPE,1,1,node_idx_2,3,VELOCITY,ERR,ERROR,*999)
2197 
2198 ! !do the major part as moving average
2199 ! DO n1=1+offset,nodes_in_Xi_1-offset
2200 ! node_idx_2=node_idx_2+1
2201 ! VELOCITY=0.0_DP
2202 ! DO n4=node_idx_2-offset,node_idx_2+offset
2203 ! dof_idx=FIELD_VAR_IND_M_2%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(n4)% &
2204 ! & DERIVATIVES(1)%VERSIONS(1)
2205 ! CALL FIELD_PARAMETER_SET_GET_LOCAL_DOF(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2206 ! & FIELD_VALUES_SET_TYPE,dof_idx,VALUE,ERR,ERROR,*999)
2207 ! VELOCITY=VELOCITY+VALUE
2208 ! ENDDO !n4
2209 ! VELOCITY=VELOCITY/(2*offset+1)
2210 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_NODE(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2211 ! & FIELD_VALUES_SET_TYPE,1,1,node_idx_2,3,VELOCITY,ERR,ERROR,*999)
2212 ! ENDDO !n1
2213 !
2214 ! !do the last three nodes of a fibre manually - arithmetic mean
2215 ! VELOCITY=0.0_DP
2216 
2217 ! node_idx_2=node_idx_2+1
2218 ! dof_idx=FIELD_VAR_IND_M_2%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx_2)% &
2219 ! & DERIVATIVES(1)%VERSIONS(1)
2220 ! CALL FIELD_PARAMETER_SET_GET_LOCAL_DOF(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2221 ! & FIELD_VALUES_SET_TYPE,dof_idx,VALUE,ERR,ERROR,*999)
2222 ! VELOCITY=VELOCITY+VALUE
2223 
2224 ! node_idx_2=node_idx_2+1
2225 ! dof_idx=FIELD_VAR_IND_M_2%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx_2)% &
2226 ! & DERIVATIVES(1)%VERSIONS(1)
2227 ! CALL FIELD_PARAMETER_SET_GET_LOCAL_DOF(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2228 ! & FIELD_VALUES_SET_TYPE,dof_idx,VALUE,ERR,ERROR,*999)
2229 ! VELOCITY=VELOCITY+VALUE
2230 
2231 ! node_idx_2=node_idx_2+1
2232 ! dof_idx=FIELD_VAR_IND_M_2%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx_2)% &
2233 ! & DERIVATIVES(1)%VERSIONS(1)
2234 ! CALL FIELD_PARAMETER_SET_GET_LOCAL_DOF(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2235 ! & FIELD_VALUES_SET_TYPE,dof_idx,VALUE,ERR,ERROR,*999)
2236 ! VELOCITY=VELOCITY+VALUE
2237 !
2238 ! VELOCITY=VELOCITY/offset
2239 
2240 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_NODE(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2241 ! & FIELD_VALUES_SET_TYPE,1,1,node_idx_2-2,3,VELOCITY,ERR,ERROR,*999)
2242 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_NODE(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2243 ! & FIELD_VALUES_SET_TYPE,1,1,node_idx_2-1,3,VELOCITY,ERR,ERROR,*999)
2244 ! CALL FIELD_PARAMETER_SET_UPDATE_LOCAL_NODE(INDEPENDENT_FIELD_MONODOMAIN,FIELD_U2_VARIABLE_TYPE, &
2245 ! & FIELD_VALUES_SET_TYPE,1,1,node_idx_2,3,VELOCITY,ERR,ERROR,*999)
2246 !!tomo new end
2247 
2248 
2249  !if there is not an adjacent element in positive XI_1 direction, go to the next FE element
2250  IF(elements_topology%ELEMENTS(my_element_idx)%ADJACENT_ELEMENTS(1)%NUMBER_OF_ADJACENT_ELEMENTS==0) EXIT
2251 
2252  !consider the adjacent element in positive XI_1 direction
2253  ne=elements_topology%ELEMENTS(my_element_idx)%ADJACENT_ELEMENTS(1)%ADJACENT_ELEMENTS(1)
2254 
2255  !if a fibre starts in the next element, go to the next FE elem
2256  dof_idx=field_var_ind_fe%COMPONENTS(4)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
2257  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type, &
2258  & field_values_set_type,dof_idx,start_elem,err,error,*999)
2259  !beginning of a fibre in this element: 1=yes, 0=no
2260  IF (start_elem==1) THEN
2261  ne=elements_topology%ELEMENTS(element_idx)%LOCAL_NUMBER
2262  EXIT
2263  ENDIF
2264 
2265  !find the element_idx that corresponds to ne
2266  my_element_idx=0
2267  DO idx=1,elements_topology%NUMBER_OF_ELEMENTS
2268  IF(ne==elements_topology%ELEMENTS(idx)%ADJACENT_ELEMENTS(0)%ADJACENT_ELEMENTS(1)) THEN
2269  my_element_idx=idx
2270  EXIT
2271  ENDIF
2272  ENDDO
2273  IF(my_element_idx==0) CALL flagerror("my_element_idx not found.",err,error,*999)
2274 
2275  dof_idx=field_var_ind_fe%COMPONENTS(1)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
2276  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type, &
2277  & field_values_set_type,dof_idx,nodes_in_xi_1,err,error,*999)
2278 
2279  start_elem=0 !fibres don't start in this element
2280 
2281  xi(1)=1.0_dp/(REAL(nodes_in_xi_1))
2282 
2283  ENDDO !
2284  !for the beginning of the next fibre, go back to the element in which the last fibre started
2285  ne=start_element
2286 
2287  dof_idx=field_var_ind_fe%COMPONENTS(1)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
2288  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type, &
2289  & field_values_set_type,dof_idx,nodes_in_xi_1,err,error,*999)
2290 
2291  my_element_idx=start_element_idx
2292  start_elem=1 !fibres start in this element
2293  xi(1)=0.0_dp
2294  xi(2)=xi(2)+1.0_dp/(REAL(nodes_in_xi_2))
2295  ENDDO !n2
2296  xi(1)=0.0_dp
2297  xi(2)=1.0_dp/(REAL(2*nodes_in_xi_2))
2298  xi(3)=xi(3)+1.0_dp/(REAL(nodes_in_xi_3))
2299  ENDDO !n3
2300 
2301  ENDDO !element_idx
2302 
2303  CASE(problem_monodomain_elasticity_velocity_subtype)
2304 
2305  control_loop_root=>problem%CONTROL_LOOP
2306  CALL control_loop_get(control_loop_root,control_loop_node,control_loop_parent,err,error,*999)
2307  !get the monodomain sub loop, solvers, solver, and finally geometric field and dependent field
2308  CALL control_loop_sub_loop_get(control_loop_parent,1,control_loop_monodomain,err,error,*999)
2309  CALL control_loop_solvers_get(control_loop_monodomain,solvers,err,error,*999)
2310  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
2311  solver_equations=>solver%SOLVER_EQUATIONS
2312  IF(ASSOCIATED(solver_equations)) THEN
2313  solver_mapping=>solver_equations%SOLVER_MAPPING
2314  IF(ASSOCIATED(solver_mapping)) THEN
2315  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
2316  IF(ASSOCIATED(equations_set)) THEN
2317  geometric_field_monodomain=>equations_set%GEOMETRY%GEOMETRIC_FIELD
2318  IF(.NOT.ASSOCIATED(geometric_field_monodomain)) THEN
2319  CALL flagerror("Geometric field is not associated.",err,error,*999)
2320  ENDIF
2321  ! the Field_V_Variable_Type contains the 3D nodal positions
2322  dependent_field_monodomain=>equations_set%DEPENDENT%DEPENDENT_FIELD
2323  IF(.NOT.ASSOCIATED(dependent_field_monodomain)) THEN
2324  CALL flagerror("Dependent field is not associated.",err,error,*999)
2325  ENDIF
2326  independent_field_monodomain=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
2327  IF(.NOT.ASSOCIATED(independent_field_monodomain)) THEN
2328  CALL flagerror("Independent field is not associated.",err,error,*999)
2329  ENDIF
2330  ELSE
2331  CALL flagerror("Equations set is not associated.",err,error,*999)
2332  ENDIF
2333  ELSE
2334  CALL flagerror("Solver mapping is not associated.",err,error,*999)
2335  ENDIF
2336  ELSE
2337  CALL flagerror("Solver equations is not associated.",err,error,*999)
2338  ENDIF
2339  NULLIFY(solvers)
2340  NULLIFY(solver)
2341  NULLIFY(solver_mapping)
2342  NULLIFY(equations_set)
2343  NULLIFY(solver_equations)
2344  !get the finite elasticity sub loop, solvers, solver, and finally the dependent and independent fields
2345  CALL control_loop_sub_loop_get(control_loop_parent,2,control_loop_elasticity,err,error,*999)
2346  CALL control_loop_solvers_get(control_loop_elasticity,solvers,err,error,*999)
2347  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
2348  solver_equations=>solver%SOLVER_EQUATIONS
2349  IF(ASSOCIATED(solver_equations)) THEN
2350  solver_mapping=>solver_equations%SOLVER_MAPPING
2351  IF(ASSOCIATED(solver_mapping)) THEN
2352  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
2353  IF(ASSOCIATED(equations_set)) THEN
2354  independent_field_elasticity=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
2355  IF(.NOT.ASSOCIATED(independent_field_elasticity)) THEN
2356  CALL flagerror("Independent field is not associated.",err,error,*999)
2357  ENDIF
2358  geometric_field_elasticity=>equations_set%GEOMETRY%GEOMETRIC_FIELD
2359  IF(.NOT.ASSOCIATED(geometric_field_elasticity)) THEN
2360  CALL flagerror("Dependent field is not associated.",err,error,*999)
2361  ENDIF
2362  ELSE
2363  CALL flagerror("Equations set is not associated.",err,error,*999)
2364  ENDIF
2365  ELSE
2366  CALL flagerror("Solver mapping is not associated.",err,error,*999)
2367  ENDIF
2368  ELSE
2369  CALL flagerror("Solver equations is not associated.",err,error,*999)
2370  ENDIF
2371 
2372 
2373  node_idx=0
2374  node_idx_2=0
2375  fibre_idx=0
2376  CALL field_variable_get(dependent_field_monodomain,field_v_variable_type,field_var_dep_m,err,error,*999)
2377  CALL field_variable_get(geometric_field_monodomain,field_u_variable_type,field_var_geo_m,err,error,*999)
2378  CALL field_variable_get(independent_field_elasticity,field_v_variable_type,field_var_ind_fe,err,error,*999)
2379 
2380  nodes_mapping=>geometric_field_monodomain%DECOMPOSITION%DOMAIN(geometric_field_monodomain%DECOMPOSITION% &
2381  & mesh_component_number)%PTR%MAPPINGS%NODES
2382 
2383  elements_topology=>geometric_field_elasticity%DECOMPOSITION%TOPOLOGY%ELEMENTS
2384 
2385  !loop over the elements of the finite elasticity mesh (internal and boundary elements)
2386  !no need to consider ghost elements here since only bioelectrical fields are changed
2387  DO element_idx=1,elements_topology%NUMBER_OF_ELEMENTS
2388  ne=elements_topology%ELEMENTS(element_idx)%LOCAL_NUMBER
2389  my_element_idx=element_idx
2390 
2391  !the Field_V_Variable_Type of the FE independent field contains the number of nodes in each Xi-direction of the bioelectrics grid
2392  dof_idx=field_var_ind_fe%COMPONENTS(1)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
2393  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type,field_values_set_type, &
2394  & dof_idx,nodes_in_xi_1,err,error,*999)
2395  dof_idx=field_var_ind_fe%COMPONENTS(2)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
2396  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type,field_values_set_type, &
2397  & dof_idx,nodes_in_xi_2,err,error,*999)
2398  dof_idx=field_var_ind_fe%COMPONENTS(3)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
2399  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type,field_values_set_type, &
2400  & dof_idx,nodes_in_xi_3,err,error,*999)
2401  !beginning of a fibre in this element: 1=yes, 0=no
2402  dof_idx=field_var_ind_fe%COMPONENTS(4)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
2403  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type,field_values_set_type, &
2404  & dof_idx,start_elem,err,error,*999)
2405 
2406  !if there is no bioelectrics grid in this finite elasticity element, or the fibres don't begin in this element, jump to the next element
2407  IF((nodes_in_xi_1==0).OR.(nodes_in_xi_2==0).OR.(nodes_in_xi_3==0).OR.(start_elem==0)) cycle
2408 
2409  start_element=ne
2410  start_element_idx=my_element_idx
2411 
2412  !assume Xi(1) to be normal to the seed surface, i.e. the seed points have Xi(1)=0
2413  xi=[0.0_dp,1.0_dp/(REAL(2*nodes_in_xi_2)),1.0_dp/(REAL(2*nodes_in_xi_3))]
2414 
2415  !assume that the bioelectrics node numbers are increased in order Xi(1), Xi(2), Xi(3)
2416  DO n3=1,nodes_in_xi_3
2417  DO n2=1,nodes_in_xi_2
2418  fibre_idx=fibre_idx+1
2419 
2420  !loop over the FE elements that contain nodes of the very same fibres
2421  DO
2422  !get the finite elasticity dependent field interpolation parameters of this element
2423  interpolation_parameters=>equations_set%EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS &
2424  & (field_u_variable_type)%PTR
2425  CALL field_interpolation_parameters_element_get(field_values_set_type,ne,interpolation_parameters, &
2426  & err,error,*999)
2427  interpolated_point=>equations_set%EQUATIONS%INTERPOLATION%DEPENDENT_INTERP_POINT(field_u_variable_type)%PTR
2428 
2429  !get the positions of the Gauss points of the Finite Elasticity element, GAUSS_POSITIONS(components,number_of_Gauss_points)
2430  gauss_positions=>geometric_field_elasticity%DECOMPOSITION%DOMAIN(geometric_field_elasticity%DECOMPOSITION% &
2431  & mesh_component_number)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(ne)%BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP( &
2432  & basis_default_quadrature_scheme)%PTR%GAUSS_POSITIONS
2433 
2434  DO n1=1,nodes_in_xi_1
2435  node_idx=node_idx+1
2436 
2437  !store the fibre number this bioelectrics node belongs to.
2438  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_v_variable_type, &
2439  & field_values_set_type,1,1,node_idx,3,fibre_idx,err,error,*999)
2440 
2441  !find the interpolated position of the bioelectric grid node from the FE dependent field
2442  CALL field_interpolate_xi(no_part_deriv,xi,interpolated_point,err,error,*999)
2443  !update the bioelectrics dependent field Field_V_Variable_Type
2444  !the Field_V_Variable_Type of the monodomain dependent field contains the nodal positions in 3D
2445  CALL field_parameter_set_update_local_node(dependent_field_monodomain,field_v_variable_type, &
2446  & field_values_set_type,1,1,node_idx,1,interpolated_point%VALUES(1,1),err,error,*999)
2447  CALL field_parameter_set_update_local_node(dependent_field_monodomain,field_v_variable_type, &
2448  & field_values_set_type,1,1,node_idx,2,interpolated_point%VALUES(2,1),err,error,*999)
2449  CALL field_parameter_set_update_local_node(dependent_field_monodomain,field_v_variable_type, &
2450  & field_values_set_type,1,1,node_idx,3,interpolated_point%VALUES(3,1),err,error,*999)
2451 
2452  IF((n1==1).AND.(ne==start_element)) THEN
2453  !a new line of bioelectrics grid nodes begins
2454  CALL field_parameter_set_update_local_node(geometric_field_monodomain,field_u_variable_type, &
2455  & field_values_set_type,1,1,node_idx,1,0.0_dp,err,error,*999)
2456  ELSE
2457  !get the position in 3D of the previous node
2458  dof_idx2=field_var_dep_m%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx-1)% &
2459  & derivatives(1)%VERSIONS(1)
2460  CALL field_parameter_set_get_local_dof(dependent_field_monodomain,field_v_variable_type, &
2461  & field_values_set_type,dof_idx2,previous_node(1),err,error,*999)
2462  dof_idx2=field_var_dep_m%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx-1)% &
2463  & derivatives(1)%VERSIONS(1)
2464  CALL field_parameter_set_get_local_dof(dependent_field_monodomain,field_v_variable_type, &
2465  & field_values_set_type,dof_idx2,previous_node(2),err,error,*999)
2466  dof_idx2=field_var_dep_m%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx-1)% &
2467  & derivatives(1)%VERSIONS(1)
2468  CALL field_parameter_set_get_local_dof(dependent_field_monodomain,field_v_variable_type, &
2469  & field_values_set_type,dof_idx2,previous_node(3),err,error,*999)
2470 
2471  !compute the distance between the previous node and the actual node
2472  dist=sqrt( &
2473  & (interpolated_point%VALUES(1,1)-previous_node(1))*(interpolated_point%VALUES(1,1)-previous_node(1))+ &
2474  & (interpolated_point%VALUES(2,1)-previous_node(2))*(interpolated_point%VALUES(2,1)-previous_node(2))+ &
2475  & (interpolated_point%VALUES(3,1)-previous_node(3))*(interpolated_point%VALUES(3,1)-previous_node(3)))
2476 
2477  !get the position in 1D of the previous node
2478  dof_idx2=field_var_geo_m%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx-1)% &
2479  & derivatives(1)%VERSIONS(1)
2480  CALL field_parameter_set_get_local_dof(geometric_field_monodomain,field_u_variable_type, &
2481  & field_values_set_type,dof_idx2,value_left,err,error,*999)
2482  !update the current 1D node position
2483  CALL field_parameter_set_update_local_node(geometric_field_monodomain,field_u_variable_type, &
2484  & field_values_set_type,1,1,node_idx,1,value_left+dist,err,error,*999)
2485 
2486  ENDIF !((n1==1).AND.(ne==START_ELEMENT))
2487 
2488  IF(calc_closest_gauss_point) THEN
2489  !calculate the closest finite elasticity Gauss point of each bioelectrics node
2490  distance=1000000.0_dp
2491  gauss_point=0
2492  DO gauss_idx=1,SIZE(gauss_positions,2)
2493  !compute the distance between the bioelectrics node and the Gauss point
2494  VALUE=sqrt( &
2495  & (xi(1)-gauss_positions(1,gauss_idx))*(xi(1)-gauss_positions(1,gauss_idx))+ &
2496  & (xi(2)-gauss_positions(2,gauss_idx))*(xi(2)-gauss_positions(2,gauss_idx))+ &
2497  & (xi(3)-gauss_positions(3,gauss_idx))*(xi(3)-gauss_positions(3,gauss_idx)))
2498  IF(VALUE<distance) THEN
2499  distance=VALUE
2500  gauss_point=gauss_idx
2501  ENDIF
2502  ENDDO !gauss_idx
2503  IF(gauss_point==0) CALL flag_warning("Closest Gauss Point not found",err,error,*999)
2504  !store the nearest Gauss Point info and the inElement info (local element number!!!)
2505  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_v_variable_type, &
2506  & field_values_set_type,1,1,node_idx,4,gauss_point,err,error,*999)
2507  CALL field_parameter_set_update_local_node(independent_field_monodomain,field_v_variable_type, &
2508  & field_values_set_type,1,1,node_idx,5,ne,err,error,*999)
2509  ENDIF !CALC_CLOSEST_GAUSS_POINT
2510 
2511  IF(start_elem==1) THEN
2512  !fibres start in this element
2513  xi(1)=xi(1)+1.0_dp/(REAL(nodes_in_xi_1-1))
2514  ELSEIF(start_elem==0) THEN
2515  !fibres don't start in this element
2516  xi(1)=xi(1)+1.0_dp/(REAL(nodes_in_xi_1))
2517  ELSE
2518  local_error="The start element index is incorrect. The index is "// &
2519  & trim(number_to_vstring(start_elem,"*",err,error))// &
2520  & " and should be zero or one."
2521  CALL flagerror(local_error,err,error,*999)
2522  ENDIF
2523 
2524  ENDDO !n1
2525 
2526 
2527  !if there is not an adjacent element in positive XI_1 direction, go to the next FE element
2528  IF(elements_topology%ELEMENTS(my_element_idx)%ADJACENT_ELEMENTS(1)%NUMBER_OF_ADJACENT_ELEMENTS==0) EXIT
2529 
2530  !consider the adjacent element in positive XI_1 direction
2531  ne=elements_topology%ELEMENTS(my_element_idx)%ADJACENT_ELEMENTS(1)%ADJACENT_ELEMENTS(1)
2532 
2533  !if a fibre starts in the next element, go to the next FE elem
2534  dof_idx=field_var_ind_fe%COMPONENTS(4)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
2535  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type, &
2536  & field_values_set_type,dof_idx,start_elem,err,error,*999)
2537  !beginning of a fibre in this element: 1=yes, 0=no
2538  IF (start_elem==1) THEN
2539  ne=elements_topology%ELEMENTS(element_idx)%LOCAL_NUMBER
2540  EXIT
2541  ENDIF
2542 
2543  !find the element_idx that corresponds to ne
2544  my_element_idx=0
2545  DO idx=1,elements_topology%NUMBER_OF_ELEMENTS
2546  IF(ne==elements_topology%ELEMENTS(idx)%ADJACENT_ELEMENTS(0)%ADJACENT_ELEMENTS(1)) THEN
2547  my_element_idx=idx
2548  EXIT
2549  ENDIF
2550  ENDDO
2551  IF(my_element_idx==0) CALL flagerror("my_element_idx not found.",err,error,*999)
2552 
2553  dof_idx=field_var_ind_fe%COMPONENTS(1)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
2554  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type, &
2555  & field_values_set_type,dof_idx,nodes_in_xi_1,err,error,*999)
2556 
2557  start_elem=0 !fibres don't start in this element
2558 
2559  xi(1)=1.0_dp/(REAL(nodes_in_xi_1))
2560 
2561  ENDDO !
2562  !for the beginning of the next fibre, go back to the element in which the last fibre started
2563  ne=start_element
2564 
2565  dof_idx=field_var_ind_fe%COMPONENTS(1)%PARAM_TO_DOF_MAP%ELEMENT_PARAM2DOF_MAP%ELEMENTS(ne)
2566  CALL field_parameter_set_get_local_dof(independent_field_elasticity,field_v_variable_type, &
2567  & field_values_set_type,dof_idx,nodes_in_xi_1,err,error,*999)
2568 
2569  my_element_idx=start_element_idx
2570  start_elem=1 !fibres start in this element
2571  xi(1)=0.0_dp
2572  xi(2)=xi(2)+1.0_dp/(REAL(nodes_in_xi_2))
2573  ENDDO !n2
2574  xi(1)=0.0_dp
2575  xi(2)=1.0_dp/(REAL(2*nodes_in_xi_2))
2576  xi(3)=xi(3)+1.0_dp/(REAL(nodes_in_xi_3))
2577  ENDDO !n3
2578 
2579  ENDDO !element_idx
2580 
2581  CASE DEFAULT
2582  local_error="The third problem specification of "// &
2583  & trim(number_to_vstring(problem%specification(2),"*",err,error))// &
2584  & " is not valid for a bioelectrics finite elasticity of a multi physics problem."
2585  CALL flagerror(local_error,err,error,*999)
2586  END SELECT
2587  CASE DEFAULT
2588  local_error="The second problem specification of "// &
2589  & trim(number_to_vstring(problem%specification(2),"*",err,error))// &
2590  & " is not valid for a multi physics problem."
2591  CALL flagerror(local_error,err,error,*999)
2592  END SELECT
2593  ELSE
2594  CALL flagerror("Control loop problem is not associated.",err,error,*999)
2595  ENDIF
2596  ELSE
2597  !the main time loop - do nothing!
2598  ENDIF
2599  ELSE
2600  CALL flagerror("Control loop is not associated.",err,error,*999)
2601  ENDIF
2602 
2603  exits("BioelectricFiniteElasticity_UpdateGeometricField")
2604  RETURN
2605 999 errors("BioelectricFiniteElasticity_UpdateGeometricField",err,error)
2606  exits("BioelectricFiniteElasticity_UpdateGeometricField")
2607  RETURN 1
2608 
2609  END SUBROUTINE bioelectricfiniteelasticity_updategeometricfield
2610 
2611  !
2612  !================================================================================================================================
2613  !
2614 
2617  SUBROUTINE bioelectricfiniteelasticity_independentfieldinterpolate(CONTROL_LOOP,ERR,ERROR,*)
2618 
2619  !Argument variables
2620  TYPE(control_loop_type), POINTER :: CONTROL_LOOP
2621  INTEGER(INTG), INTENT(OUT) :: ERR
2622  TYPE(varying_string), INTENT(OUT) :: ERROR
2623  !Local Variables
2624  TYPE(control_loop_type), POINTER :: CONTROL_LOOP_ROOT,CONTROL_LOOP_PARENT,CONTROL_LOOP_ELASTICITY,CONTROL_LOOP_MONODOMAIN
2625  TYPE(problem_type), POINTER :: PROBLEM
2626  TYPE(solvers_type), POINTER :: SOLVERS
2627  TYPE(solver_type), POINTER :: SOLVER
2628  TYPE(field_type), POINTER :: INDEPENDENT_FIELD_MONODOMAIN,INDEPENDENT_FIELD_ELASTICITY
2629  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
2630  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
2631  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
2632  TYPE(domain_mapping_type), POINTER :: ELEMENTS_MAPPING,NODES_MAPPING
2633  TYPE(varying_string) :: LOCAL_ERROR
2634  TYPE(field_variable_type), POINTER :: FIELD_VARIABLE_U,FIELD_VARIABLE_V,FIELD_VARIABLE_FE
2635  INTEGER(INTG) :: node_idx,element_idx,gauss_idx,ne
2636  INTEGER(INTG) :: nearestGP,inElement,dof_idx
2637  INTEGER(INTG) :: NUMBER_OF_GAUSS_POINTS
2638  REAL(DP) :: ACTIVE_STRESS
2639  REAL(DP) :: TITIN_STRESS_UNBOUND,TITIN_STRESS_BOUND,TITIN_STRESS_CROSS_FIBRE_UNBOUND,TITIN_STRESS_CROSS_FIBRE_BOUND,ACTIVATION
2640  INTEGER(INTG), PARAMETER :: MAX_NUMBER_OF_GAUSS_POINTS=64
2641  INTEGER(INTG) :: NUMBER_OF_NODES(max_number_of_gauss_points)
2642  REAL(DP) :: ACTIVE_STRESS_VALUES(max_number_of_gauss_points)
2643  REAL(DP) :: TITIN_STRESS_VALUES_UNBOUND(max_number_of_gauss_points),TITIN_STRESS_VALUES_BOUND(max_number_of_gauss_points)
2644  REAL(DP) :: TITIN_STRESS_VALUES_CROSS_FIBRE_UNBOUND(max_number_of_gauss_points)
2645  REAL(DP) :: TITIN_STRESS_VALUES_CROSS_FIBRE_BOUND(max_number_of_gauss_points)
2646  REAL(DP) :: ACTIVATION_VALUES(max_number_of_gauss_points)
2647  REAL(DP) :: A_1,A_2,x_1,x_2
2648  REAL(DP) :: A_1_VALUES(max_number_of_gauss_points),A_2_VALUES(max_number_of_gauss_points), &
2649  & x_1_VALUES(MAX_NUMBER_OF_GAUSS_POINTS),x_2_VALUES(MAX_NUMBER_OF_GAUSS_POINTS)
2650 
2651  NULLIFY(control_loop_parent)
2652  NULLIFY(control_loop_monodomain)
2653  NULLIFY(control_loop_elasticity)
2654  NULLIFY(solvers)
2655  NULLIFY(solver)
2656  NULLIFY(field_variable_u)
2657  NULLIFY(field_variable_v)
2658  NULLIFY(field_variable_fe)
2659 
2660  enters("BioelectricFiniteElasticity_IndependentFieldInterpolate",err,error,*999)
2661 
2662  IF(ASSOCIATED(control_loop)) THEN
2663  problem=>control_loop%PROBLEM
2664  IF(ASSOCIATED(problem)) THEN
2665  IF(.NOT.ALLOCATED(problem%specification)) THEN
2666  CALL flagerror("Problem specification is not allocated.",err,error,*999)
2667  ELSE IF(SIZE(problem%specification,1)<3) THEN
2668  CALL flagerror("Problem specification must have three entries for a bioelectric-finite elasticity problem.", &
2669  & err,error,*999)
2670  END IF
2671  SELECT CASE(problem%SPECIFICATION(3))
2672  CASE(problem_gudunov_monodomain_1d3d_elasticity_subtype,problem_monodomain_elasticity_w_titin_subtype, &
2673  & problem_monodomain_elasticity_velocity_subtype,problem_monodomain_1d3d_active_strain_subtype)
2674  IF(control_loop%NUMBER_OF_SUB_LOOPS==0) THEN
2675  control_loop_root=>problem%CONTROL_LOOP
2676  CALL control_loop_get(control_loop_root,control_loop_node,control_loop_parent,err,error,*999)
2677  !--- MONODOMAIN ---
2678  CALL control_loop_sub_loop_get(control_loop_parent,1,control_loop_monodomain,err,error,*999)
2679  CALL control_loop_solvers_get(control_loop_monodomain,solvers,err,error,*999)
2680  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
2681  solver_equations=>solver%SOLVER_EQUATIONS
2682  IF(ASSOCIATED(solver_equations)) THEN
2683  solver_mapping=>solver_equations%SOLVER_MAPPING
2684  IF(ASSOCIATED(solver_mapping)) THEN
2685  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
2686  IF(ASSOCIATED(equations_set)) THEN
2687  independent_field_monodomain=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
2688  IF(.NOT.ASSOCIATED(independent_field_monodomain)) CALL flagerror("Independent field is not associated.", &
2689  & err,error,*999)
2690  ELSE
2691  CALL flagerror("Equations set is not associated.",err,error,*999)
2692  ENDIF
2693  ELSE
2694  CALL flagerror("Solver mapping is not associated.",err,error,*999)
2695  ENDIF
2696  ELSE
2697  CALL flagerror("Solver equations is not associated.",err,error,*999)
2698  ENDIF
2699 
2700  !--- FINITE ELASTICITY ---
2701  NULLIFY(solvers)
2702  NULLIFY(solver)
2703  CALL control_loop_sub_loop_get(control_loop_parent,2,control_loop_elasticity,err,error,*999)
2704  CALL control_loop_solvers_get(control_loop_elasticity,solvers,err,error,*999)
2705  CALL solvers_solver_get(solvers,1,solver,err,error,*999)
2706  solver_equations=>solver%SOLVER_EQUATIONS
2707  IF(ASSOCIATED(solver_equations)) THEN
2708  solver_mapping=>solver_equations%SOLVER_MAPPING
2709  IF(ASSOCIATED(solver_mapping)) THEN
2710  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
2711  IF(ASSOCIATED(equations_set)) THEN
2712  independent_field_elasticity=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
2713  IF(.NOT.ASSOCIATED(independent_field_elasticity)) CALL flagerror("Independent field is not associated.",err, &
2714  & error,*999)
2715  ELSE
2716  CALL flagerror("Equations set is not associated.",err,error,*999)
2717  ENDIF
2718  ELSE
2719  CALL flagerror("Solver mapping is not associated.",err,error,*999)
2720  ENDIF
2721  ELSE
2722  CALL flagerror("Solver equations is not associated.",err,error,*999)
2723  ENDIF
2724 
2725  !--- NOW INTERPOLATE ---
2726  elements_mapping=>independent_field_elasticity%DECOMPOSITION%DOMAIN(independent_field_elasticity%DECOMPOSITION% &
2727  & mesh_component_number)%PTR%MAPPINGS%ELEMENTS
2728  nodes_mapping=>independent_field_monodomain%DECOMPOSITION%DOMAIN(independent_field_monodomain%DECOMPOSITION% &
2729  & mesh_component_number)%PTR%MAPPINGS%NODES
2730 
2731  CALL field_variable_get(independent_field_monodomain,field_u_variable_type,field_variable_u,err,error,*999)
2732  CALL field_variable_get(independent_field_monodomain,field_v_variable_type,field_variable_v,err,error,*999)
2733  CALL field_variable_get(independent_field_elasticity,field_u_variable_type,field_variable_fe,err,error,*999)
2734 
2735  !loop over the finite elasticity elements
2736  !first process the internal and boundary elements
2737  DO element_idx=elements_mapping%INTERNAL_START,elements_mapping%BOUNDARY_FINISH
2738  ne=elements_mapping%DOMAIN_LIST(element_idx)
2739 
2740  number_of_gauss_points=independent_field_elasticity%DECOMPOSITION%DOMAIN(independent_field_elasticity% &
2741  & decomposition%MESH_COMPONENT_NUMBER)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(ne)%BASIS%QUADRATURE%QUADRATURE_SCHEME_MAP &
2742  & (basis_default_quadrature_scheme)%PTR%NUMBER_OF_GAUSS
2743 
2744  IF(number_of_gauss_points>max_number_of_gauss_points) CALL flagerror( &
2745  & "NUMBER_OF_GAUSS_POINTS is greater than MAX_NUMBER_OF_GAUSS_POINTS.",err,error,*999)
2746  number_of_nodes=0
2747  active_stress_values=0.0_dp
2748  titin_stress_values_unbound=0.0_dp
2749  titin_stress_values_bound=0.0_dp
2750  titin_stress_values_cross_fibre_unbound=0.0_dp
2751  titin_stress_values_cross_fibre_bound=0.0_dp
2752  activation_values=0.0_dp
2753  a_1_values=0.0_dp
2754  a_2_values=0.0_dp
2755  x_1_values=0.0_dp
2756  x_2_values=0.0_dp
2757 
2758  !loop over the bioelectrics nodes
2759  DO node_idx=1,nodes_mapping%NUMBER_OF_LOCAL
2760  dof_idx=field_variable_v%COMPONENTS(5)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2761  & versions(1)
2762  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_v_variable_type,field_values_set_type, &
2763  & dof_idx,inelement,err,error,*999) !component 5 of variable V contains inElem info (LOCAL NUMBERING!!!)
2764 
2765  !check if the bioelectrics node is located within the finite elasticity element
2766  IF(inelement==ne) THEN
2767  !component 4 of variable V contains Nearest Gauss Point info
2768  dof_idx=field_variable_v%COMPONENTS(4)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2769  & versions(1)
2770  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_v_variable_type,field_values_set_type, &
2771  & dof_idx,nearestgp,err,error,*999)
2772  IF(nearestgp>max_number_of_gauss_points) CALL flagerror( &
2773  & "Nearest Gauss Point is greater than MAX_NUMBER_OF_GAUSS_POINTS.",err,error,*999)
2774 
2775  SELECT CASE(problem%SPECIFICATION(3))
2776  CASE(problem_gudunov_monodomain_1d3d_elasticity_subtype,problem_monodomain_elasticity_velocity_subtype)
2777 
2778  !component 1 of variable U contains the active stress
2779  dof_idx=field_variable_u%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2780  & versions(1)
2781  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2782  & field_values_set_type,dof_idx,active_stress,err,error,*999)
2783 
2784  !count the number of bioelectrics nodes that are closest to each finite elasticity Gauss point
2785  number_of_nodes(nearestgp)=number_of_nodes(nearestgp)+1
2786  !add up the active stress value
2787  active_stress_values(nearestgp)=active_stress_values(nearestgp)+active_stress
2788 
2789  CASE(problem_monodomain_elasticity_w_titin_subtype)
2790 
2791  !component 1 of variable U contains the active stress
2792  dof_idx=field_variable_u%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2793  & versions(1)
2794  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2795  & field_values_set_type,dof_idx,active_stress,err,error,*999)
2796 
2797  !count the number of bioelectrics nodes that are closest to each finite elasticity Gauss point
2798  number_of_nodes(nearestgp)=number_of_nodes(nearestgp)+1
2799  !add up the active stress value
2800  active_stress_values(nearestgp)=active_stress_values(nearestgp)+active_stress
2801 
2802  !component 2 of variable U contains the titin stress unbound
2803  dof_idx=field_variable_u%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2804  & versions(1)
2805  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2806  & field_values_set_type,dof_idx,titin_stress_unbound,err,error,*999)
2807  !component 3 of variable U contains the titin stress bound
2808  dof_idx=field_variable_u%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2809  & versions(1)
2810  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2811  & field_values_set_type,dof_idx,titin_stress_bound,err,error,*999)
2812  !component 4 of variable U contains the titin XF-stress (cross-fibre directions) unbound
2813  dof_idx=field_variable_u%COMPONENTS(4)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2814  & versions(1)
2815  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2816  & field_values_set_type,dof_idx,titin_stress_cross_fibre_unbound,err,error,*999)
2817  !component 5 of variable U contains the titin XF-stress (cross-fibre directions) bound
2818  dof_idx=field_variable_u%COMPONENTS(5)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2819  & versions(1)
2820  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2821  & field_values_set_type,dof_idx,titin_stress_cross_fibre_bound,err,error,*999)
2822  !component 6 of variable U contains the titin activation
2823  dof_idx=field_variable_u%COMPONENTS(6)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2824  & versions(1)
2825  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2826  & field_values_set_type,dof_idx,activation,err,error,*999)
2827 
2828  titin_stress_values_unbound(nearestgp)=titin_stress_values_unbound(nearestgp)+titin_stress_unbound
2829  titin_stress_values_bound(nearestgp)=titin_stress_values_bound(nearestgp)+titin_stress_bound
2830  titin_stress_values_cross_fibre_unbound(nearestgp)=titin_stress_values_cross_fibre_unbound(nearestgp) + &
2831  & titin_stress_cross_fibre_unbound
2832  titin_stress_values_cross_fibre_bound(nearestgp)=titin_stress_values_cross_fibre_bound(nearestgp) + &
2833  & titin_stress_cross_fibre_bound
2834  activation_values(nearestgp)=activation_values(nearestgp)+activation
2835 
2836  CASE(problem_monodomain_1d3d_active_strain_subtype)
2837 
2838  !count the number of bioelectrics nodes that are closest to each finite elasticity Gauss point
2839  number_of_nodes(nearestgp)=number_of_nodes(nearestgp)+1
2840 
2841  !component 1 of variable U contains A_1
2842  dof_idx=field_variable_u%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2843  & versions(1)
2844  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2845  & field_values_set_type,dof_idx,a_1,err,error,*999)
2846  !component 2 of variable U contains A_2
2847  dof_idx=field_variable_u%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2848  & versions(1)
2849  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2850  & field_values_set_type,dof_idx,a_2,err,error,*999)
2851  !component 3 of variable U contains x_1
2852  dof_idx=field_variable_u%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2853  & versions(1)
2854  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2855  & field_values_set_type,dof_idx,x_1,err,error,*999)
2856  !component 4 of variable U contains x_2
2857  dof_idx=field_variable_u%COMPONENTS(4)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
2858  & versions(1)
2859  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u_variable_type, &
2860  & field_values_set_type,dof_idx,x_2,err,error,*999)
2861 
2862  a_1_values(nearestgp)=a_1_values(nearestgp)+a_1
2863  a_2_values(nearestgp)=a_2_values(nearestgp)+a_2
2864  x_1_values(nearestgp)=x_1_values(nearestgp)+x_1
2865  x_2_values(nearestgp)=x_2_values(nearestgp)+x_2
2866  END SELECT
2867 
2868  ENDIF
2869  ENDDO
2870 
2871  !loop over the finite elasticity Gauss points
2872  DO gauss_idx=1,number_of_gauss_points
2873  !make sure we don't divide by zero
2874  IF(number_of_nodes(gauss_idx)<=0) THEN
2875  active_stress=0.0_dp
2876  titin_stress_unbound=0.0_dp
2877  titin_stress_bound=0.0_dp
2878  titin_stress_cross_fibre_unbound=0.0_dp
2879  titin_stress_cross_fibre_bound=0.0_dp
2880  activation=0.0_dp
2881  a_1=0.0_dp
2882  a_2=0.0_dp
2883  x_1=0.0_dp
2884  x_2=0.0_dp
2885  ELSE
2886  active_stress=active_stress_values(gauss_idx)/number_of_nodes(gauss_idx)
2887  titin_stress_unbound=titin_stress_values_unbound(gauss_idx)/number_of_nodes(gauss_idx)
2888  titin_stress_bound=titin_stress_values_bound(gauss_idx)/number_of_nodes(gauss_idx)
2889  titin_stress_cross_fibre_unbound=titin_stress_values_cross_fibre_unbound(gauss_idx)/number_of_nodes(gauss_idx)
2890  titin_stress_cross_fibre_bound=titin_stress_values_cross_fibre_bound(gauss_idx)/ &
2891  & number_of_nodes(gauss_idx)
2892  activation=activation_values(gauss_idx)/number_of_nodes(gauss_idx)
2893  a_1=a_1_values(gauss_idx)/number_of_nodes(gauss_idx)
2894  a_2=a_2_values(gauss_idx)/number_of_nodes(gauss_idx)
2895  x_1=x_1_values(gauss_idx)/number_of_nodes(gauss_idx)
2896  x_2=x_2_values(gauss_idx)/number_of_nodes(gauss_idx)
2897  ENDIF
2898 
2899  SELECT CASE(problem%SPECIFICATION(3))
2900  CASE(problem_gudunov_monodomain_1d3d_elasticity_subtype,problem_monodomain_elasticity_velocity_subtype)
2901 
2902  dof_idx=field_variable_fe%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2903  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2904  & field_values_set_type,dof_idx,active_stress,err,error,*999)
2905 
2906  CASE(problem_monodomain_elasticity_w_titin_subtype)
2907 
2908  dof_idx=field_variable_fe%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2909  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2910  & field_values_set_type,dof_idx,active_stress,err,error,*999)
2911  dof_idx=field_variable_fe%COMPONENTS(2)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2912  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2913  & field_values_set_type,dof_idx,titin_stress_unbound,err,error,*999)
2914  dof_idx=field_variable_fe%COMPONENTS(3)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2915  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2916  & field_values_set_type,dof_idx,titin_stress_bound,err,error,*999)
2917  dof_idx=field_variable_fe%COMPONENTS(4)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2918  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2919  & field_values_set_type,dof_idx,titin_stress_cross_fibre_unbound,err,error,*999)
2920  dof_idx=field_variable_fe%COMPONENTS(5)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2921  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2922  & field_values_set_type,dof_idx,titin_stress_cross_fibre_bound,err,error,*999)
2923  dof_idx=field_variable_fe%COMPONENTS(6)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2924  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2925  & field_values_set_type,dof_idx,activation,err,error,*999)
2926 
2927  CASE(problem_monodomain_1d3d_active_strain_subtype)
2928 
2929  dof_idx=field_variable_fe%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2930  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2931  & field_values_set_type,dof_idx,a_1,err,error,*999)
2932  dof_idx=field_variable_fe%COMPONENTS(2)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2933  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2934  & field_values_set_type,dof_idx,a_2,err,error,*999)
2935  dof_idx=field_variable_fe%COMPONENTS(3)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2936  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2937  & field_values_set_type,dof_idx,x_1,err,error,*999)
2938  dof_idx=field_variable_fe%COMPONENTS(4)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_idx,ne)
2939  CALL field_parameter_set_update_local_dof(independent_field_elasticity,field_u_variable_type, &
2940  & field_values_set_type,dof_idx,x_2,err,error,*999)
2941 
2942  END SELECT
2943 
2944  ENDDO !gauss_idx
2945  ENDDO !element_idx
2946 
2947  !now the ghost elements -- get the relevant info from the other computational nodes
2948  CALL field_parameter_set_update_start(independent_field_elasticity, &
2949  & field_u_variable_type,field_values_set_type,err,error,*999)
2950  CALL field_parameter_set_update_finish(independent_field_elasticity, &
2951  & field_u_variable_type,field_values_set_type,err,error,*999)
2952 
2953  ENDIF
2954  CASE DEFAULT
2955  local_error="Independent field interpolation is not implemented for problem subtype " &
2956  & //trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))
2957  CALL flagerror(local_error,err,error,*999)
2958  END SELECT
2959  ELSE
2960  CALL flagerror("Problem is not associated.",err,error,*999)
2961  ENDIF
2962  ELSE
2963  CALL flagerror("Control loop is not associated.",err,error,*999)
2964  ENDIF
2965 
2966  exits("BioelectricFiniteElasticity_IndependentFieldInterpolate")
2967  RETURN
2968 999 errors("BioelectricFiniteElasticity_IndependentFieldInterpolate",err,error)
2969  exits("BioelectricFiniteElasticity_IndependentFieldInterpolate")
2970  RETURN 1
2971 
2972  END SUBROUTINE bioelectricfiniteelasticity_independentfieldinterpolate
2973 
2974  !
2975  !================================================================================================================================
2976  !
2977 
2980  SUBROUTINE bioelectric_finite_elasticity_compute_titin(CONTROL_LOOP,ERR,ERROR,*)
2981 
2982  !Argument variables
2983  TYPE(control_loop_type), POINTER :: CONTROL_LOOP !A pointer to the time control loop
2984  INTEGER(INTG), INTENT(OUT) :: ERR !The error code
2985  TYPE(varying_string), INTENT(OUT) :: ERROR !The error string
2986  !Local Variables
2987  TYPE(control_loop_type), POINTER :: CONTROL_LOOP_ROOT,CONTROL_LOOP_PARENT,CONTROL_LOOP_MONODOMAIN
2988  TYPE(problem_type), POINTER :: PROBLEM
2989  TYPE(solvers_type), POINTER :: SOLVERS
2990  TYPE(solver_type), POINTER :: SOLVER
2991  TYPE(field_type), POINTER :: INDEPENDENT_FIELD_MONODOMAIN
2992  TYPE(field_variable_type), POINTER :: FIELD_VAR_IND_M
2993  TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
2994  TYPE(solver_mapping_type), POINTER :: SOLVER_MAPPING
2995  TYPE(equations_set_type), POINTER :: EQUATIONS_SET
2996  TYPE(domain_mapping_type), POINTER :: NODES_MAPPING
2997  INTEGER(INTG) :: node_idx,dof_idx
2998  REAL(DP), PARAMETER :: LENGTH_ACTIN=1.04_dp,length_mband=0.0625_dp,length_myosin=0.7375_dp
2999  REAL(DP), PARAMETER :: LENGTH_ZERO=0.635_dp,length_zdisc=0.05_dp
3000  REAL(DP) :: F0,FORCE,TITIN_BOUND,TITIN_XF_BOUND,TITIN_UNBOUND,TITIN_XF_UNBOUND
3001  REAL(DP) :: ELONGATION,ELONGATION_NEW
3002  REAL(DP) :: ELONGATION_DIST_IG,ELONGATION_PEVK
3003  REAL(DP) :: LENGTH_TITIN,LENGTH_INIT_TITIN,LENGTH_DIST_IG_F0,LENGTH_DIST_IG
3004  REAL(DP) :: STIFFNESS_PEVK
3005  REAL(DP) :: SARCO_LENGTH_AT_ACTIVATION,SARCO_LENGTH
3006  REAL(DP) :: ACTIN_MYOSIN_DISTANCE,d10
3007  REAL(DP) :: SLOPE,STIFFNESS_DIST,FORCE_DISTAL_IG,DELTA_F,DIFF_QUOT
3008  REAL(DP), PARAMETER, DIMENSION(5) :: COEFF_MATRIX=[5.0239_dp,-0.6717_dp,-2.5841_dp,-5.0128_dp,-5.0239_dp]
3009  REAL(DP), DIMENSION(250) :: LENGTHS_DIST_IG,FORCES_DIST_IG
3010  REAL(DP), PARAMETER :: DX=0.001_dp
3011  REAL(DP), PARAMETER :: FORCE_INCREMENT=1.e-5_dp,tol=1.e-5_dp
3012  INTEGER(INTG) :: INDEX_REF,INDEX_PSEUDO,INDEX_I
3013  INTEGER(INTG), PARAMETER :: DIM_DATA=250
3014  INTEGER(INTG) :: SWITCH_MODEL
3015 
3016  NULLIFY(control_loop_root)
3017  NULLIFY(control_loop_parent)
3018  NULLIFY(control_loop_monodomain)
3019  NULLIFY(problem)
3020  NULLIFY(solvers)
3021  NULLIFY(solver)
3022  NULLIFY(independent_field_monodomain)
3023  NULLIFY(solver_equations)
3024  NULLIFY(solver_mapping)
3025  NULLIFY(equations_set)
3026  NULLIFY(field_var_ind_m)
3027 
3028  enters("BIOELECTRIC_FINITE_ELASTICITY_COMPUTE_TITIN",err,error,*999)
3029 
3030  !the realtion between length and force of the distal_Ig region is very nonlinear. Linear interpolation of Rode's data is used here.
3031  forces_dist_ig= &
3032  & [0.0_dp,0.001_dp,0.002_dp,0.003_dp,0.004_dp,0.005_dp,0.006_dp,0.007_dp,0.008_dp,0.009_dp,0.01_dp,0.011_dp,0.012_dp, &
3033  & 0.013_dp,0.014_dp,0.015_dp,0.016_dp,0.017_dp,0.018_dp,0.019_dp,0.02_dp,0.021_dp,0.022_dp,0.023_dp,0.024_dp, &
3034  & 0.025_dp,0.026_dp,0.027_dp,0.028_dp,0.029_dp,0.03_dp,0.031_dp,0.032_dp,0.033_dp,0.034_dp,0.035_dp,0.036_dp, &
3035  & 0.037_dp,0.038_dp,0.039_dp,0.04_dp,0.041_dp,0.042_dp,0.043_dp,0.044_dp,0.045_dp,0.046_dp,0.047_dp,0.048_dp, &
3036  & 0.049_dp,0.05_dp,0.051_dp,0.052_dp,0.053_dp,0.054_dp,0.055_dp,0.056_dp,0.057_dp,0.058_dp,0.059_dp,0.06_dp, &
3037  & 0.061_dp,0.062_dp,0.063_dp,0.064_dp,0.065_dp,0.066_dp,0.067_dp,0.068_dp,0.069_dp,0.070_dp,0.071_dp,0.072_dp, &
3038  & 0.073_dp,0.074_dp,0.075_dp,0.076_dp,0.077_dp,0.078_dp,0.079_dp,0.080_dp,0.081_dp,0.082_dp,0.083_dp,0.084_dp, &
3039  & 0.085_dp,0.086_dp,0.087_dp,0.088_dp,0.089_dp,0.090_dp,0.091_dp,0.092_dp,0.093_dp,0.094_dp,0.095_dp,0.096_dp, &
3040  & 0.097_dp,0.098_dp,0.099_dp,0.1_dp,0.101_dp,0.102_dp,0.103_dp,0.104_dp,0.105_dp,0.106_dp,0.107_dp,0.108_dp, &
3041  & 0.109_dp,0.11_dp,0.111_dp,0.112_dp,0.113_dp,0.114_dp,0.115_dp,0.116_dp,0.117_dp,0.118_dp,0.119_dp,0.12_dp, &
3042  & 0.121_dp,0.122_dp,0.123_dp,0.124_dp,0.125_dp,0.126_dp,0.127_dp,0.128_dp,0.129_dp,0.13_dp,0.131_dp,0.132_dp, &
3043  & 0.133_dp,0.134_dp,0.135_dp,0.136_dp,0.137_dp,0.138_dp,0.139_dp,0.14_dp,0.141_dp,0.142_dp,0.143_dp,0.144_dp, &
3044  & 0.145_dp,0.146_dp,0.147_dp,0.148_dp,0.149_dp,0.15_dp,0.151_dp,0.152_dp,0.153_dp,0.154_dp,0.155_dp,0.156_dp, &
3045  & 0.157_dp,0.158_dp,0.159_dp,0.16_dp,0.161_dp,0.162_dp,0.163_dp,0.164_dp,0.165_dp,0.166_dp,0.167_dp, 0.168_dp, &
3046  & 0.169_dp,0.17_dp,0.171_dp,0.172_dp,0.173_dp,0.174_dp,0.175_dp,0.176_dp,0.177_dp,0.178_dp,0.179_dp,0.18_dp, &
3047  & 0.181_dp,0.182_dp,0.183_dp,0.184_dp,0.185_dp,0.186_dp,0.187_dp,0.188_dp,0.189_dp,0.19_dp,0.191_dp,0.192_dp, &
3048  & 0.193_dp,0.194_dp,0.195_dp,0.196_dp,0.197_dp,0.198_dp,0.199_dp,0.2_dp,0.201_dp,0.202_dp,0.203_dp,0.204_dp, &
3049  & 0.205_dp,0.206_dp,0.207_dp,0.208_dp,0.209_dp,0.21_dp,0.211_dp,0.212_dp,0.213_dp,0.214_dp,0.215_dp,0.216_dp, &
3050  & 0.217_dp,0.218_dp,0.219_dp,0.22_dp,0.221_dp,0.222_dp,0.223_dp,0.224_dp,0.225_dp,0.226_dp,0.227_dp,0.228_dp, &
3051  & 0.229_dp,0.23_dp,0.231_dp,0.232_dp,0.233_dp,0.234_dp,0.235_dp,0.236_dp,0.237_dp,0.238_dp,0.239_dp,0.24_dp, &
3052  & 0.241_dp,0.242_dp,0.243_dp,0.244_dp,0.245_dp,0.246_dp,0.247_dp,0.248_dp,0.249_dp]
3053  lengths_dist_ig= &
3054  & [0.0_dp,0.03461753545561_dp,0.049729169766010_dp,0.058506390323323_dp,0.064606296848594_dp,0.06922519775133_dp, &
3055  & 0.0729080120998386_dp,0.0759458896446241_dp,0.0785230355395668_dp,0.0807314335143191_dp,0.0826674161660979_dp, &
3056  & 0.0843819721302_dp,0.0859161360822_dp,0.087300738288_dp,0.0885510536196_dp,0.08970061165_dp,0.090751366_dp, &
3057  & 0.0917285714001_dp,0.0926271710799_dp,0.093467026018_dp,0.094254010845_dp,0.094992014919_dp,0.09568255451_dp, &
3058  & 0.0963346932312_dp,0.0969518155718_dp,0.097537000419_dp,0.098093047899_dp,0.098622503780_dp,0.09912768169_dp, &
3059  & 0.0996103583026_dp,0.1000734170008_dp,0.100517614158_dp,0.100942907967_dp,0.101351270601_dp,0.101745244913_dp, &
3060  & 0.1021260518375_dp,0.1024947934706_dp,0.102851281365_dp,0.103194317381_dp,0.103528086731_dp,0.103853341524_dp, &
3061  & 0.1041686065999_dp,0.1044739635997_dp,0.104772842609_dp,0.105064758806_dp,0.105347078318_dp,0.105624524436_dp, &
3062  & 0.1058959740402_dp,0.1061596374590_dp,0.106419674124_dp,0.106673222977_dp,0.106921779223_dp,0.107167251003_dp, &
3063  & 0.1074055929211_dp,0.1076418938850_dp,0.107872798106_dp,0.108100580266_dp,0.108324935479_dp,0.108545154704_dp, &
3064  & 0.1087634145225_dp,0.1089769308376_dp,0.109189523632_dp,0.109397109133_dp,0.109604335645_dp,0.109806785604_dp, &
3065  & 0.1100090293896_dp,0.1102069598668_dp,0.110404790711_dp,0.110598542577_dp,0.110792294444_dp,0.110982362315_dp, &
3066  & 0.1111722575399_dp,0.1113591719379_dp,0.111545514634_dp,0.111729654458_dp,0.111912731875_dp,0.112094428489_dp, &
3067  & 0.1122745119339_dp,0.1124540532647_dp,0.112631398979_dp,0.112808744694_dp,0.112983883284_dp,0.113158733268_dp, &
3068  & 0.1133324054841_dp,0.1135049882670_dp,0.113677360515_dp,0.113847891889_dp,0.114018423263_dp,0.114187784974_dp, &
3069  & 0.1143564686814_dp,0.1145249703398_dp,0.114691998719_dp,0.114859027099_dp,0.115025270473_dp,0.115190825075_dp, &
3070  & 0.1153563796784_dp,0.1155207617063_dp,0.115685013868_dp,0.115849024045_dp,0.116012135436_dp,0.116175246828_dp, &
3071  & 0.1163378963393_dp,0.1165000194746_dp,0.116662142610_dp,0.116823704015_dp,0.116984982740_dp,0.117146261465_dp, &
3072  & 0.1173069696799_dp,0.1174675396300_dp,0.117628109580_dp,0.117788165045_dp,0.117948154078_dp,0.118108143111_dp, &
3073  & 0.1182677147299_dp,0.1184272433357_dp,0.118586771941_dp,0.118745999791_dp,0.118905181478_dp,0.119064363164_dp, &
3074  & 0.1192233610196_dp,0.1193823026790_dp,0.119541244338_dp,0.119700101992_dp,0.119858904246_dp,0.120017706500_dp, &
3075  & 0.1201764919257_dp,0.1203352494533_dp,0.120494006981_dp,0.120652768322_dp,0.120811570168_dp,0.120970372014_dp, &
3076  & 0.1211291738603_dp,0.1212880693042_dp,0.121446999172_dp,0.121605929041_dp,0.121764923098_dp,0.121924059629_dp, &
3077  & 0.1220831961613_dp,0.1222423326927_dp,0.122401700265_dp,0.122561117298_dp,0.122720534331_dp,0.122880045624_dp, &
3078  & 0.1230398124447_dp,0.1231995792652_dp,0.123359346085_dp,0.123519381317_dp,0.123679562894_dp,0.123839744470_dp, &
3079  & 0.1239999260467_dp,0.1241605622478_dp,0.124321219453_dp,0.124481876659_dp,0.124642637543_dp,0.124803827368_dp, &
3080  & 0.1249650171945_dp,0.1251262070202_dp,0.125287609380_dp,0.125449385133_dp,0.125611160886_dp,0.125772936638_dp, &
3081  & 0.1259350043157_dp,0.1260974158090_dp,0.126259827302_dp,0.126422238795_dp,0.126584979901_dp,0.126748073636_dp, &
3082  & 0.1269111673704_dp,0.1270742611047_dp,0.127237669513_dp,0.127401488846_dp,0.127565308178_dp,0.127729127511_dp, &
3083  & 0.1278931842348_dp,0.1280577695422_dp,0.128222354849_dp,0.128386940157_dp,0.128551614623_dp,0.128717003455_dp, &
3084  & 0.1288823922862_dp,0.1290477811173_dp,0.129213169948_dp,0.129379259581_dp,0.129545486803_dp,0.129711714025_dp, &
3085  & 0.1298779412472_dp,0.1300445897140_dp,0.130211687650_dp,0.130378785586_dp,0.130545883522_dp,0.130713029866_dp, &
3086  & 0.1308810284274_dp,0.1310490269884_dp,0.131217025549_dp,0.131385024110_dp,0.131553528414_dp,0.131722455223_dp, &
3087  & 0.1318913820322_dp,0.1320603088411_dp,0.132229235650_dp,0.132399074312_dp,0.132568954822_dp,0.132738835332_dp, &
3088  & 0.1329087158420_dp,0.1330788764544_dp,0.133249734060_dp,0.133420591667_dp,0.133591449273_dp,0.133762306879_dp, &
3089  & 0.1339336992411_dp,0.1341055553883_dp,0.134277411535_dp,0.134449267682_dp,0.134621123829_dp,0.134793694483_dp, &
3090  & 0.1349665687658_dp,0.1351394430487_dp,0.135312317331_dp,0.135485191614_dp,0.135658878561_dp,0.135832788820_dp, &
3091  & 0.1360066990800_dp,0.1361806093395_dp,0.136354519599_dp,0.136529253243_dp,0.136704215658_dp,0.136879178072_dp, &
3092  & 0.1370541404878_dp,0.1372291029027_dp,0.137404806924_dp,0.137580836098_dp,0.137756865271_dp,0.137932894445_dp, &
3093  & 0.1381089236188_dp,0.1382855157880_dp,0.138462624830_dp,0.138639733872_dp,0.138816842915_dp,0.138993951957_dp, &
3094  & 0.1391713448843_dp,0.1393495454909_dp,0.139527746097_dp,0.139705946704_dp,0.139884147310_dp,0.140062347917_dp, &
3095  & 0.1402415516727_dp,0.1404208541990_dp,0.140600156725270_dp,0.140779459251523_dp,0.140958761777775_dp]
3096 
3097  IF(ASSOCIATED(control_loop)) THEN
3098  problem=>control_loop%PROBLEM
3099  IF(ASSOCIATED(problem)) THEN
3100  SELECT CASE(problem%SPECIFICATION(3))
3101  CASE(problem_monodomain_elasticity_w_titin_subtype)
3102  IF(control_loop%NUMBER_OF_SUB_LOOPS==0) THEN
3103  control_loop_root=>problem%CONTROL_LOOP
3104  CALL control_loop_get(control_loop_root,control_loop_node,control_loop_parent,err,error,*999)
3105  !The first control_loop is the one for monodomain
3106  CALL control_loop_sub_loop_get(control_loop_parent,1,control_loop_monodomain,err,error,*999)
3107  CALL control_loop_solvers_get(control_loop_monodomain,solvers,err,error,*999)
3108  !The second solver is associated with the diffusion part of the monodomain equation
3109  CALL solvers_solver_get(solvers,2,solver,err,error,*999)
3110  solver_equations=>solver%SOLVER_EQUATIONS
3111  IF(ASSOCIATED(solver_equations)) THEN
3112  solver_mapping=>solver_equations%SOLVER_MAPPING
3113  IF(ASSOCIATED(solver_mapping)) THEN
3114  equations_set=>solver_mapping%EQUATIONS_SETS(1)%PTR
3115  IF(ASSOCIATED(equations_set)) THEN
3116  independent_field_monodomain=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
3117  IF(.NOT.ASSOCIATED(independent_field_monodomain)) THEN
3118  CALL flagerror("Independent field is not associated.",err,error,*999)
3119  ENDIF
3120 
3121  CALL field_variable_get(independent_field_monodomain,field_u1_variable_type,field_var_ind_m,err,error,*999)
3122 
3123  nodes_mapping=>independent_field_monodomain%DECOMPOSITION%DOMAIN(independent_field_monodomain%DECOMPOSITION% &
3124  & mesh_component_number)%PTR%MAPPINGS%NODES
3125 
3126  ! Initialization
3127  index_ref=1
3128 
3129  DO node_idx=1,nodes_mapping%NUMBER_OF_LOCAL
3130 
3131  !the fourth component of the U1 variable contains the half sarcomere length at activation
3132  dof_idx=field_var_ind_m%COMPONENTS(4)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
3133  & versions(1)
3134  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u1_variable_type, &
3135  & field_values_set_type,dof_idx,sarco_length_at_activation,err,error,*999)
3136 
3137  !the first component of the U1 variable contains the actual half sarcomere length
3138  dof_idx=field_var_ind_m%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node_idx)%DERIVATIVES(1)% &
3139  & versions(1)
3140  CALL field_parameter_set_get_local_dof(independent_field_monodomain,field_u1_variable_type, &
3141  & field_values_set_type,dof_idx,sarco_length,err,error,*999)
3142 
3143  elongation=sarco_length-sarco_length_at_activation
3144 
3145  IF(elongation.LT.0) THEN
3146  length_titin=sarco_length-length_myosin-length_mband-length_zdisc
3147  !function to approximate the relation between the initial titin length and the initial passive force F0
3148  titin_unbound=coeff_matrix(1)*exp(length_titin)+coeff_matrix(2)*length_titin**3+coeff_matrix(3)* &
3149  & length_titin**2+coeff_matrix(4)*length_titin+coeff_matrix(5)
3150  CALL field_parameter_set_update_local_node(independent_field_monodomain, &
3151  & field_u_variable_type,field_values_set_type,1,1,node_idx,2,titin_unbound,err,error,*999)
3152  CALL field_parameter_set_update_local_node(independent_field_monodomain, &
3153  & field_u_variable_type,field_values_set_type,1,1,node_idx,3,titin_unbound,err,error,*999)
3154 
3155  ! calculate x-fibre titin stress (trigonometry):
3156  ! fitting function to calculate the filament-lattice parameter d10 in nanometer (from Elliot 1963 -mammalian muscle)
3157  d10=-13.39_dp*sarco_length+58.37_dp
3158  ! d10=-13.39_DP*1.0_DP+58.37_DP
3159  ! calculate the distance between actin and myosin filament in micro meter (geometrical relationship)
3160  actin_myosin_distance=0.001_dp*(2.0_dp/3.0_dp*d10)
3161  ! calculate x-fibre stress with tangens-function (unbound titin)
3162  titin_xf_unbound=0.5_dp*titin_unbound*actin_myosin_distance/length_titin
3163  CALL field_parameter_set_update_local_node(independent_field_monodomain, &
3164  & field_u_variable_type,field_values_set_type,1,1,node_idx,4,titin_xf_unbound,err,error,*999)
3165  CALL field_parameter_set_update_local_node(independent_field_monodomain, &
3166  & field_u_variable_type,field_values_set_type,1,1,node_idx,5,titin_xf_unbound,err,error,*999)
3167 
3168  ELSE ! Force enhancement
3169 
3170  ! Calculate Titin-Force for unbound titin filaments
3171  length_titin=sarco_length-length_myosin-length_mband-length_zdisc
3172  !function to approximate the relation between the initial titin length and the initial passive force F0
3173  titin_unbound=coeff_matrix(1)*exp(length_titin)+coeff_matrix(2)*length_titin**3+coeff_matrix(3)* &
3174  & length_titin**2+coeff_matrix(4)*length_titin+coeff_matrix(5)
3175 
3176  ! Calculate Titin-Force for bound titin filaments
3177  ! Switch between different force-enhancent models/implementations
3178  ! 1=linear approximation 2=square approximation 3=simple iterative solution of Rode's model
3179  ! 4=Newton's-method to solve Rode's model
3180  switch_model=4
3181  !linear approximation
3182  IF(switch_model.EQ.1) THEN
3183  length_init_titin=sarco_length_at_activation-length_myosin-length_mband-length_zdisc
3184  !function to approximate the relation between the initial titin length and the initial passive force F0
3185  f0=coeff_matrix(1)*exp(length_init_titin)+coeff_matrix(2)*length_init_titin**3+coeff_matrix(3)* &
3186  & length_init_titin**2+coeff_matrix(4)*length_init_titin+coeff_matrix(5)
3187  !function to approximate the relation between the initial sarcomere length and the stiffness of the PEVK region.
3188  stiffness_pevk=1000.0_dp*(0.1880_dp*sarco_length_at_activation**4-0.8694_dp*sarco_length_at_activation**3+ &
3189  & 1.5084_dp*sarco_length_at_activation**2-1.1577_dp*sarco_length_at_activation+0.3345_dp)
3190  IF(elongation.LT.0.02) THEN ! 0.02 -> offset value
3191  force=0.0_dp
3192  ELSE
3193  force=(elongation-0.02_dp)*stiffness_pevk
3194  ENDIF
3195  titin_bound=f0+force
3196 
3197  !square approximation
3198  ELSEIF(switch_model.EQ.2) THEN
3199  length_init_titin=sarco_length_at_activation-length_myosin-length_mband-length_zdisc
3200  !function to approximate the relation between the initial titin length and the initial passive force F0
3201  f0=coeff_matrix(1)*exp(length_init_titin)+coeff_matrix(2)*length_init_titin**3+coeff_matrix(3)* &
3202  & length_init_titin**2+coeff_matrix(4)*length_init_titin+coeff_matrix(5)
3203  force=2.0_dp*elongation**2
3204  titin_bound=f0+force
3205 
3206  !Rode's model -> Ekin's implementation
3207  ELSEIF(switch_model.EQ.3) THEN
3208  length_init_titin=sarco_length_at_activation-length_myosin-length_mband-length_zdisc
3209  !function to approximate the relation between the initial titin length and the initial passive force F0
3210  f0=coeff_matrix(1)*exp(length_init_titin)+coeff_matrix(2)*length_init_titin**3+coeff_matrix(3)* &
3211  & length_init_titin**2+coeff_matrix(4)*length_init_titin+coeff_matrix(5)
3212  !function to approximate the relation between the initial sarcomere length and the stiffness of the PEVK region.
3213  stiffness_pevk=1000.0_dp*(0.1880_dp*sarco_length_at_activation**4-0.8694_dp*sarco_length_at_activation**3+ &
3214  & 1.5084_dp*sarco_length_at_activation**2-1.1577_dp*sarco_length_at_activation+0.3345_dp)
3215 
3216  IF(f0.LE.0) THEN
3217  length_dist_ig_f0=-35.63_dp+39.58889_dp*(f0+0.9_dp)
3218  ELSEIF(f0.GE.0.24_dp) THEN
3219  length_dist_ig_f0=0.1411_dp+0.196576763485477_dp*(f0-0.2495_dp)
3220  ELSE
3221  index_pseudo=ceiling(f0/dx)
3222  index_i=index_ref+index_pseudo-1
3223  length_dist_ig_f0=lengths_dist_ig(index_i)-(lengths_dist_ig(index_i+1)-lengths_dist_ig(index_i))* &
3224  & (forces_dist_ig(index_i)-f0)/(forces_dist_ig(index_i+1)-forces_dist_ig(index_i))
3225  ENDIF
3226 
3227  elongation_new=-1.0_dp
3228  force=0.0_dp
3229  DO WHILE (elongation_new.LT.elongation)
3230  force=force+force_increment
3231  titin_bound=force+f0
3232  elongation_pevk=force/stiffness_pevk
3233  IF (titin_bound.LE.0.0_dp) THEN
3234  length_dist_ig=-35.63_dp+39.58889_dp*(titin_bound+0.9_dp)
3235  ELSE IF (titin_bound.GE.0.24_dp) THEN
3236  length_dist_ig=0.1411_dp+0.196576763485477_dp*(titin_bound-0.2495_dp)
3237  ELSE
3238  index_pseudo=ceiling(titin_bound/dx)
3239  index_i=index_ref+index_pseudo-1
3240  length_dist_ig=lengths_dist_ig(index_i)-(lengths_dist_ig(index_i+1)-lengths_dist_ig( &
3241  & index_i))*(forces_dist_ig(index_i)-titin_bound)/(forces_dist_ig(index_i+1)-forces_dist_ig(index_i))
3242  END IF
3243  elongation_dist_ig=length_dist_ig-length_dist_ig_f0
3244  elongation_new=elongation_pevk+elongation_dist_ig
3245  ENDDO
3246 
3247  !Rode's titin model -> solve with Newton's method
3248  ELSEIF(switch_model.EQ.4) THEN
3249  length_init_titin=sarco_length_at_activation-length_myosin-length_mband-length_zdisc
3250  !function to approximate the relation between the initial titin length and the initial passive force F0
3251  f0=coeff_matrix(1)*exp(length_init_titin)+coeff_matrix(2)*length_init_titin**3+coeff_matrix(3)* &
3252  & length_init_titin**2+coeff_matrix(4)*length_init_titin+coeff_matrix(5)
3253  !function to approximate the relation between the initial sarcomere length and the stiffness of the PEVK region.
3254  stiffness_pevk=1000.0_dp*(0.1880_dp*sarco_length_at_activation**4-0.8694_dp*sarco_length_at_activation**3+ &
3255  & 1.5084_dp*sarco_length_at_activation**2-1.1577_dp*sarco_length_at_activation+0.3345_dp)
3256 
3257  !calculate LENGTH_DIST_IG_F0 with linear inter- or extrapolation
3258  index_i=2
3259  slope=0.0_dp
3260  length_dist_ig_f0=0.0_dp
3261  IF(f0.GE.forces_dist_ig(dim_data)) THEN
3262  index_i=dim_data
3263  ELSE
3264  DO WHILE (f0.GT.forces_dist_ig(index_i))
3265  index_i=index_i+1
3266  ENDDO
3267  ENDIF
3268  slope=(forces_dist_ig(index_i)-forces_dist_ig(index_i-1))/ &
3269  & (lengths_dist_ig(index_i)-lengths_dist_ig(index_i-1))
3270  length_dist_ig_f0=lengths_dist_ig(index_i-1)+slope*(f0-forces_dist_ig(index_i-1))
3271 
3272  !initialize Newton-method to calculate the titin force
3273  index_i=2
3274  stiffness_dist=1.0_dp
3275  force_distal_ig=f0
3276  force=0.0_dp ! delta P
3277  titin_bound=0.0_dp ! total titin force (= P_PEVK)
3278  delta_f=10.0_dp ! start value for iteration
3279  diff_quot=1.0_dp ! numerical derivative
3280  elongation_pevk=elongation/2.0_dp
3281  length_dist_ig=length_dist_ig_f0+elongation-elongation_pevk
3282 
3283  DO WHILE(abs(delta_f).GT.tol) !Newton-method (solve FORCE_DISTAL_IG-FORCE_0=0)
3284 
3285  IF(length_dist_ig.GE.lengths_dist_ig(dim_data)) THEN !Extrapolation if LENGTHS_DIST_IG(end)>LENGTH_DIST_IG
3286  index_i=dim_data
3287  ELSE
3288  DO WHILE(length_dist_ig.GT.lengths_dist_ig(index_i))
3289  index_i=index_i+1
3290  ENDDO
3291  ENDIF
3292  stiffness_dist=(forces_dist_ig(index_i)-forces_dist_ig(index_i-1))/ &
3293  & (lengths_dist_ig(index_i)-lengths_dist_ig(index_i-1))
3294  force_distal_ig=forces_dist_ig(index_i-1)+stiffness_dist* &
3295  & (length_dist_ig-lengths_dist_ig(index_i-1))
3296 
3297  force=stiffness_pevk*elongation_pevk
3298  titin_bound=force+f0
3299 
3300  delta_f=titin_bound-force_distal_ig !new iteration for FORCE_DISTAL_IG-TITIN_BOUND
3301  diff_quot=stiffness_pevk+stiffness_dist
3302  elongation_pevk=elongation_pevk-delta_f/diff_quot
3303  length_dist_ig=length_dist_ig_f0+elongation-elongation_pevk
3304  ENDDO
3305  ENDIF ! switch model
3306 
3307  ! calculate x-fibre titin stress
3308  ! fitting function to calculate the filament-lattice parameter d10 in nanometer (from Elliot 1963 -mammalian muscle)
3309  d10=-13.39_dp*sarco_length+58.37_dp
3310  ! d10=-13.39_DP*1.0_DP+58.37_DP
3311  ! calculate the distance between actin and myosin filament in micro meter (geometrical relationship)
3312  actin_myosin_distance=0.001_dp*(2.0_dp/3.0_dp*d10)
3313  ! calculate x-fibre stress with tangent (unbound titin)
3314  titin_xf_unbound=0.5_dp*titin_unbound*actin_myosin_distance/length_titin
3315  ! calculate x-fibre stress with tangent (for bound titin filaments)
3316  titin_xf_bound=0.5_dp*titin_bound*actin_myosin_distance/length_dist_ig
3317 
3318  CALL field_parameter_set_update_local_node(independent_field_monodomain, &
3319  & field_u_variable_type,field_values_set_type,1,1,node_idx,2,titin_unbound,err,error,*999)
3320  CALL field_parameter_set_update_local_node(independent_field_monodomain, &
3321  & field_u_variable_type,field_values_set_type,1,1,node_idx,3,titin_bound,err,error,*999)
3322  CALL field_parameter_set_update_local_node(independent_field_monodomain, &
3323  & field_u_variable_type,field_values_set_type,1,1,node_idx,4,titin_xf_unbound,err,error,*999)
3324  CALL field_parameter_set_update_local_node(independent_field_monodomain, &
3325  & field_u_variable_type,field_values_set_type,1,1,node_idx,5,titin_xf_bound,err,error,*999)
3326 
3327  ENDIF ! Check if elongation is positive or not
3328  ENDDO ! Over the nodes
3329 
3330  !now the ghost elements -- get the relevant info from the other computational nodes
3331  CALL field_parameter_set_update_start(independent_field_monodomain, &
3332  & field_u_variable_type,field_values_set_type,err,error,*999)
3333  CALL field_parameter_set_update_finish(independent_field_monodomain, &
3334  & field_u_variable_type,field_values_set_type,err,error,*999)
3335 
3336  ELSE
3337  CALL flagerror("Equations set is not associated.",err,error,*999)
3338  ENDIF
3339  ELSE
3340  CALL flagerror("Solver mapping is not associated.",err,error,*999)
3341  ENDIF
3342  ELSE
3343  CALL flagerror("Solver equations is not associated.",err,error,*999)
3344  ENDIF
3345  ENDIF
3346  CASE DEFAULT
3347  CALL flagerror("Problem subtype not implemented for titin",err,error,*999)
3348  END SELECT
3349  ELSE
3350  CALL flagerror("Problem not associated",err,error,*999)
3351  ENDIF
3352  ELSE
3353  CALL flagerror("Control loop is not associated.",err,error,*999)
3354  ENDIF
3355 
3356  exits("BIOELECTRIC_FINITE_ELASTICITY_COMPUTE_TITIN")
3357  RETURN
3358 999 errorsexits("BIOELECTRIC_FINITE_ELASTICITY_COMPUTE_TITIN",err,error)
3359  RETURN 1
3360 
3361  END SUBROUTINE bioelectric_finite_elasticity_compute_titin
3362 
3363  !
3364  !================================================================================================================================
3365  !
3366 
3367 END MODULE bioelectric_finite_elasticity_routines
3368 
equations_set_constants::equations_set_fem_solution_method
integer(intg), parameter equations_set_fem_solution_method
Finite Element Method solution method.
Definition: equations_set_constants.f90:348
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_updategeometricfield
subroutine, public bioelectricfiniteelasticity_updategeometricfield(CONTROL_LOOP, CALC_CLOSEST_GAUSS_POINT, ERR, ERROR,)
Update the the bioelectric equation geometric field from the finite elasticity dependent field (defor...
Definition: bioelectric_finite_elasticity_routines.f90:1658
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
bioelectric_finite_elasticity_routines::bioelectric_finite_elasticity_post_solve
subroutine, public bioelectric_finite_elasticity_post_solve(CONTROL_LOOP, SOLVER, ERR, ERROR,)
Sets up the bioelectrics finite elasticity problem post solve.
Definition: bioelectric_finite_elasticity_routines.f90:609
solver_routines::solvers_create_finish
subroutine, public solvers_create_finish(SOLVERS, ERR, ERROR,)
Finish the creation of solvers.
Definition: solver_routines.f90:22598
control_loop_routines::control_loop_progress_output
integer(intg), parameter, public control_loop_progress_output
Progress output from control loop.
Definition: control_loop_routines.f90:75
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_equationssetsetup
subroutine, public bioelectricfiniteelasticity_equationssetsetup(EQUATIONS_SET, EQUATIONS_SET_SETUP, ERR, ERROR,)
Sets up the bioelectrics finite elasticity equation.
Definition: bioelectric_finite_elasticity_routines.f90:163
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_equationssetsolutionmethodset
subroutine, public bioelectricfiniteelasticity_equationssetsolutionmethodset(EQUATIONS_SET, SOLUTION_METHOD, ERR, ERROR,)
Sets/changes the solution method for a bioelectrics finite elasticity equation type of a multi physic...
Definition: bioelectric_finite_elasticity_routines.f90:99
types::cellml_equations_type
Contains information about the CellML equations for a solver.
Definition: types.f90:2475
types::equations_type
Contains information about the equations in an equations set.
Definition: types.f90:1735
finite_elasticity_routines::finite_elasticity_post_solve
subroutine, public finite_elasticity_post_solve(CONTROL_LOOP, SOLVER, ERR, ERROR,)
Sets up the finite elasticity problem post solve.
Definition: finite_elasticity_routines.f90:9326
equations_set_constants::equations_set_gfem_solution_method
integer(intg), parameter equations_set_gfem_solution_method
Grid-based Finite Element Method solution method.
Definition: equations_set_constants.f90:352
types::region_type
Contains information for a region.
Definition: types.f90:3252
problem_constants::problem_control_time_loop_type
integer(intg), parameter problem_control_time_loop_type
Time control loop.
Definition: problem_constants.f90:283
types::control_loop_time_type
Contains information on a time iteration control loop.
Definition: types.f90:3148
problem_constants::problem_setup_control_type
integer(intg), parameter problem_setup_control_type
Solver setup for a problem.
Definition: problem_constants.f90:263
problem_constants
This module handles all problem wide constants.
Definition: problem_constants.f90:45
problem_constants::solver_equations_first_order_dynamic
integer(intg), parameter solver_equations_first_order_dynamic
Solver equations are first order dynamic.
Definition: problem_constants.f90:311
control_loop_routines::control_loop_node
integer(intg), parameter, public control_loop_node
The identifier for a each "leaf" node in a control loop.
Definition: control_loop_routines.f90:67
maths::matrix_transpose
Returns the transpose of a matrix A in A^T.
Definition: maths.f90:191
solver_routines::solver_dynamic_order_set
subroutine, public solver_dynamic_order_set(SOLVER, ORDER, ERR, ERROR,)
Sets/changes the order for a dynamic solver.
Definition: solver_routines.f90:6327
strings::number_to_vstring
Converts a number to its equivalent varying string representation.
Definition: strings.f90:161
types::decomposition_type
Contains information on the mesh decomposition.
Definition: types.f90:1063
equations_set_constants::equations_set_standard_monodomain_elasticity_subtype
integer(intg), parameter equations_set_standard_monodomain_elasticity_subtype
Definition: equations_set_constants.f90:309
control_loop_routines::control_loop_sub_loop_get
subroutine, public control_loop_sub_loop_get(CONTROL_LOOP, SUB_LOOP_INDEX, SUB_LOOP, ERR, ERROR,)
Gets/returns a pointer to the sub loops as specified by the sub loop index for a control loop...
Definition: control_loop_routines.f90:1417
types::solver_type
Contains information on the type of solver to be used.
Definition: types.f90:2777
solver_routines::solver_petsc_library
integer(intg), parameter, public solver_petsc_library
PETSc solver library.
Definition: solver_routines.f90:104
solver_routines::solvers_number_set
subroutine, public solvers_number_set(SOLVERS, NUMBER_OF_SOLVERS, ERR, ERROR,)
Sets/changes the number of solvers.
Definition: solver_routines.f90:22805
constants::no_part_deriv
integer(intg), parameter no_part_deriv
No partial derivative i.e., u.
Definition: constants.f90:177
solver_routines::solver_dynamic_crank_nicolson_scheme
integer(intg), parameter, public solver_dynamic_crank_nicolson_scheme
Crank-Nicolson dynamic solver.
Definition: solver_routines.f90:282
problem_constants::problem_bioelectric_finite_elasticity_type
integer(intg), parameter problem_bioelectric_finite_elasticity_type
Definition: problem_constants.f90:110
solver_routines::solver_dynamic_degree_set
subroutine, public solver_dynamic_degree_set(SOLVER, DEGREE, ERR, ERROR,)
Sets/changes the degree of the polynomial used to interpolate time for a dynamic solver.
Definition: solver_routines.f90:5484
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_controllooppreloop
subroutine, public bioelectricfiniteelasticity_controllooppreloop(CONTROL_LOOP, ERR, ERROR,)
Sets up the bioelectrics finite elasticity problem pre-control loop.
Definition: bioelectric_finite_elasticity_routines.f90:674
equations_matrices_routines
This module handles all equations matrix and rhs routines.
Definition: equations_matrices_routines.f90:45
solver_routines::solver_dynamic_first_order
integer(intg), parameter, public solver_dynamic_first_order
Dynamic solver has first order terms.
Definition: solver_routines.f90:255
solver_routines::solver_type_set
subroutine, public solver_type_set(SOLVER, SOLVE_TYPE, ERR, ERROR,)
Sets/changes the type for a solver.
Definition: solver_routines.f90:21345
problem_constants::problem_monodomain_elasticity_w_titin_subtype
integer(intg), parameter problem_monodomain_elasticity_w_titin_subtype
Definition: problem_constants.f90:253
types::equations_set_type
Contains information on an equations set.
Definition: types.f90:1941
equations_routines
This module handles all equations routines.
Definition: equations_routines.f90:45
solver_routines::solver_dae_type
integer(intg), parameter, public solver_dae_type
A differential-algebraic equation solver.
Definition: solver_routines.f90:91
strings
This module contains all string manipulation and transformation routines.
Definition: strings.f90:45
solver_routines::solvers_create_start
subroutine, public solvers_create_start(CONTROL_LOOP, SOLVERS, ERR, ERROR,)
Start the creation of a solvers for the control loop.
Definition: solver_routines.f90:22651
types::solvers_type
Contains information on the solvers to be used in a control loop.
Definition: types.f90:2805
bioelectric_finite_elasticity_routines::bioelectric_finite_elasticity_compute_titin
subroutine bioelectric_finite_elasticity_compute_titin(CONTROL_LOOP, ERR, ERROR,)
Computes force enhancement based on the titin model of C Rode et al. (2009). Force depression is not ...
Definition: bioelectric_finite_elasticity_routines.f90:2981
constants::first_part_deriv
integer(intg), parameter first_part_deriv
First partial derivative i.e., du/ds.
Definition: constants.f90:178
problem_constants::solver_equations_static
integer(intg), parameter solver_equations_static
Solver equations are static.
Definition: problem_constants.f90:309
solver_routines::solver_equations_sparsity_type_set
subroutine, public solver_equations_sparsity_type_set(SOLVER_EQUATIONS, SPARSITY_TYPE, ERR, ERROR,)
Sets/changes the sparsity type for solver equations.
Definition: solver_routines.f90:7734
maths
This module contains all mathematics support routines.
Definition: maths.f90:45
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_controllooppostloop
subroutine, public bioelectricfiniteelasticity_controllooppostloop(CONTROL_LOOP, ERR, ERROR,)
Sets up the bioelectrics finite elasticity problem post-control loop.
Definition: bioelectric_finite_elasticity_routines.f90:1347
solver_routines::solvers_solver_get
subroutine, public solvers_solver_get(SOLVERS, SOLVER_INDEX, SOLVER, ERR, ERROR,)
Returns a pointer to the specified solver in the list of solvers.
Definition: solver_routines.f90:22876
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_computefibrestretch
subroutine bioelectricfiniteelasticity_computefibrestretch(CONTROL_LOOP, ERR, ERROR,)
Computes the fibre stretch for bioelectrics finite elasticity problems.
Definition: bioelectric_finite_elasticity_routines.f90:764
types::field_type
Contains information for a field defined on a region.
Definition: types.f90:1346
equations_set_constants::equations_set_1d3d_monodomain_elasticity_subtype
integer(intg), parameter equations_set_1d3d_monodomain_elasticity_subtype
Definition: equations_set_constants.f90:310
problem_constants::solver_equations_linear
integer(intg), parameter solver_equations_linear
Solver equations are linear.
Definition: problem_constants.f90:301
types::control_loop_type
Contains information on a control loop.
Definition: types.f90:3185
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
solver_routines::solver_equations_create_finish
subroutine, public solver_equations_create_finish(SOLVER_EQUATIONS, ERR, ERROR,)
Finishes the process of creating solver equations.
Definition: solver_routines.f90:7166
solver_routines::solver_sparse_matrices
integer(intg), parameter, public solver_sparse_matrices
Use sparse solver matrices.
Definition: solver_routines.f90:360
solver_routines::solver_equations_create_start
subroutine, public solver_equations_create_start(SOLVER, SOLVER_EQUATIONS, ERR, ERROR,)
Starts the process of creating solver equations.
Definition: solver_routines.f90:7209
solver_routines::solver_dynamic_type
integer(intg), parameter, public solver_dynamic_type
A dynamic solver.
Definition: solver_routines.f90:90
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
problem_constants::problem_setup_solvers_type
integer(intg), parameter problem_setup_solvers_type
Solver setup for a problem.
Definition: problem_constants.f90:264
iso_varying_string::varying_string
Definition: iso_varying_string.f90:75
constants
This module contains all program wide constants.
Definition: constants.f90:45
problem_constants::solver_equations_nonlinear
integer(intg), parameter solver_equations_nonlinear
Solver equations are nonlinear.
Definition: problem_constants.f90:302
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_finiteelementcalculate
subroutine, public bioelectricfiniteelasticity_finiteelementcalculate(EQUATIONS_SET, ELEMENT_NUMBER, ERR, ERROR,)
Calculates the element stiffness matrices and RHS for a bioelectrics finite elasticity equation finit...
Definition: bioelectric_finite_elasticity_routines.f90:188
solver_routines::solver_library_type_set
subroutine, public solver_library_type_set(SOLVER, SOLVER_LIBRARY_TYPE, ERR, ERROR,)
Sets/changes the type of library type to use for the solver.
Definition: solver_routines.f90:9067
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_problemspecificationset
subroutine, public bioelectricfiniteelasticity_problemspecificationset(problem, problemSpecification, err, error,)
Sets the problem specification for a bioelectric finite elasticity problem type . ...
Definition: bioelectric_finite_elasticity_routines.f90:213
base_routines::flag_warning
Flags a warning to the user.
Definition: base_routines.f90:232
problem_constants::problem_setup_initial_type
integer(intg), parameter problem_setup_initial_type
Initial setup for a problem.
Definition: problem_constants.f90:262
types::field_interpolated_point_metrics_type
Contains the interpolated point coordinate metrics. Old CMISS name GL,GU,RG.
Definition: types.f90:1112
problem_constants::problem_monodomain_elasticity_velocity_subtype
integer(intg), parameter problem_monodomain_elasticity_velocity_subtype
Definition: problem_constants.f90:254
types::problem_setup_type
Definition: types.f90:3215
solver_routines::solver_equations_linearity_type_set
subroutine, public solver_equations_linearity_type_set(SOLVER_EQUATIONS, LINEARITY_TYPE, ERR, ERROR,)
Sets/changes the linearity type for solver equations.
Definition: solver_routines.f90:7577
solver_routines::cellml_equations_create_start
subroutine, public cellml_equations_create_start(SOLVER, CELLML_EQUATIONS, ERR, ERROR,)
Starts the process of creating CellML equations.
Definition: solver_routines.f90:818
base_routines::exits
subroutine, public exits(NAME)
Records the exit out of the named procedure.
Definition: base_routines.f90:504
control_loop_routines::control_loop_solvers_get
recursive subroutine, public control_loop_solvers_get(CONTROL_LOOP, SOLVERS, ERR, ERROR,)
Returns a pointer to the solvers for a control loop.
Definition: control_loop_routines.f90:1333
types
This module contains all type definitions in order to avoid cyclic module references.
Definition: types.f90:70
solver_routines::solver_cellml_equations_get
subroutine, public solver_cellml_equations_get(SOLVER, CELLML_EQUATIONS, ERR, ERROR,)
Returns a pointer to the CellML equations for a solver.
Definition: solver_routines.f90:947
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
equations_set_constants::equations_set_1d3d_monodomain_active_strain_subtype
integer(intg), parameter equations_set_1d3d_monodomain_active_strain_subtype
Definition: equations_set_constants.f90:313
iso_varying_string::trim
Definition: iso_varying_string.f90:204
control_loop_routines::control_loop_get
Returns the specified control loop as indexed by the control loop identifier from the control loop ro...
Definition: control_loop_routines.f90:86
control_loop_routines::control_loop_type_set
subroutine, public control_loop_type_set(CONTROL_LOOP, LOOP_TYPE, ERR, ERROR,)
Sets/changes the control loop type.
Definition: control_loop_routines.f90:1750
problem_constants::problem_multi_physics_class
integer(intg), parameter problem_multi_physics_class
Definition: problem_constants.f90:63
solver_routines::solver_nonlinear_type
integer(intg), parameter, public solver_nonlinear_type
A nonlinear solver.
Definition: solver_routines.f90:89
solver_routines::solver_solution_update
subroutine, public solver_solution_update(SOLVER, ERR, ERROR,)
Updates the solver solution from the field variables.
Definition: solver_routines.f90:21086
bioelectric_routines::bioelectric_post_solve
subroutine, public bioelectric_post_solve(SOLVER, ERR, ERROR,)
Performs post solve actions for a bioelectrics problem class.
Definition: bioelectric_routines.f90:386
finite_elasticity_routines::finiteelasticity_gaussdeformationgradienttensor
subroutine, public finiteelasticity_gaussdeformationgradienttensor(dependentInterpPointMetrics, geometricInterpPointMetrics, fibreInterpolatedPoint, dZdNu, err, error,)
Evaluates the deformation gradient tensor at a given Gauss point.
Definition: finite_elasticity_routines.f90:3594
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_independentfieldinterpolate
subroutine bioelectricfiniteelasticity_independentfieldinterpolate(CONTROL_LOOP, ERR, ERROR,)
Interpolates the finite elasticity independent field from the biolectrics independent field...
Definition: bioelectric_finite_elasticity_routines.f90:2618
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_convergencecheck
subroutine bioelectricfiniteelasticity_convergencecheck(CONTROL_LOOP, ERR, ERROR,)
Check for the convergence of the bioelectric finite elasticity while loop, i.e., if the force-length ...
Definition: bioelectric_finite_elasticity_routines.f90:1514
problem_constants::problem_setup_finish_action
integer(intg), parameter problem_setup_finish_action
Finish setup action.
Definition: problem_constants.f90:274
equations_mapping_routines
This module handles all equations mapping routines.
Definition: equations_mapping_routines.f90:45
types::solver_equations_type
Contains information about the solver equations for a solver.
Definition: types.f90:2452
types::decomposition_elements_type
Contains the topology information for the elements of a decomposition.
Definition: types.f90:1017
equations_set_constants::equations_set_gfv_solution_method
integer(intg), parameter equations_set_gfv_solution_method
Grid-based Finite Volume solution method.
Definition: equations_set_constants.f90:354
types::problem_type
Contains information for a problem.
Definition: types.f90:3221
problem_constants::problem_setup_cellml_equations_type
integer(intg), parameter problem_setup_cellml_equations_type
CellML equations setup for a problem.
Definition: problem_constants.f90:266
bioelectric_finite_elasticity_routines::bioelectric_finite_elasticity_problem_setup
subroutine, public bioelectric_finite_elasticity_problem_setup(PROBLEM, PROBLEM_SETUP, ERR, ERROR,)
Sets up the bioelectric finite elasticity problem.
Definition: bioelectric_finite_elasticity_routines.f90:269
bioelectric_finite_elasticity_routines::bioelectricfiniteelasticity_forcelengthvelocityrelation
subroutine bioelectricfiniteelasticity_forcelengthvelocityrelation(CONTROL_LOOP, ERR, ERROR,)
Computes the bioelectrics finite elasticity force-length and force_velocity relations.
Definition: bioelectric_finite_elasticity_routines.f90:979
bioelectric_finite_elasticity_routines
This module handles all routines pertaining to bioelectrics coupled with finite elasticity.
Definition: bioelectric_finite_elasticity_routines.f90:47
solver_routines
This module handles all solver routines.
Definition: solver_routines.f90:45
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
solver_routines::solver_dynamic_restart_set
subroutine, public solver_dynamic_restart_set(SOLVER, RESTART, ERR, ERROR,)
Sets/changes the restart value for a dynamic solver.
Definition: solver_routines.f90:6249
field_io_routines
Implements lists of Field IO operation.
Definition: field_IO_routines.f90:46
control_loop_routines::control_loop_create_start
subroutine, public control_loop_create_start(PROBLEM, CONTROL_LOOP, ERR, ERROR,)
Start the process of creating a control loop for a problem.
Definition: control_loop_routines.f90:194
problem_constants::problem_setup_solver_equations_type
integer(intg), parameter problem_setup_solver_equations_type
Solver equations setup for a problem.
Definition: problem_constants.f90:265
equations_set_constants::equations_set_monodomain_elasticity_velocity_subtype
integer(intg), parameter equations_set_monodomain_elasticity_velocity_subtype
Definition: equations_set_constants.f90:312
control_loop_routines::control_loop_label_set
Definition: control_loop_routines.f90:96
biodomain_equation_routines::biodomain_pre_solve
subroutine, public biodomain_pre_solve(SOLVER, ERR, ERROR,)
Performs pre-solve actions for mono- and bi-domain problems.
Definition: biodomain_equation_routines.f90:1611
types::solver_mapping_type
Contains information on the solver mapping between the global equation sets and the solver matrices...
Definition: types.f90:3091
control_loop_routines::control_loop_output_type_set
subroutine, public control_loop_output_type_set(CONTROL_LOOP, OUTPUT_TYPE, ERR, ERROR,)
Sets/changes the output type for a control loop.
Definition: control_loop_routines.f90:1190
solver_routines::solver_dynamic_scheme_set
subroutine, public solver_dynamic_scheme_set(SOLVER, SCHEME, ERR, ERROR,)
Sets/changes the scheme for a dynamic solver.
Definition: solver_routines.f90:6384
finite_elasticity_routines::finite_elasticity_pre_solve
subroutine, public finite_elasticity_pre_solve(CONTROL_LOOP, SOLVER, ERR, ERROR,)
Sets up the finite elasticity problem pre-solve.
Definition: finite_elasticity_routines.f90:9695
types::field_variable_type
Contains information for a field variable defined on a field.
Definition: types.f90:1289
solver_routines::cellml_equations_create_finish
subroutine, public cellml_equations_create_finish(CELLML_EQUATIONS, ERR, ERROR,)
Finishes the process of creating CellML equations.
Definition: solver_routines.f90:780
solver_routines::solver_label_set
Definition: solver_routines.f90:405
equations_set_constants::equations_set_fd_solution_method
integer(intg), parameter equations_set_fd_solution_method
Finite Difference solution method.
Definition: equations_set_constants.f90:350
problem_constants::problem_control_load_increment_loop_type
integer(intg), parameter problem_control_load_increment_loop_type
Load increment control loop.
Definition: problem_constants.f90:285
types::domain_mapping_type
Contains information on the domain mappings (i.e., local and global numberings).
Definition: types.f90:904
types::field_interpolation_parameters_type
Contains the parameters required to interpolate a field variable within an element. Old CMISS name XE.
Definition: types.f90:1141
types::equations_set_setup_type
Contains information on the setup information for an equations set.
Definition: types.f90:1866
problem_constants::problem_setup_start_action
integer(intg), parameter problem_setup_start_action
Start setup action.
Definition: problem_constants.f90:273
solver_routines::solver_equations_time_dependence_type_set
subroutine, public solver_equations_time_dependence_type_set(SOLVER_EQUATIONS, TIME_DEPENDENCE_TYPE, ERR, ERROR,)
Sets/changes the time dependence type for solver equations.
Definition: solver_routines.f90:7788
control_loop_routines
This module handles all control loop routines.
Definition: control_loop_routines.f90:45
solver_routines::solver_cmiss_library
integer(intg), parameter, public solver_cmiss_library
CMISS (internal) solver library.
Definition: solver_routines.f90:103
biodomain_equation_routines
This module handles all bioelectric domain equation routines.
Definition: biodomain_equation_routines.f90:46
maths::matrix_product
Calculates and returns the matrix-product A*B in the matrix C.
Definition: maths.f90:167
base_routines::errors
subroutine, public errors(NAME, ERR, ERROR)
Records the exiting error of the subroutine.
Definition: base_routines.f90:481
equations_set_constants
This module defines all constants shared across equations set routines.
Definition: equations_set_constants.f90:45
bioelectric_routines
This module handles all bioelectric class routines.
Definition: bioelectric_routines.f90:45
equations_set_constants::equations_set_bem_solution_method
integer(intg), parameter equations_set_bem_solution_method
Boundary Element Method solution method.
Definition: equations_set_constants.f90:349
solver_routines::solver_solver_equations_get
subroutine, public solver_solver_equations_get(SOLVER, SOLVER_EQUATIONS, ERR, ERROR,)
Returns a pointer to the solver equations for a solver.
Definition: solver_routines.f90:21230
equations_set_constants::equations_set_monodomain_elasticity_w_titin_subtype
integer(intg), parameter equations_set_monodomain_elasticity_w_titin_subtype
Definition: equations_set_constants.f90:311
types::basis_type
Contains all information about a basis .
Definition: types.f90:184
equations_set_constants::equations_set_fv_solution_method
integer(intg), parameter equations_set_fv_solution_method
Finite Volume solution method.
Definition: equations_set_constants.f90:351
solver_routines::solver_dynamic_first_degree
integer(intg), parameter, public solver_dynamic_first_degree
Dynamic solver uses a first degree polynomial for time interpolation.
Definition: solver_routines.f90:271
control_loop_routines::control_loop_create_finish
recursive subroutine, public control_loop_create_finish(CONTROL_LOOP, ERR, ERROR,)
Finish the process of creating a control loop.
Definition: control_loop_routines.f90:153
problem_constants::problem_monodomain_1d3d_active_strain_subtype
integer(intg), parameter problem_monodomain_1d3d_active_strain_subtype
Definition: problem_constants.f90:255
control_loop_routines::control_loop_number_of_sub_loops_set
subroutine, public control_loop_number_of_sub_loops_set(CONTROL_LOOP, NUMBER_OF_SUB_LOOPS, ERR, ERROR,)
Sets/changes the number of sub loops in a control loop.
Definition: control_loop_routines.f90:1073
bioelectric_finite_elasticity_routines::bioelectric_finite_elasticity_pre_solve
subroutine, public bioelectric_finite_elasticity_pre_solve(CONTROL_LOOP, SOLVER, ERR, ERROR,)
Sets up the bioelectrics finite elasticity problem pre-solve.
Definition: bioelectric_finite_elasticity_routines.f90:536
base_routines::flagerror
Flags an error condition.
Definition: base_routines.f90:240
finite_elasticity_routines::finiteelasticity_controltimelooppreloop
subroutine, public finiteelasticity_controltimelooppreloop(CONTROL_LOOP, ERR, ERROR,)
Runs before each time loop for a finite elasticity problem.
Definition: finite_elasticity_routines.f90:9545
problem_constants::problem_control_while_loop_type
integer(intg), parameter problem_control_while_loop_type
While control loop.
Definition: problem_constants.f90:284
finite_elasticity_routines
This module handles all finite elasticity routines.
Definition: finite_elasticity_routines.f90:45
problem_constants::problem_gudunov_monodomain_1d3d_elasticity_subtype
integer(intg), parameter problem_gudunov_monodomain_1d3d_elasticity_subtype
Definition: problem_constants.f90:252
kinds
This module contains all kind definitions.
Definition: kinds.f90:45
field_io_routines::field_io_nodes_export
subroutine, public field_io_nodes_export(FIELDS, FILE_NAME, METHOD, ERR, ERROR,)
Export nodal information.
Definition: field_IO_routines.f90:6163
biodomain_equation_routines::biodomain_control_loop_post_loop
subroutine, public biodomain_control_loop_post_loop(CONTROL_LOOP, ERR, ERROR,)
Runs after each control loop iteration.
Definition: biodomain_equation_routines.f90:109
problem_constants::problem_gudunov_monodomain_simple_elasticity_subtype
integer(intg), parameter problem_gudunov_monodomain_simple_elasticity_subtype
Definition: problem_constants.f90:251
input_output
This module handles all formating and input and output.
Definition: input_output.f90:45