63 USE generated_mesh_routines
143 INTEGER(INTG),
INTENT(OUT) :: ERR
146 INTEGER(INTG) :: node_idx,component_idx,deriv_idx,variable_idx,dim_idx,local_ny,variable_type
147 INTEGER(INTG) :: NUMBER_OF_DIMENSIONS,user_node,global_node,local_node
148 REAL(DP) :: X(3),DEFORMED_X(3),P
149 REAL(DP),
POINTER :: GEOMETRIC_PARAMETERS(:)
150 TYPE(
domain_type),
POINTER :: DOMAIN,DOMAIN_PRESSURE
156 TYPE(
field_type),
POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD
159 INTEGER(INTG),
ALLOCATABLE :: INNER_SURFACE_NODES(:),OUTER_SURFACE_NODES(:),TOP_SURFACE_NODES(:),BOTTOM_SURFACE_NODES(:)
160 INTEGER(INTG) :: INNER_NORMAL_XI,OUTER_NORMAL_XI,TOP_NORMAL_XI,BOTTOM_NORMAL_XI,MESH_COMPONENT
161 INTEGER(INTG) :: MY_COMPUTATIONAL_NODE_NUMBER, DOMAIN_NUMBER, MPI_IERROR
162 REAL(DP) :: PIN,POUT,LAMBDA,DEFORMED_Z
163 LOGICAL :: X_FIXED,Y_FIXED,NODE_EXISTS, X_OKAY,Y_OKAY
166 NULLIFY(geometric_parameters)
168 enters(
"FiniteElasticity_BoundaryConditionsAnalyticCalculate",err,error,*999)
172 IF(
ASSOCIATED(equations_set))
THEN 173 IF(
ASSOCIATED(equations_set%ANALYTIC))
THEN 174 dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
175 IF(
ASSOCIATED(dependent_field))
THEN 176 geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
177 IF(
ASSOCIATED(geometric_field))
THEN 178 CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
180 NULLIFY(geometric_variable)
181 CALL field_variable_get(geometric_field,field_u_variable_type,geometric_variable,err,error,*999)
182 mesh_component=geometric_variable%COMPONENTS(1)%MESH_COMPONENT_NUMBER
183 CALL field_parameter_set_data_get(geometric_field,field_u_variable_type,field_values_set_type,geometric_parameters, &
186 IF(
ASSOCIATED(boundary_conditions))
THEN 187 decomposition=>dependent_field%DECOMPOSITION
188 IF(
ASSOCIATED(decomposition))
THEN 189 mesh=>decomposition%MESH
190 IF(
ASSOCIATED(mesh))
THEN 191 generated_mesh=>mesh%GENERATED_MESH
192 IF(
ASSOCIATED(generated_mesh))
THEN 193 nodes_mapping=>decomposition%DOMAIN(1)%PTR%MAPPINGS%NODES
194 IF(
ASSOCIATED(nodes_mapping))
THEN 196 CALL generated_mesh_surface_get(generated_mesh,mesh_component,1_intg, &
197 & inner_surface_nodes,inner_normal_xi,err,error,*999)
198 CALL generated_mesh_surface_get(generated_mesh,mesh_component,2_intg, &
199 & outer_surface_nodes,outer_normal_xi,err,error,*999)
200 CALL generated_mesh_surface_get(generated_mesh,mesh_component,3_intg, &
201 & top_surface_nodes,top_normal_xi,err,error,*999)
202 CALL generated_mesh_surface_get(generated_mesh,mesh_component,4_intg, &
203 & bottom_surface_nodes,bottom_normal_xi,err,error,*999)
206 DO node_idx=1,
SIZE(inner_surface_nodes,1)
207 user_node=inner_surface_nodes(node_idx)
209 CALL decomposition_node_domain_get(decomposition,user_node,1,domain_number,err,error,*999)
210 IF(domain_number==my_computational_node_number)
THEN 218 DO node_idx=1,
SIZE(outer_surface_nodes,1)
219 user_node=outer_surface_nodes(node_idx)
221 CALL decomposition_node_domain_get(decomposition,user_node,1,domain_number,err,error,*999)
222 IF(domain_number==my_computational_node_number)
THEN 230 DO node_idx=1,
SIZE(top_surface_nodes,1)
231 user_node=top_surface_nodes(node_idx)
233 CALL decomposition_node_domain_get(decomposition,user_node,1,domain_number,err,error,*999)
234 IF(domain_number==my_computational_node_number)
THEN 235 CALL meshtopologynodecheckexists(mesh,1,user_node,node_exists,global_node,err,error,*999)
236 IF(.NOT.node_exists) cycle
239 local_ny=geometric_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(local_node)% &
240 & derivatives(1)%VERSIONS(1)
241 deformed_z=geometric_parameters(local_ny)*lambda
248 DO node_idx=1,
SIZE(bottom_surface_nodes,1)
249 user_node=bottom_surface_nodes(node_idx)
251 CALL decomposition_node_domain_get(decomposition,user_node,1,domain_number,err,error,*999)
252 IF(domain_number==my_computational_node_number)
THEN 262 DO node_idx=1,
SIZE(bottom_surface_nodes,1)
263 user_node=bottom_surface_nodes(node_idx)
264 CALL decomposition_node_domain_get(decomposition,user_node,1,domain_number,err,error,*999)
265 IF(domain_number==my_computational_node_number)
THEN 266 CALL meshtopologynodecheckexists(mesh,1,user_node,node_exists,global_node,err,error,*999)
267 IF(.NOT.node_exists) cycle
270 local_ny=geometric_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(local_node)% &
271 & derivatives(1)%VERSIONS(1)
272 x(1)=geometric_parameters(local_ny)
273 CALL meshtopologynodecheckexists(mesh,1,user_node,node_exists,global_node,err,error,*999)
274 IF(.NOT.node_exists) cycle
278 local_ny=geometric_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(local_node)% &
279 & derivatives(1)%VERSIONS(1)
280 x(2)=geometric_parameters(local_ny)
281 IF(abs(x(1))<1e-7_dp)
THEN 288 IF(abs(x(2))<1e-7_dp)
THEN 298 CALL mpi_reduce(x_fixed,x_okay,1,mpi_logical,mpi_lor,0,mpi_comm_world,mpi_ierror)
299 CALL mpi_reduce(y_fixed,y_okay,1,mpi_logical,mpi_lor,0,mpi_comm_world,mpi_ierror)
300 IF(my_computational_node_number==0)
THEN 301 IF(.NOT.(x_okay.AND.y_okay))
THEN 302 CALL flagerror(
"Could not fix nodes to prevent rigid body motion",err,error,*999)
306 CALL flagerror(
"Domain nodes mapping is not associated.",err,error,*999)
309 CALL flagerror(
"Generated mesh is not associated. For the Cylinder analytic solution, "// &
310 &
"it must be available for automatic boundary condition assignment",err,error,*999)
313 CALL flagerror(
"Mesh is not associated",err,error,*999)
316 CALL flagerror(
"Decomposition is not associated",err,error,*999)
320 DO variable_idx=1,dependent_field%NUMBER_OF_VARIABLES
321 variable_type=dependent_field%VARIABLES(variable_idx)%VARIABLE_TYPE
322 field_variable=>dependent_field%VARIABLE_TYPE_MAP(variable_type)%PTR
324 & field_variable%NUMBER_OF_GLOBAL_DOFS,err,error,*999)
325 IF(
ASSOCIATED(field_variable))
THEN 326 CALL field_parameter_set_create(dependent_field,variable_type,field_analytic_values_set_type,err,error,*999)
328 IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE==field_node_based_interpolation)
THEN 329 domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
330 IF(
ASSOCIATED(domain))
THEN 331 IF(
ASSOCIATED(domain%TOPOLOGY))
THEN 332 domain_nodes=>domain%TOPOLOGY%NODES
333 IF(
ASSOCIATED(domain_nodes))
THEN 335 IF(field_variable%COMPONENTS(4)%INTERPOLATION_TYPE==field_node_based_interpolation)
THEN 336 domain_pressure=>field_variable%COMPONENTS(4)%DOMAIN
337 IF(
ASSOCIATED(domain_pressure))
THEN 338 IF(
ASSOCIATED(domain_pressure%TOPOLOGY))
THEN 339 domain_pressure_nodes=>domain_pressure%TOPOLOGY%NODES
340 IF(
ASSOCIATED(domain_pressure_nodes))
THEN 343 DO node_idx=1,domain_nodes%NUMBER_OF_NODES
345 DO dim_idx=1,number_of_dimensions
347 local_ny=geometric_variable%COMPONENTS(dim_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% &
348 & nodes(node_idx)%DERIVATIVES(1)%VERSIONS(1)
349 x(dim_idx)=geometric_parameters(local_ny)
352 DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
353 SELECT CASE(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE)
356 SELECT CASE(variable_type)
357 CASE(field_u_variable_type)
358 SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
362 & equations_set%ANALYTIC%ANALYTIC_USER_PARAMS,deformed_x,p,err,error,*999)
364 CALL flagerror(
"Not implemented.",err,error,*999)
366 CALL flagerror(
"Not implemented.",err,error,*999)
368 CALL flagerror(
"Not implemented.",err,error,*999)
371 domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,
"*", &
372 & err,error))//
" is invalid." 373 CALL flagerror(local_error,err,error,*999)
375 CASE(field_deludeln_variable_type)
376 SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
380 CALL flagerror(
"Not implemented.",err,error,*999)
382 CALL flagerror(
"Not implemented.",err,error,*999)
384 CALL flagerror(
"Not implemented.",err,error,*999)
387 domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,
"*", &
388 & err,error))//
" is invalid." 389 CALL flagerror(local_error,err,error,*999)
394 CALL flagerror(local_error,err,error,*999)
397 local_error=
"The analytic function type of "// &
400 CALL flagerror(local_error,err,error,*999)
403 DO component_idx=1,number_of_dimensions
405 local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
406 & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
407 CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
408 & field_analytic_values_set_type,local_ny,deformed_x(component_idx),err,error,*999)
411 user_node=domain_nodes%NODES(node_idx)%USER_NUMBER
412 CALL meshtopologynodecheckexists(mesh,domain_pressure%MESH_COMPONENT_NUMBER,user_node, &
413 & node_exists,global_node,err,error,*999)
415 CALL decomposition_node_domain_get(decomposition,user_node, &
416 & domain_pressure%MESH_COMPONENT_NUMBER,domain_number,err,error,*999)
417 IF(domain_number==my_computational_node_number)
THEN 419 local_node=domain_pressure%mappings%nodes%global_to_local_map(global_node)%local_number(1)
421 local_ny=field_variable%COMPONENTS(4)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% &
422 & nodes(local_node)%DERIVATIVES(deriv_idx)%VERSIONS(1)
425 CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
426 & field_analytic_values_set_type,local_ny,p/2.0_dp,err,error,*999)
433 CALL flagerror(
"Domain for pressure topology node is not associated",err,error,*999)
436 CALL flagerror(
"Domain for pressure topology is not associated",err,error,*999)
439 CALL flagerror(
"Domain for pressure component is not associated",err,error,*999)
442 CALL flagerror(
"Non-nodal based interpolation of pressure cannot be used with analytic solutions", &
446 CALL flagerror(
"Domain topology nodes is not associated.",err,error,*999)
449 CALL flagerror(
"Domain topology is not associated.",err,error,*999)
452 CALL flagerror(
"Domain is not associated.",err,error,*999)
455 CALL flagerror(
"Only node based interpolation is implemented.",err,error,*999)
457 CALL field_parameter_set_update_start(dependent_field,variable_type,field_analytic_values_set_type, &
459 CALL field_parameter_set_update_finish(dependent_field,variable_type,field_analytic_values_set_type, &
462 CALL flagerror(
"Field variable is not associated.",err,error,*999)
466 CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type,field_values_set_type, &
467 & geometric_parameters,err,error,*999)
469 CALL flagerror(
"Boundary conditions is not associated.",err,error,*999)
472 CALL flagerror(
"Equations set geometric field is not associated.",err,error,*999)
475 CALL flagerror(
"Equations set dependent field is not associated.",err,error,*999)
478 CALL flagerror(
"Equations set analytic is not associated.",err,error,*999)
481 CALL flagerror(
"Equations set is not associated.",err,error,*999)
485 exits(
"FiniteElasticity_BoundaryConditionsAnalyticCalculate")
487 999
errors(
"FiniteElasticity_BoundaryConditionsAnalyticCalculate",err,error)
488 exits(
"FiniteElasticity_BoundaryConditionsAnalyticCalculate")
500 REAL(DP),
INTENT(IN) :: X(:)
501 REAL(DP),
INTENT(IN) :: ANALYTIC_USER_PARAMS(:)
502 REAL(DP),
INTENT(OUT) :: DEFORMED_X(3)
503 REAL(DP),
INTENT(OUT) :: P
504 INTEGER(INTG),
INTENT(OUT) :: ERR
507 REAL(DP) :: PIN,POUT,LAMBDA,TSI,A1,A2,C1,C2
508 REAL(DP) :: MU1,MU2,MU,K
511 REAL(DP) :: DEFORMED_R,DEFORMED_THETA
512 REAL(DP) :: DELTA,RES
513 REAL(DP),
PARAMETER :: STEP=1e-5_dp, reltol=1e-12_dp
516 enters(
"FiniteElasticity_CylinderAnalyticCalculate",err,error,*999)
547 ELSEIF (delta<1e-3_dp)
THEN 548 CALL flagerror(
"FiniteElasticity_CylinderAnalyticCalculate failed to converge.",err,error,*999)
556 res=delta/(1.0_dp+mu1)
561 mu2=sqrt(((a1/a2)**2*(lambda*mu1**2-1.0_dp)+1.0_dp)/lambda)
564 r=sqrt(x(1)**2+x(2)**2)
565 theta=atan2(x(2),x(1))
568 k=a1**2*(lambda*mu1**2-1.0_dp)
569 mu=sqrt(1.0_dp/lambda*(1.0_dp+k/r**2))
571 deformed_theta=theta+tsi*lambda*x(3)
572 deformed_x(1)=deformed_r*cos(deformed_theta)
573 deformed_x(2)=deformed_r*sin(deformed_theta)
574 deformed_x(3)=lambda*x(3)
577 p=pout-(c1/lambda+c2*lambda)*(1.0_dp/lambda/mu1**2-r**2/(r**2+k)+log(mu**2/mu1**2))+c1*tsi**2*lambda*(r**2-a1**2) &
578 & -2.0_dp*(c1/lambda**2/mu**2+c2*(1.0_dp/lambda**2+1.0_dp/mu**2+tsi**2*r**2))
580 exits(
"FiniteElasticity_CylinderAnalyticCalculate")
582 999 errorsexits(
"FiniteElasticity_CylinderAnalyticCalculate",err,error)
594 REAL(DP) :: FINITE_ELASTICITY_CYLINDER_ANALYTIC_FUNC_EVALUATE
595 REAL(DP) :: MU1,PIN,POUT,LAMBDA,TSI,A1,A2,C1,C2
599 k=a1**2*(lambda*mu1**2-1.0_dp)
600 mu=sqrt(1.0_dp/lambda*(1.0_dp+k/a2**2))
602 finite_elasticity_cylinder_analytic_func_evaluate= &
603 & 2.0_dp*(c1/lambda**2/mu**2 + c2*(1.0_dp/lambda**2+1.0_dp/mu**2+tsi**2*a2**2))+ &
604 & pout-(c1/lambda+c2*lambda)*(1.0_dp/lambda/mu1**2-a2**2/(a2**2+k)+2*log(mu/mu1))+ &
605 & c1*tsi**2*lambda*(a2**2-a1**2)-2.0_dp*(c1/lambda**2/mu**2+c2*(1.0_dp/lambda**2+ &
606 & 1.0_dp/mu**2+tsi**2*a2**2))+pin
617 & materials_interpolated_point,elasticity_tensor,hydro_elasticity_voigt,stress_tensor,dzdnu, &
622 REAL(DP),
INTENT(OUT) :: ELASTICITY_TENSOR(:,:)
623 REAL(DP),
INTENT(OUT) :: HYDRO_ELASTICITY_VOIGT(:)
624 REAL(DP),
INTENT(OUT) :: STRESS_TENSOR(:)
625 REAL(DP),
INTENT(IN) :: DZDNU(:,:)
626 REAL(DP),
INTENT(IN) :: Jznu
627 INTEGER(INTG),
INTENT(IN) :: ELEMENT_NUMBER,GAUSS_POINT_NUMBER
628 INTEGER(INTG),
INTENT(OUT) :: ERR
631 INTEGER(INTG) :: PRESSURE_COMPONENT,i,j,dof_idx
632 REAL(DP) :: P, I1, I3
633 REAL(DP) :: DZDNUT(3,3),AZL(3,3),AZU(3,3),TEMP(3,3)
634 REAL(DP) :: AZLv(6), AZUv(6)
635 REAL(DP) :: TEMPTERM1,TEMPTERM2,VALUE
636 REAL(DP),
POINTER :: C(:)
637 REAL(DP) :: B(6),E(6),DQ_DE(6),Q
638 REAL(DP) :: I3EE(6,6)
639 REAL(DP) :: ADJCC(6,6)
640 REAL(DP) :: AZUE(6,6)
644 enters(
"FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR",err,error,*999)
646 NULLIFY(field_variable,c)
656 CALL invert(azl,azu,i3,err,error,*999)
671 i3ee = reshape([0.0_dp, 4.0_dp*azlv(3), 4.0_dp*azlv(2), 0.0_dp, 0.0_dp,-4.0_dp*azlv(6), &
672 & 4.0_dp*azlv(3), 0.0_dp, 4.0_dp*azlv(1), 0.0_dp,-4.0_dp*azlv(5), 0.0_dp, &
673 & 4.0_dp*azlv(2), 4.0_dp*azlv(1), 0.0_dp, -2.0_dp*azlv(4), 0.0_dp, 0.0_dp, &
674 & 0.0_dp, 0.0_dp, -4.0_dp*azlv(4), -2.0_dp*azlv(3), 2.0_dp*azlv(6), 2.0_dp*azlv(5), &
675 & 0.0_dp, -4.0_dp*azlv(5), 0.0_dp, 2.0_dp*azlv(6), -2.0_dp*azlv(2), 2.0_dp*azlv(4), &
676 & -4.0_dp*azlv(6), 0.0_dp, 0.0_dp, 2.0_dp*azlv(5), 2.0_dp*azlv(4), -2.0_dp*azlv(1)], [6,6])
677 adjcc = reshape([0.0_dp, azlv(3), azlv(2), 0.0_dp, 0.0_dp,-azlv(6), &
678 & azlv(3), 0.0_dp, azlv(1), 0.0_dp,-azlv(5), 0.0_dp, &
679 & azlv(2), azlv(1), 0.0_dp, -azlv(4), 0.0_dp, 0.0_dp, &
680 & 0.0_dp, 0.0_dp, -azlv(4), -0.5_dp*azlv(3), 0.5_dp*azlv(6), 0.5_dp*azlv(5), &
681 & 0.0_dp, -azlv(5), 0.0_dp,0.5_dp*azlv(6), -0.5_dp*azlv(2), 0.5_dp*azlv(4), &
682 & -azlv(6), 0.0_dp, 0.0_dp, 0.5_dp*azlv(5), 0.5_dp*azlv(4), -0.5_dp*azlv(1)], [6,6])
691 azue(i,j) = -2.0_dp*azuv(i)*azuv(j) + 0.5_dp*i3ee(i,j)/i3
697 elasticity_tensor=0.0_dp
699 SELECT CASE(equations_set%specification(3))
701 local_error=
"Analytic Jacobian has not been validated for the Mooney-Rivlin equations, please use finite differences instead." 702 CALL flagerror(local_error,err,error,*999)
703 pressure_component=dependent_interpolated_point%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
704 p=dependent_interpolated_point%VALUES(pressure_component,
no_part_deriv)
709 i1=azl(1,1)+azl(2,2)+azl(3,3)
710 tempterm1=-2.0_dp*c(2)
711 tempterm2=2.0_dp*(c(1)+
i1*c(2))
712 stress_tensor(1)=tempterm1*azl(1,1)+tempterm2
713 stress_tensor(2)=tempterm1*azl(2,2)+tempterm2
714 stress_tensor(3)=tempterm1*azl(3,3)+tempterm2
715 stress_tensor(4)=tempterm1*azl(2,1)
716 stress_tensor(5)=tempterm1*azl(3,1)
717 stress_tensor(6)=tempterm1*azl(3,2)
723 CALL field_variable_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,field_variable,err,error,*999)
724 DO i=1,field_variable%NUMBER_OF_COMPONENTS
725 dof_idx=field_variable%COMPONENTS(i)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP% &
726 & gauss_points(gauss_point_number,element_number)
727 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
728 & field_values_set_type,dof_idx,
VALUE,err,error,*999)
729 stress_tensor(i)=stress_tensor(i)+
VALUE 735 tempterm1=4.0_dp*c(2)
736 tempterm2=-2.0_dp*c(2)
737 elasticity_tensor(2,1)=tempterm1
738 elasticity_tensor(3,1)=tempterm1
739 elasticity_tensor(1,2)=tempterm1
740 elasticity_tensor(3,2)=tempterm1
741 elasticity_tensor(1,3)=tempterm1
742 elasticity_tensor(2,3)=tempterm1
743 elasticity_tensor(4,4)=tempterm2
744 elasticity_tensor(5,5)=tempterm2
745 elasticity_tensor(6,6)=tempterm2
747 elasticity_tensor=elasticity_tensor + p*azue
750 hydro_elasticity_voigt = azuv
758 stress_tensor(1:3)=stress_tensor(1:3)+p
761 pressure_component=dependent_interpolated_point%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
762 p=dependent_interpolated_point%VALUES(pressure_component,
no_part_deriv)
763 b=[2.0_dp*c(2),2.0_dp*c(3),2.0_dp*c(3),c(4),c(4),c(3)]
764 e=[0.5_dp*(azl(1,1)-1.0_dp),0.5_dp*(azl(2,2)-1.0_dp),0.5_dp*(azl(3,3)-1.0_dp),azl(2,1),azl(3,1),azl(3,2)]
766 tempterm1=0.5_dp*c(1)*exp(0.5_dp*dot_product(e,dq_de))
768 stress_tensor=tempterm1*dq_de + p*azuv
771 CALL field_variable_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,field_variable,err,error,*999)
772 DO i=1,field_variable%NUMBER_OF_COMPONENTS
773 dof_idx=field_variable%COMPONENTS(i)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP% &
774 & gauss_points(gauss_point_number,element_number)
775 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
776 & field_values_set_type,dof_idx,
VALUE,err,error,*999)
777 stress_tensor(i)=stress_tensor(i)+
VALUE 787 elasticity_tensor(i,j)=tempterm1*dq_de(i)*dq_de(j)
790 b=[2.0_dp*c(2),2.0_dp*c(3),2.0_dp*c(3),c(4),c(4),c(3)]
792 elasticity_tensor(i,i)=elasticity_tensor(i,i)+tempterm1*b(i)
797 elasticity_tensor(i,j)=elasticity_tensor(j,i)
802 elasticity_tensor=elasticity_tensor + p*azue
805 hydro_elasticity_voigt = azuv
812 local_error=
"Analytic Jacobian has not been implemented for the third equations set specification of "// &
814 CALL flagerror(local_error,err,error,*999)
817 exits(
"FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR")
819 999 errorsexits(
"FINITE_ELASTICITY_GAUSS_ELASTICITY_TENSOR",err,error)
834 INTEGER(INTG),
INTENT(IN) :: ELEMENT_NUMBER
835 INTEGER(INTG),
INTENT(OUT) :: ERR
838 INTEGER(INTG) :: FIELD_VAR_TYPE,ng,nh,ns,nhs,ni,mh,ms,mhs,oh
839 INTEGER(INTG) :: PRESSURE_COMPONENT
840 INTEGER(INTG) :: SUM_ELEMENT_PARAMETERS,TOTAL_NUMBER_OF_SURFACE_PRESSURE_CONDITIONS
841 INTEGER(INTG) :: NUMBER_OF_DIMENSIONS,NUMBER_OF_XI
842 INTEGER(INTG) :: ELEMENT_BASE_DOF_INDEX(4),component_idx,component_idx2
843 INTEGER(INTG),
PARAMETER :: OFF_DIAG_COMP(3)=[0,1,3],off_diag_dep_var1(3)=[1,1,2],off_diag_dep_var2(3)=[2,3,3]
844 INTEGER(INTG) :: MESH_COMPONENT_NUMBER,NUMBER_OF_ELEMENT_PARAMETERS(4)
845 REAL(DP) :: DZDNU(3,3),CAUCHY_TENSOR(3,3),HYDRO_ELASTICITY_TENSOR(3,3)
846 REAL(DP) :: JGW_SUB_MAT(3,3)
847 REAL(DP) :: TEMPVEC(3)
848 REAL(DP) :: STRESS_TENSOR(6),ELASTICITY_TENSOR(6,6),HYDRO_ELASTICITY_VOIGT(6)
849 REAL(DP) :: DPHIDZ(3,64,3),DJDZ(64,3)
850 REAL(DP) :: JGW_DPHINS_DZ,JGW_DPHIMS_DZ,PHIMS,PHINS,TEMPTERM
851 REAL(DP) :: Jznu,JGW,SUM1,SUM2
857 & MATERIALS_INTERP_POINT,DEPENDENT_INTERP_POINT
859 & DEPENDENT_INTERP_POINT_METRICS
866 TYPE(
field_type),
POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD,FIBRE_FIELD
870 enters(
"FiniteElasticity_FiniteElementJacobianEvaluate",err,error,*999)
872 IF(
ASSOCIATED(equations_set))
THEN 873 equations=>equations_set%EQUATIONS
874 IF(
ASSOCIATED(equations))
THEN 875 equations_matrices=>equations%EQUATIONS_MATRICES
876 nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
877 jacobian_matrix=>nonlinear_matrices%JACOBIANS(1)%PTR
878 IF(jacobian_matrix%UPDATE_JACOBIAN)
THEN 879 dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
880 geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
881 materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
882 fibre_field=>equations%INTERPOLATION%FIBRE_FIELD
884 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
885 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
888 number_of_dimensions=equations_set%REGION%COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS
889 number_of_xi=dependent_basis%NUMBER_OF_XI
891 equations_mapping=>equations%EQUATIONS_MAPPING
892 nonlinear_mapping=>equations_mapping%NONLINEAR_MAPPING
894 field_variable=>nonlinear_mapping%RESIDUAL_VARIABLES(1)%PTR
895 field_var_type=field_variable%VARIABLE_TYPE
897 pressure_component=field_variable%NUMBER_OF_COMPONENTS
899 boundary_conditions=>equations_set%BOUNDARY_CONDITIONS
901 & rhs_variable,boundary_conditions_variable,err,error,*999)
905 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
906 & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
907 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
908 & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
909 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
910 & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
911 IF(
ASSOCIATED(fibre_field))
THEN 912 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
913 & fibre_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
917 geometric_interp_point=>equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR
918 geometric_interp_point_metrics=>equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR
919 IF(
ASSOCIATED(fibre_field))
THEN 920 fibre_interp_point=>equations%INTERPOLATION%FIBRE_INTERP_POINT(field_u_variable_type)%PTR
922 materials_interp_point=>equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR
923 dependent_interp_point=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_type)%PTR
924 dependent_interp_point_metrics=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT_METRICS(field_var_type)%PTR
926 sum_element_parameters=0
928 DO nh=1,field_variable%NUMBER_OF_COMPONENTS
929 mesh_component_number=field_variable%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
930 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR% &
931 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
933 IF(field_variable%COMPONENTS(nh)%INTERPOLATION_TYPE==field_node_based_interpolation)
THEN 934 number_of_element_parameters(nh)=dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
935 ELSEIF(field_variable%COMPONENTS(nh)%INTERPOLATION_TYPE==field_element_based_interpolation)
THEN 936 number_of_element_parameters(nh)=1
938 element_base_dof_index(nh)=sum_element_parameters
939 sum_element_parameters=sum_element_parameters+number_of_element_parameters(nh)
943 DO ng=1,dependent_quadrature_scheme%NUMBER_OF_GAUSS
945 & dependent_interp_point,err,error,*999)
947 & geometric_interp_point,err,error,*999)
949 & geometric_interp_point_metrics,err,error,*999)
951 & dependent_interp_point_metrics,err,error,*999)
953 & materials_interp_point,err,error,*999)
954 IF(
ASSOCIATED(fibre_field))
THEN 956 & fibre_interp_point,err,error,*999)
959 jznu=dependent_interp_point_metrics%JACOBIAN/geometric_interp_point_metrics%JACOBIAN
960 jgw=dependent_interp_point_metrics%JACOBIAN*dependent_quadrature_scheme%GAUSS_WEIGHTS(ng)
963 DO nh=1,number_of_dimensions
964 DO ns=1,number_of_element_parameters(nh)
967 DO mh=1,number_of_dimensions
971 & dependent_interp_point_metrics%DXI_DX(ni,mh)
973 & dependent_interp_point_metrics%DXI_DX(ni,mh)*dependent_interp_point_metrics%GU(ni,mh)
975 dphidz(mh,ns,nh)=sum1
977 djdz(ns,nh)=sum2*dependent_interp_point_metrics%JACOBIAN
982 & geometric_interp_point_metrics,fibre_interp_point,dzdnu,err,error,*999)
985 & materials_interp_point,elasticity_tensor,hydro_elasticity_voigt,stress_tensor, &
986 & dzdnu,jznu,element_number,ng,err,error,*999)
989 DO nh=1,number_of_dimensions
990 DO mh=1,number_of_dimensions
992 hydro_elasticity_tensor(mh,nh)=hydro_elasticity_voigt(
tensor_to_voigt3(mh,nh))
999 DO nh=1,number_of_dimensions
1001 DO ns=1,number_of_element_parameters(nh)
1002 tempvec=matmul(jgw_sub_mat,dphidz(:,ns,nh))
1006 DO ms=ns,number_of_element_parameters(nh)
1008 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)+ &
1009 & dot_product(dphidz(:,ms,nh),tempvec)
1010 DO component_idx=1,number_of_dimensions
1011 DO component_idx2=1,number_of_dimensions
1012 tempterm=cauchy_tensor(component_idx,component_idx2)* &
1013 & dphidz(component_idx2,ms,component_idx)
1024 DO oh=1,off_diag_comp(number_of_dimensions)
1025 nh=off_diag_dep_var1(oh)
1026 mh=off_diag_dep_var2(oh)
1027 nhs=element_base_dof_index(nh)
1029 DO ns=1,number_of_element_parameters(nh)
1031 tempvec=matmul(jgw_sub_mat,dphidz(:,ns,nh))
1033 mhs=element_base_dof_index(mh)
1034 DO ms=1,number_of_element_parameters(mh)
1036 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)+ &
1037 & dot_product(dphidz(:,ms,mh),tempvec)
1038 DO component_idx=1,number_of_dimensions
1039 DO component_idx2=1,number_of_dimensions
1040 tempterm=cauchy_tensor(component_idx,component_idx2)* &
1041 & dphidz(component_idx2,ms,component_idx)
1052 IF(field_variable%COMPONENTS(pressure_component)%INTERPOLATION_TYPE==field_node_based_interpolation)
THEN 1054 DO nh=1,number_of_dimensions
1055 DO ns=1,number_of_element_parameters(nh)
1056 jgw_dphins_dz=jgw*dphidz(nh,ns,nh)
1059 mhs=element_base_dof_index(pressure_component)
1060 DO ms=1,number_of_element_parameters(pressure_component)
1062 phims=quadrature_schemes(pressure_component)%PTR%GAUSS_BASIS_FNS(ms,
no_part_deriv,ng)
1063 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)+ &
1064 & jgw_dphins_dz*phims
1068 ELSEIF(field_variable%COMPONENTS(pressure_component)%INTERPOLATION_TYPE==field_element_based_interpolation)
THEN 1070 DO nh=1,number_of_dimensions
1071 DO ns=1,number_of_element_parameters(nh)
1072 jgw_dphins_dz=jgw*dphidz(nh,ns,nh)
1075 mhs=element_base_dof_index(pressure_component)+1
1076 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)+ &
1084 IF(field_variable%COMPONENTS(pressure_component)%INTERPOLATION_TYPE==field_node_based_interpolation)
THEN 1086 DO mh=1,number_of_dimensions
1087 DO ms=1,number_of_element_parameters(mh)
1088 tempvec=matmul(hydro_elasticity_tensor,dphidz(:,ms,mh))
1089 jgw_dphims_dz=jgw*tempvec(mh)
1092 nhs=element_base_dof_index(pressure_component)
1093 DO ns=1,number_of_element_parameters(pressure_component)
1095 phins=quadrature_schemes(pressure_component)%PTR%GAUSS_BASIS_FNS(ns,
no_part_deriv,ng)
1096 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)+ &
1097 & jgw_dphims_dz*phins
1101 ELSEIF(field_variable%COMPONENTS(pressure_component)%INTERPOLATION_TYPE==field_element_based_interpolation)
THEN 1103 DO mh=1,number_of_dimensions
1104 DO ms=1,number_of_element_parameters(mh)
1105 tempvec=matmul(hydro_elasticity_tensor,dphidz(:,ms,mh))
1106 jgw_dphims_dz=jgw*tempvec(mh)
1109 nhs=element_base_dof_index(pressure_component)+1
1110 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs) + &
1119 IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling)
THEN 1122 CALL field_interpolationparametersscalefactorselementget(element_number, &
1123 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR,err,error,*999)
1126 DO nh=1,number_of_dimensions
1127 DO ns=1,number_of_element_parameters(nh)
1131 DO ms=ns,number_of_element_parameters(nh)
1133 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)* &
1134 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,nh)* &
1135 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
1139 DO oh=1,off_diag_comp(number_of_dimensions)
1140 nh=off_diag_dep_var1(oh)
1141 mh=off_diag_dep_var2(oh)
1142 nhs=element_base_dof_index(nh)
1143 DO ns=1,number_of_element_parameters(nh)
1145 mhs=element_base_dof_index(mh)
1147 DO ms=1,number_of_element_parameters(mh)
1149 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)* &
1150 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
1151 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
1157 IF(field_variable%COMPONENTS(pressure_component)%INTERPOLATION_TYPE==field_node_based_interpolation)
THEN 1159 DO nh=1,number_of_dimensions
1160 DO ns=1,number_of_element_parameters(nh)
1163 mhs=element_base_dof_index(pressure_component)
1164 DO ms=1,number_of_element_parameters(pressure_component)
1166 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(nhs,mhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(nhs,mhs)* &
1167 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR% &
1168 & scale_factors(ms,pressure_component)* &
1169 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
1170 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)* &
1171 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR% &
1172 & scale_factors(ms,pressure_component)* &
1173 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
1177 ELSEIF(field_variable%COMPONENTS(pressure_component)%INTERPOLATION_TYPE==field_element_based_interpolation)
THEN 1179 DO nh=1,number_of_dimensions
1180 DO ns=1,number_of_element_parameters(nh)
1183 mhs=element_base_dof_index(pressure_component)+1
1184 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)* &
1185 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
1192 DO nhs=2,element_base_dof_index(pressure_component)
1194 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(nhs,mhs)
1200 IF(dependent_field%DECOMPOSITION%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)%BOUNDARY_ELEMENT.AND. &
1201 & total_number_of_surface_pressure_conditions>0)
THEN 1206 CALL flagerror(
"Equations set equations is not associated.",err,error,*999)
1209 CALL flagerror(
"Equations set is not associated.",err,error,*999)
1212 exits(
"FiniteElasticity_FiniteElementJacobianEvaluate")
1214 999
errors(
"FiniteElasticity_FiniteElementJacobianEvaluate",err,error)
1215 exits(
"FiniteElasticity_FiniteElementJacobianEvaluate")
1228 REAL(DP),
INTENT(INOUT) :: ELASTICITY_TENSOR(6,6)
1229 REAL(DP),
INTENT(IN) :: DZDNU(3,3)
1230 REAL(DP),
INTENT(IN) :: Jznu
1231 INTEGER(INTG),
INTENT(OUT) :: ERR
1234 INTEGER(INTG) :: i,j
1237 enters(
"FINITE_ELASTICITY_PUSH_ELASTICITY_TENSOR",err,error,*999)
1251 elasticity_tensor=matmul(matmul(t,elasticity_tensor),transpose(t))/jznu
1253 exits(
"FINITE_ELASTICITY_PUSH_ELASTICITY_TENSOR")
1255 999 errorsexits(
"FINITE_ELASTICITY_PUSH_ELASTICITY_TENSOR",err,error)
1267 REAL(DP),
INTENT(INOUT) :: STRESS_TENSOR(6)
1268 REAL(DP),
INTENT(IN) :: DZDNU(3,3)
1269 REAL(DP),
INTENT(IN) :: Jznu
1270 INTEGER(INTG),
INTENT(OUT) :: ERR
1273 INTEGER(INTG) :: i,j
1276 enters(
"FINITE_ELASTICITY_PUSH_STRESS_TENSOR",err,error,*999)
1290 stress_tensor=matmul(t,stress_tensor)/jznu
1292 exits(
"FINITE_ELASTICITY_PUSH_STRESS_TENSOR")
1294 999 errorsexits(
"FINITE_ELASTICITY_PUSH_STRESS_TENSOR",err,error)
1304 & deformationgradienttensor,growthtensor,elasticdeformationgradienttensor,jg,je,err,error,*)
1308 INTEGER(INTG),
INTENT(IN) :: numberOfDimensions
1309 INTEGER(INTG),
INTENT(IN) :: gaussPointNumber
1310 INTEGER(INTG),
INTENT(IN) :: elementNumber
1312 REAL(DP),
INTENT(IN) :: deformationGradientTensor(3,3)
1313 REAL(DP),
INTENT(OUT) :: growthTensor(3,3)
1314 REAL(DP),
INTENT(OUT) :: elasticDeformationGradientTensor(3,3)
1315 REAL(DP),
INTENT(OUT) :: Jg
1316 REAL(DP),
INTENT(OUT) :: Je
1317 INTEGER(INTG),
INTENT(OUT) :: err
1320 REAL(DP) :: growthTensorInverse(3,3),J
1322 enters(
"FiniteElasticity_GaussGrowthTensor",err,error,*999)
1324 IF(
ASSOCIATED(equationsset))
THEN 1327 elasticdeformationgradienttensor=deformationgradienttensor
1328 je=
determinant(elasticdeformationgradienttensor,err,error)
1330 CALL flagerror(
"Equations set is not associated.",err,error,*999)
1337 CALL writestringmatrix(
diagnostic_output_type,1,1,3,1,1,3,3,3,deformationgradienttensor, &
1338 &
write_string_matrix_name_and_indices,
'(" F',
'(",I1,",:)',
' :",3(X,E13.6))',
'(13X,3(X,E13.6))',err,error,*999)
1339 j=
determinant(deformationgradienttensor,err,error)
1342 CALL writestringmatrix(
diagnostic_output_type,1,1,3,1,1,3,3,3,elasticdeformationgradienttensor, &
1343 &
write_string_matrix_name_and_indices,
'(" Fe',
'(",I1,",:)',
' :",3(X,E13.6))',
'(13X,3(X,E13.6))',err,error,*999)
1346 CALL writestringmatrix(
diagnostic_output_type,1,1,3,1,1,3,3,3,growthtensor, &
1347 &
write_string_matrix_name_and_indices,
'(" Fg',
'(",I1,",:)',
' :",3(X,E13.6))',
'(13X,3(X,E13.6))',err,error,*999)
1351 exits(
"FiniteElasticity_GaussGrowthTensor")
1353 999 errorsexits(
"FiniteElasticity_GaussGrowthTensor",err,error)
1364 jacobian,greenstraintensor,err,error,*)
1367 REAL(DP),
INTENT(IN) :: deformationGradientTensor(3,3)
1368 REAL(DP) :: deformationGradientTensorT(3,3)
1369 REAL(DP),
INTENT(OUT) :: rightCauchyDeformationTensor(3,3)
1370 REAL(DP),
INTENT(OUT) :: fingerDeformationTensor(3,3)
1371 REAL(DP),
INTENT(OUT) :: Jacobian
1372 REAL(DP),
INTENT(OUT) :: greenStrainTensor(3,3)
1373 INTEGER(INTG),
INTENT(OUT) :: err
1379 enters(
"FiniteElasticity_StrainTensor",err,error,*999)
1381 CALL matrixtranspose(deformationgradienttensor, deformationgradienttensort,err,error,*999)
1382 CALL matrixproduct(deformationgradienttensort, deformationgradienttensor, rightcauchydeformationtensor,err,error,*999)
1384 CALL invert(rightcauchydeformationtensor,fingerdeformationtensor,i3,err,error,*999)
1385 jacobian=
determinant(deformationgradienttensor,err,error)
1387 greenstraintensor=0.5_dp*rightcauchydeformationtensor
1389 greenstraintensor(i,i)=greenstraintensor(i,i)-0.5_dp
1398 &
' :",3(X,E13.6))',
'(12X,3(X,E13.6))',err,error,*999)
1402 &
' :",3(X,E13.6))',
'(12X,3(X,E13.6))',err,error,*999)
1407 &
' :",3(X,E13.6))',
'(12X,3(X,E13.6))',err,error,*999)
1410 exits(
"FiniteElasticity_StrainTensor")
1412 999 errorsexits(
"FiniteElasticity_StrainTensor",err,error)
1426 INTEGER(INTG),
INTENT(IN) :: ELEMENT_NUMBER
1427 INTEGER(INTG),
INTENT(OUT) :: ERR
1430 TYPE(
basis_type),
POINTER :: DEPENDENT_BASIS,COMPONENT_BASIS
1438 TYPE(
field_type),
POINTER :: DEPENDENT_FIELD,FIBRE_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD,EQUATIONS_SET_FIELD,SOURCE_FIELD
1439 TYPE(
field_type),
POINTER :: INDEPENDENT_FIELD
1443 & FIBRE_INTERPOLATION_PARAMETERS,MATERIALS_INTERPOLATION_PARAMETERS,DEPENDENT_INTERPOLATION_PARAMETERS, &
1444 & DARCY_DEPENDENT_INTERPOLATION_PARAMETERS,SOURCE_INTERPOLATION_PARAMETERS,DARCY_MATERIALS_INTERPOLATION_PARAMETERS, &
1445 & DENSITY_INTERPOLATION_PARAMETERS,INDEPENDENT_INTERPOLATION_PARAMETERS
1447 & MATERIALS_INTERPOLATED_POINT,DEPENDENT_INTERPOLATED_POINT,DARCY_DEPENDENT_INTERPOLATED_POINT,SOURCE_INTERPOLATED_POINT, &
1448 & DENSITY_INTERPOLATED_POINT,INDEPENDENT_INTERPOLATED_POINT,DARCY_MATERIALS_INTERPOLATED_POINT
1450 & DEPENDENT_INTERPOLATED_POINT_METRICS
1451 TYPE(
basis_type),
POINTER :: DEPENDENT_BASIS_1,GEOMETRIC_BASIS
1455 LOGICAL :: DARCY_DENSITY,DARCY_DEPENDENT
1456 INTEGER(INTG) :: component_idx,component_idx2,parameter_idx,gauss_idx,element_dof_idx,FIELD_VAR_TYPE,DARCY_FIELD_VAR_TYPE
1457 INTEGER(INTG) :: imatrix,Ncompartments
1458 INTEGER(INTG) :: i,j,numberOfXDimensions,numberOfXiDimensions
1459 INTEGER(INTG) :: NDOFS,mh,ms,mhs,mi,nh,ns
1460 INTEGER(INTG) :: DEPENDENT_NUMBER_OF_COMPONENTS
1461 INTEGER(INTG) :: NUMBER_OF_DIMENSIONS,NUMBER_OF_XI,HYDROSTATIC_PRESSURE_COMPONENT
1462 INTEGER(INTG) :: NUMBER_OF_FIELD_COMPONENT_INTERPOLATION_PARAMETERS
1463 INTEGER(INTG) :: DEPENDENT_COMPONENT_INTERPOLATION_TYPE
1464 INTEGER(INTG) :: DEPENDENT_NUMBER_OF_GAUSS_POINTS
1465 INTEGER(INTG) :: MESH_COMPONENT_1,MESH_COMPONENT_NUMBER
1466 INTEGER(INTG) :: TOTAL_NUMBER_OF_SURFACE_PRESSURE_CONDITIONS
1467 INTEGER(INTG) :: var1
1468 INTEGER(INTG) :: var2
1469 INTEGER(INTG),
POINTER :: EQUATIONS_SET_FIELD_DATA(:)
1470 REAL(DP) :: DZDNU(3,3),DZDNUT(3,3),AZL(3,3),AZU(3,3),Fe(3,3),FeT(3,3),Fg(3,3),C(3,3),f(3,3),E(3,3),I3,P, &
1471 & piolaTensor(3,3),TEMP(3,3)
1472 REAL(DP) :: cauchyTensor(3,3),JGW_CAUCHY_TENSOR(3,3),kirchoffTensor(3,3),STRESS_TENSOR(6)
1473 REAL(DP) :: deformationGradientTensor(3,3),growthTensor(3,3),growthTensorInverse(3,3),growthTensorInverseTranspose(3,3), &
1474 & fibreGrowth,sheetGrowth,normalGrowth,fibreVector(3),sheetVector(3),normalVector(3)
1475 REAL(DP) :: dNudXi(3,3),dXidNu(3,3)
1476 REAL(DP) :: DFDZ(64,3,3)
1477 REAL(DP) :: DPHIDZ(3,64,3)
1478 REAL(DP) :: GAUSS_WEIGHT,Jznu,Jxxi,Jzxi,Je,Jg,JGW
1479 REAL(DP) :: SUM1,TEMPTERM1
1480 REAL(DP) :: THICKNESS
1481 REAL(DP) :: DARCY_MASS_INCREASE,DARCY_VOL_INCREASE,DARCY_RHO_0_F,DENSITY
1482 REAL(DP) :: Mfact, bfact, p0fact
1483 INTEGER(INTG) :: EQUATIONS_SET_SUBTYPE
1485 enters(
"FiniteElasticity_FiniteElementResidualEvaluate",err,error,*999)
1487 NULLIFY(boundary_conditions,boundary_conditions_variable)
1488 NULLIFY(dependent_basis,component_basis)
1489 NULLIFY(equations,equations_mapping,equations_matrices,nonlinear_matrices,rhs_vector)
1490 NULLIFY(dependent_field,fibre_field,geometric_field,materials_field,source_field,independent_field)
1491 NULLIFY(field_variable)
1492 NULLIFY(dependent_quadrature_scheme,component_quadrature_scheme)
1494 NULLIFY(materials_interpolation_parameters,dependent_interpolation_parameters)
1495 NULLIFY(independent_interpolation_parameters,darcy_materials_interpolation_parameters)
1496 NULLIFY(darcy_dependent_interpolation_parameters,density_interpolation_parameters)
1498 NULLIFY(geometric_interpolated_point_metrics,dependent_interpolated_point_metrics)
1499 NULLIFY(materials_interpolated_point,dependent_interpolated_point,darcy_dependent_interpolated_point)
1500 NULLIFY(density_interpolated_point,independent_interpolated_point)
1501 NULLIFY(dependent_basis_1)
1502 NULLIFY(decomposition)
1503 NULLIFY(equations_set_field_data)
1505 IF(
ASSOCIATED(equations_set))
THEN 1506 IF(.NOT.
ALLOCATED(equations_set%SPECIFICATION))
THEN 1507 CALL flagerror(
"Equations set specification is not allocated.",err,error,*999)
1508 ELSE IF(
SIZE(equations_set%SPECIFICATION,1)/=3)
THEN 1509 CALL flagerror(
"Equations set specification must have three entries for a finite elasticity type equations set.", &
1512 equations_set_subtype = equations_set%SPECIFICATION(3)
1513 equations=>equations_set%EQUATIONS
1514 IF(
ASSOCIATED(equations))
THEN 1517 var1=equations%EQUATIONS_MAPPING%NONLINEAR_MAPPING%RESIDUAL_VARIABLES(1)%PTR%VARIABLE_NUMBER
1518 var2=equations%EQUATIONS_MAPPING%RHS_MAPPING%RHS_VARIABLE%VARIABLE_NUMBER
1523 boundary_conditions=>equations_set%BOUNDARY_CONDITIONS
1525 & rhs_variable,boundary_conditions_variable,err,error,*999)
1529 equations_matrices=>equations%EQUATIONS_MATRICES
1530 nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
1531 rhs_vector=>equations_matrices%RHS_VECTOR
1532 equations_mapping =>equations%EQUATIONS_MAPPING
1534 fibre_field =>equations%INTERPOLATION%FIBRE_FIELD
1535 geometric_field =>equations%INTERPOLATION%GEOMETRIC_FIELD
1536 materials_field =>equations%INTERPOLATION%MATERIALS_FIELD
1537 dependent_field =>equations%INTERPOLATION%DEPENDENT_FIELD
1538 source_field =>equations%INTERPOLATION%SOURCE_FIELD
1539 independent_field=>equations%INTERPOLATION%INDEPENDENT_FIELD
1541 decomposition =>dependent_field%DECOMPOSITION
1542 mesh_component_number = decomposition%MESH_COMPONENT_NUMBER
1544 domain_element_mapping=>decomposition%DOMAIN(1)%PTR%MAPPINGS%ELEMENTS
1546 dependent_basis=>decomposition%DOMAIN(mesh_component_number)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)%BASIS
1548 dependent_number_of_gauss_points=dependent_quadrature_scheme%NUMBER_OF_GAUSS
1549 dependent_number_of_components=dependent_field%VARIABLES(var1)%NUMBER_OF_COMPONENTS
1550 geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
1551 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1553 number_of_dimensions=equations_set%REGION%COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS
1554 number_of_xi=decomposition%DOMAIN(mesh_component_number)%PTR%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)%BASIS%NUMBER_OF_XI
1567 darcy_density=.true.
1569 darcy_density=.false.
1578 darcy_dependent=.true.
1580 darcy_dependent=.false.
1584 field_variable=>equations_set%EQUATIONS%EQUATIONS_MAPPING%NONLINEAR_MAPPING%RESIDUAL_VARIABLES(1)%PTR
1585 field_var_type=field_variable%VARIABLE_TYPE
1586 dependent_interpolation_parameters=>equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR
1587 geometric_interpolation_parameters=>equations%INTERPOLATION%GEOMETRIC_INTERP_PARAMETERS(field_u_variable_type)%PTR
1588 IF(
ASSOCIATED(fibre_field))
THEN 1589 fibre_interpolation_parameters=>equations%INTERPOLATION%FIBRE_INTERP_PARAMETERS(field_u_variable_type)%PTR
1591 IF(
ASSOCIATED(materials_field))
THEN 1592 materials_interpolation_parameters=>equations%INTERPOLATION%MATERIALS_INTERP_PARAMETERS(field_u_variable_type)%PTR
1595 IF(darcy_dependent)
THEN 1596 darcy_dependent_interpolation_parameters=>equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_v_variable_type)%PTR
1598 independent_interpolation_parameters=>equations%INTERPOLATION%INDEPENDENT_INTERP_PARAMETERS(field_u_variable_type)%PTR
1604 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
1605 & geometric_interpolation_parameters,err,error,*999)
1606 IF(
ASSOCIATED(fibre_field))
THEN 1607 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
1608 & fibre_interpolation_parameters,err,error,*999)
1610 IF(
ASSOCIATED(materials_field))
THEN 1611 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
1612 & materials_interpolation_parameters,err,error,*999)
1618 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
1619 & dependent_interpolation_parameters,err,error,*999)
1620 IF(darcy_dependent)
THEN 1621 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
1622 & darcy_dependent_interpolation_parameters,err,error,*999)
1624 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
1625 & independent_interpolation_parameters,err,error,*999)
1633 geometric_interpolated_point=>equations%INTERPOLATION%GEOMETRIC_INTERP_POINT(field_u_variable_type)%PTR
1634 geometric_interpolated_point_metrics=>equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR
1635 IF(
ASSOCIATED(fibre_field))
THEN 1636 fibre_interpolated_point=>equations%INTERPOLATION%FIBRE_INTERP_POINT(field_u_variable_type)%PTR
1638 IF(
ASSOCIATED(materials_field))
THEN 1639 materials_interpolated_point=>equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR
1640 density_interpolated_point=>equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_v_variable_type)%PTR
1642 dependent_interpolated_point=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_type)%PTR
1643 dependent_interpolated_point_metrics=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT_METRICS(field_var_type)%PTR
1644 IF(darcy_dependent)
THEN 1645 darcy_dependent_interpolated_point=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_v_variable_type)%PTR
1647 independent_interpolated_point=>equations%INTERPOLATION%INDEPENDENT_INTERP_POINT(field_u_variable_type)%PTR
1649 IF(
ASSOCIATED(source_field))
THEN 1654 SELECT CASE(equations_set_subtype)
1660 DO gauss_idx=1,dependent_number_of_gauss_points
1663 & dependent_interpolated_point,err,error,*999)
1664 CALL field_interpolated_point_metrics_calculate(dependent_basis%NUMBER_OF_XI,dependent_interpolated_point_metrics, &
1667 & geometric_interpolated_point,err,error,*999)
1668 CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,geometric_interpolated_point_metrics, &
1670 IF(
ASSOCIATED(fibre_field))
THEN 1672 & fibre_interpolated_point,err,error,*999)
1675 & materials_interpolated_point,err,error,*999)
1678 DO nh=1,number_of_dimensions
1679 mesh_component_number=field_variable%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
1680 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR% &
1681 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1683 DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
1685 DO mh=1,number_of_dimensions
1687 DO mi=1,number_of_xi
1688 sum1=sum1+dependent_interpolated_point_metrics%DXI_DX(mi,mh)* &
1691 dphidz(mh,ns,nh)=sum1
1697 & geometric_interpolated_point_metrics,fibre_interpolated_point,dzdnu,err,error,*999)
1699 jznu=dependent_interpolated_point_metrics%JACOBIAN/geometric_interpolated_point_metrics%JACOBIAN
1700 jgw=dependent_interpolated_point_metrics%JACOBIAN*dependent_quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
1704 & materials_interpolated_point,stress_tensor,dzdnu,jznu,element_number,gauss_idx,err,error,*999)
1707 DO nh=1,number_of_dimensions
1708 DO mh=1,number_of_dimensions
1715 DO mh=1,number_of_dimensions
1716 mesh_component_number=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
1717 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR% &
1718 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1719 DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
1721 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)+ &
1722 & dot_product(dphidz(:,ms,mh),jgw_cauchy_tensor(:,mh))
1726 jgw=geometric_interpolated_point_metrics%JACOBIAN*dependent_quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
1729 mesh_component_number=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
1730 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR% &
1731 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1733 tempterm1=jgw*(jznu-1.0_dp)
1734 IF(field_variable%COMPONENTS(mh)%INTERPOLATION_TYPE==field_node_based_interpolation)
THEN 1735 DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
1737 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)+ &
1738 & tempterm1*component_quadrature_scheme%GAUSS_BASIS_FNS(ms,
no_part_deriv,gauss_idx)
1740 ELSEIF(field_variable%COMPONENTS(mh)%INTERPOLATION_TYPE==field_element_based_interpolation)
THEN 1742 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)+tempterm1
1775 IF(decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)%BOUNDARY_ELEMENT.AND. &
1776 & total_number_of_surface_pressure_conditions>0)
THEN 1781 IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling)
THEN 1784 CALL field_interpolationparametersscalefactorselementget(element_number, &
1785 & dependent_interpolation_parameters,err,error,*999)
1787 DO mh=1,number_of_dimensions
1788 mesh_component_number=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
1789 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR% &
1790 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1792 DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
1794 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)* &
1795 & dependent_interpolation_parameters%SCALE_FACTORS(ms,mh)
1798 IF(field_variable%COMPONENTS(mh)%INTERPOLATION_TYPE==field_node_based_interpolation)
THEN 1799 mesh_component_number=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
1800 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR% &
1801 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1802 DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
1804 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)* &
1805 & dependent_interpolation_parameters%SCALE_FACTORS(ms,mh)
1829 DO gauss_idx=1,dependent_number_of_gauss_points
1830 gauss_weight=dependent_quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
1833 & dependent_interpolated_point,err,error,*999)
1834 CALL field_interpolated_point_metrics_calculate(dependent_basis%NUMBER_OF_XI,dependent_interpolated_point_metrics, &
1837 & geometric_interpolated_point,err,error,*999)
1838 CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,geometric_interpolated_point_metrics, &
1840 IF(
ASSOCIATED(fibre_field))
THEN 1842 & fibre_interpolated_point,err,error,*999)
1845 & materials_interpolated_point,err,error,*999)
1846 IF(darcy_dependent)
THEN 1848 & darcy_dependent_interpolated_point,err,error,*999)
1851 & independent_interpolated_point,err,error,*999)
1856 & geometric_interpolated_point_metrics,fibre_interpolated_point,dzdnu,err,error,*999)
1858 IF(jznu<0.0_dp)
THEN 1859 local_error =
"Warning: Volume is negative for gauss point "//
trim(
number_to_vstring(gauss_idx,
"*",err,error))//&
1866 jzxi=dependent_interpolated_point_metrics%JACOBIAN
1867 jxxi=geometric_interpolated_point_metrics%JACOBIAN
1876 & dzdnu,fg,fe,jg,je,err,error,*999)
1883 & materials_interpolated_point,darcy_dependent_interpolated_point, &
1884 & independent_interpolated_point,cauchytensor,jznu,dzdnu,element_number,gauss_idx,err,error,*999)
1892 & 3,3,piolatensor,
write_string_matrix_name_and_indices,
'(" T',
'(",I1,",:)',
' :",3(X,E13.6))', &
1893 &
'(12X,3(X,E13.6))',err,error,*999)
1896 & 3,3,cauchytensor,
write_string_matrix_name_and_indices,
'(" sigma',
'(",I1,",:)',
' :",3(X,E13.6))', &
1897 &
'(12X,3(X,E13.6))',err,error,*999)
1903 darcy_mass_increase = darcy_dependent_interpolated_point%VALUES(4,
no_part_deriv)
1904 darcy_vol_increase = darcy_mass_increase / darcy_rho_0_f
1910 IF(number_of_dimensions == 3)
THEN 1911 thickness = materials_interpolated_point%VALUES(materials_interpolated_point%INTERPOLATION_PARAMETERS% &
1912 & field_variable%NUMBER_OF_COMPONENTS,1)
1917 jgw=jzxi*dependent_quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
1920 DO nh=1,number_of_dimensions
1921 mesh_component_number=field_variable%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
1922 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%ptr% &
1923 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1925 DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
1927 DO mh=1,number_of_dimensions
1929 DO mi=1,number_of_xi
1930 sum1=sum1+dependent_interpolated_point_metrics%DXI_DX(mi,mh)* &
1933 dphidz(mh,ns,nh)=sum1
1940 DO mh=1,number_of_dimensions
1941 mesh_component_number=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
1942 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%ptr% &
1943 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1944 DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
1946 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)+ &
1947 & jgw*dot_product(dphidz(1:number_of_dimensions,ms,mh),cauchytensor(1:number_of_dimensions,mh))
1953 hydrostatic_pressure_component=dependent_field%VARIABLES(var1)%NUMBER_OF_COMPONENTS
1954 dependent_component_interpolation_type=dependent_field%VARIABLES(var1)%COMPONENTS(hydrostatic_pressure_component)% &
1955 & interpolation_type
1957 tempterm1=gauss_weight*(jzxi-(jg-darcy_vol_increase)*jxxi)
1959 tempterm1=gauss_weight*(jzxi/jxxi - 1.0_dp)*jxxi
1961 IF(dependent_component_interpolation_type==field_node_based_interpolation)
THEN 1962 component_basis=>dependent_field%VARIABLES(var1)%COMPONENTS(hydrostatic_pressure_component)%DOMAIN% &
1963 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
1965 number_of_field_component_interpolation_parameters=component_basis%NUMBER_OF_ELEMENT_PARAMETERS
1966 DO parameter_idx=1,number_of_field_component_interpolation_parameters
1968 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)+ &
1969 & component_quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,1,gauss_idx)*tempterm1
1971 ELSEIF(dependent_component_interpolation_type==field_element_based_interpolation)
THEN 1973 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)+ tempterm1
1979 IF(decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)%BOUNDARY_ELEMENT.AND. &
1980 & total_number_of_surface_pressure_conditions>0)
THEN 1989 DO gauss_idx=1,dependent_number_of_gauss_points
1996 gauss_weight=dependent_quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
1999 & dependent_interpolated_point,err,error,*999)
2000 CALL field_interpolatedpointmetricscalculate(dependent_basis%NUMBER_OF_XI,dependent_interpolated_point_metrics, &
2003 & geometric_interpolated_point,err,error,*999)
2004 CALL field_interpolatedpointmetricscalculate(geometric_basis%NUMBER_OF_XI,geometric_interpolated_point_metrics, &
2006 IF(
ASSOCIATED(fibre_field))
THEN 2008 & fibre_interpolated_point,err,error,*999)
2013 & geometric_interpolated_point_metrics,fibre_interpolated_point,dzdnu,err,error,*999)
2015 jxxi=geometric_interpolated_point_metrics%JACOBIAN
2017 jzxi=dependent_interpolated_point_metrics%JACOBIAN
2019 hydrostatic_pressure_component=dependent_interpolated_point%INTERPOLATION_PARAMETERS%FIELD_VARIABLE% &
2020 & number_of_components
2021 p=dependent_interpolated_point%VALUES(hydrostatic_pressure_component,1)
2024 & dzdnu,fg,fe,jg,je,err,error,*999)
2029 IF(number_of_dimensions==3)
THEN 2030 CALL field_parametersetgetlocalgausspoint(dependent_field,field_u2_variable_type,field_values_set_type, &
2031 & gauss_idx,element_number,1,piolatensor(1,1),err,error,*999)
2032 CALL field_parametersetgetlocalgausspoint(dependent_field,field_u2_variable_type,field_values_set_type, &
2033 & gauss_idx,element_number,2,piolatensor(1,2),err,error,*999)
2034 CALL field_parametersetgetlocalgausspoint(dependent_field,field_u2_variable_type,field_values_set_type, &
2035 & gauss_idx,element_number,3,piolatensor(1,3),err,error,*999)
2036 CALL field_parametersetgetlocalgausspoint(dependent_field,field_u2_variable_type,field_values_set_type, &
2037 & gauss_idx,element_number,4,piolatensor(2,2),err,error,*999)
2038 CALL field_parametersetgetlocalgausspoint(dependent_field,field_u2_variable_type,field_values_set_type, &
2039 & gauss_idx,element_number,5,piolatensor(2,3),err,error,*999)
2040 CALL field_parametersetgetlocalgausspoint(dependent_field,field_u2_variable_type,field_values_set_type, &
2041 & gauss_idx,element_number,6,piolatensor(3,3),err,error,*999)
2043 piolatensor(1,1)=piolatensor(1,1)+p*f(1,1)
2044 piolatensor(2,2)=piolatensor(2,2)+p*f(2,2)
2045 piolatensor(3,3)=piolatensor(3,3)+p*f(3,3)
2046 piolatensor(1,2)=piolatensor(1,2)+p*f(1,2)
2047 piolatensor(1,3)=piolatensor(1,3)+p*f(1,3)
2048 piolatensor(2,3)=piolatensor(2,3)+p*f(2,3)
2049 piolatensor(2,1)=piolatensor(1,2)
2050 piolatensor(3,1)=piolatensor(1,3)
2051 piolatensor(3,2)=piolatensor(2,3)
2052 ELSE IF(number_of_dimensions==2)
THEN 2053 CALL field_parametersetgetlocalgausspoint(dependent_field,field_u2_variable_type,field_values_set_type, &
2054 & gauss_idx,element_number,1,piolatensor(1,1),err,error,*999)
2055 CALL field_parametersetgetlocalgausspoint(dependent_field,field_u2_variable_type,field_values_set_type, &
2056 & gauss_idx,element_number,2,piolatensor(1,2),err,error,*999)
2057 CALL field_parametersetgetlocalgausspoint(dependent_field,field_u2_variable_type,field_values_set_type, &
2058 & gauss_idx,element_number,3,piolatensor(2,2),err,error,*999)
2060 piolatensor(1,1)=piolatensor(1,1)+p*f(1,1)
2061 piolatensor(2,2)=piolatensor(2,2)+p*f(2,2)
2062 piolatensor(1,2)=piolatensor(1,2)+p*f(1,2)
2063 piolatensor(2,1)=piolatensor(1,2)
2065 CALL field_parametersetgetlocalgausspoint(dependent_field,field_u2_variable_type,field_values_set_type, &
2066 & gauss_idx,element_number,1,piolatensor(1,1),err,error,*999)
2067 piolatensor(1,1)=piolatensor(1,1)+p*f(1,1)
2075 cauchytensor=kirchofftensor/je
2083 & 3,3,piolatensor,
write_string_matrix_name_and_indices,
'(" T',
'(",I1,",:)',
' :",3(X,E13.6))', &
2084 &
'(12X,3(X,E13.6))',err,error,*999)
2087 & 3,3,cauchytensor,
write_string_matrix_name_and_indices,
'(" sigma',
'(",I1,",:)',
' :",3(X,E13.6))', &
2088 &
'(12X,3(X,E13.6))',err,error,*999)
2099 IF(number_of_dimensions == 3)
THEN 2100 IF(
ASSOCIATED(materials_field))
THEN 2102 & materials_interpolated_point,err,error,*999)
2103 thickness = materials_interpolated_point%VALUES(materials_interpolated_point%INTERPOLATION_PARAMETERS% &
2104 & field_variable%NUMBER_OF_COMPONENTS,1)
2133 DO nh=1,number_of_dimensions
2134 mesh_component_number=field_variable%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
2135 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR% &
2136 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2138 DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
2140 DO mh=1,number_of_dimensions
2142 DO mi=1,number_of_xi
2143 sum1=sum1+dependent_interpolated_point_metrics%DXI_DX(mi,mh)* &
2146 dphidz(mh,ns,nh)=sum1
2150 jgw=jzxi*dependent_quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
2153 DO mh=1,number_of_dimensions
2154 mesh_component_number=field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
2155 dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR% &
2156 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2157 DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
2159 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)+ &
2160 & jgw*dot_product(dphidz(:,ms,mh),cauchytensor(:,mh))
2166 hydrostatic_pressure_component=dependent_field%VARIABLES(var1)%NUMBER_OF_COMPONENTS
2167 dependent_component_interpolation_type=dependent_field%VARIABLES(var1)%COMPONENTS(hydrostatic_pressure_component)% &
2168 & interpolation_type
2170 tempterm1=gauss_weight*(jzxi-(jg-darcy_vol_increase)*jxxi)
2172 tempterm1=gauss_weight*(jzxi-jg*jxxi)
2174 IF(dependent_component_interpolation_type==field_node_based_interpolation)
THEN 2175 component_basis=>dependent_field%VARIABLES(var1)%COMPONENTS(hydrostatic_pressure_component)%DOMAIN% &
2176 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2178 number_of_field_component_interpolation_parameters=component_basis%NUMBER_OF_ELEMENT_PARAMETERS
2179 DO parameter_idx=1,number_of_field_component_interpolation_parameters
2180 element_dof_idx=element_dof_idx+1
2182 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)= &
2183 & nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)+ &
2184 & gauss_weight*jzxi*component_quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,1,gauss_idx)* &
2185 & (je-1.0_dp-darcy_vol_increase)
2187 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)= &
2188 & nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)+ &
2189 & gauss_weight*jzxi*component_quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,1,gauss_idx)* &
2193 ELSEIF(dependent_component_interpolation_type==field_element_based_interpolation)
THEN 2195 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)= &
2196 & nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)+tempterm1
2202 IF(decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)%BOUNDARY_ELEMENT.AND. &
2203 & total_number_of_surface_pressure_conditions>0)
THEN 2212 equations_set_field=>equations_set%EQUATIONS_SET_FIELD%EQUATIONS_SET_FIELD_FIELD
2213 CALL field_parameter_set_data_get(equations_set_field,field_u_variable_type, &
2214 & field_values_set_type,equations_set_field_data,err,error,*999)
2216 ncompartments = equations_set_field_data(2)
2218 DO gauss_idx=1,dependent_number_of_gauss_points
2219 gauss_weight=dependent_quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
2222 & dependent_interpolated_point,err,error,*999)
2223 CALL field_interpolated_point_metrics_calculate(dependent_basis%NUMBER_OF_XI,dependent_interpolated_point_metrics, &
2226 & geometric_interpolated_point,err,error,*999)
2227 CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,geometric_interpolated_point_metrics, &
2229 IF(
ASSOCIATED(fibre_field))
THEN 2231 & fibre_interpolated_point,err,error,*999)
2234 & materials_interpolated_point,err,error,*999)
2242 & 3,3,piolatensor,
write_string_matrix_name_and_indices,
'(" Piola Tensor',
'(",I1,",:)',
' :",3(X,E13.6))', &
2243 &
'(17X,3(X,E13.6))',err,error,*999)
2246 & 3,3,cauchytensor,
write_string_matrix_name_and_indices,
'(" Cauchy Tensor',
'(",I1,",:)',
' :",3(X,E13.6))', &
2247 &
'(17X,3(X,E13.6))',err,error,*999)
2253 darcy_mass_increase = 0.0_dp
2254 DO imatrix=1,ncompartments
2255 darcy_field_var_type=field_v_variable_type+field_number_of_variable_subtypes*(imatrix-1)
2256 darcy_dependent_interpolation_parameters=>&
2257 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(darcy_field_var_type)%PTR
2259 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
2260 & darcy_dependent_interpolation_parameters,err,error,*999)
2262 darcy_dependent_interpolated_point=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT(darcy_field_var_type)%PTR
2264 & darcy_dependent_interpolated_point,err,error,*999)
2266 darcy_mass_increase = darcy_mass_increase + darcy_dependent_interpolated_point%VALUES(4,
no_part_deriv)
2269 darcy_vol_increase = darcy_mass_increase / darcy_rho_0_f
2273 & geometric_interpolated_point_metrics,fibre_interpolated_point,dzdnu,err,error,*999)
2275 jxxi=geometric_interpolated_point_metrics%JACOBIAN
2279 & materials_interpolated_point,darcy_dependent_interpolated_point, &
2280 & independent_interpolated_point,cauchytensor,jznu,dzdnu,element_number,gauss_idx,err,error,*999)
2284 & number_of_xi,dfdz,err,error,*999)
2289 IF(number_of_dimensions == 3)
THEN 2290 thickness = materials_interpolated_point%VALUES(materials_interpolated_point%INTERPOLATION_PARAMETERS% &
2291 & field_variable%NUMBER_OF_COMPONENTS,1)
2297 DO component_idx=1,number_of_dimensions
2298 dependent_component_interpolation_type=dependent_field%VARIABLES(var1)%COMPONENTS(component_idx)%INTERPOLATION_TYPE
2299 IF(dependent_component_interpolation_type==field_node_based_interpolation)
THEN 2300 dependent_basis=>dependent_field%VARIABLES(var1)%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY% &
2301 & elements%ELEMENTS(element_number)%BASIS
2302 number_of_field_component_interpolation_parameters=dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
2303 DO parameter_idx=1,number_of_field_component_interpolation_parameters
2304 element_dof_idx=element_dof_idx+1
2305 DO component_idx2=1,number_of_dimensions
2306 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)= &
2307 & nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)+ &
2308 & gauss_weight*jxxi*jznu*thickness*cauchytensor(component_idx,component_idx2)* &
2309 & dfdz(parameter_idx,component_idx2,component_idx)
2312 ELSEIF(dependent_component_interpolation_type==field_element_based_interpolation)
THEN 2314 CALL flagerror(
"Finite elasticity with element based interpolation is not implemented.",err,error,*999)
2320 hydrostatic_pressure_component=dependent_field%VARIABLES(var1)%NUMBER_OF_COMPONENTS
2321 dependent_component_interpolation_type=dependent_field%VARIABLES(var1)%COMPONENTS(component_idx)%INTERPOLATION_TYPE
2322 IF(dependent_component_interpolation_type==field_node_based_interpolation)
THEN 2323 component_basis=>dependent_field%VARIABLES(var1)%COMPONENTS(hydrostatic_pressure_component)%DOMAIN% &
2324 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2326 number_of_field_component_interpolation_parameters=component_basis%NUMBER_OF_ELEMENT_PARAMETERS
2327 DO parameter_idx=1,number_of_field_component_interpolation_parameters
2328 element_dof_idx=element_dof_idx+1
2329 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)= &
2330 & nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)+ &
2331 & gauss_weight*jxxi*component_quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,1,gauss_idx)* &
2332 & (jznu-1.0_dp-darcy_vol_increase)
2334 ELSEIF(dependent_component_interpolation_type==field_element_based_interpolation)
THEN 2335 element_dof_idx=element_dof_idx+1
2336 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)= &
2337 & nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)+gauss_weight*jxxi* &
2338 & (jznu-1.0_dp-darcy_vol_increase)
2344 IF(decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)%BOUNDARY_ELEMENT.AND. &
2345 & total_number_of_surface_pressure_conditions>0)
THEN 2359 DO gauss_idx=1,dependent_number_of_gauss_points
2360 gauss_weight=dependent_quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
2364 & dependent_interpolated_point,err,error,*999)
2365 CALL field_interpolated_point_metrics_calculate(dependent_basis%NUMBER_OF_XI,dependent_interpolated_point_metrics, &
2368 & geometric_interpolated_point,err,error,*999)
2369 CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,geometric_interpolated_point_metrics, &
2371 IF(
ASSOCIATED(fibre_field))
THEN 2373 & fibre_interpolated_point,err,error,*999)
2376 & materials_interpolated_point,err,error,*999)
2377 IF(darcy_dependent)
THEN 2379 & darcy_dependent_interpolated_point,err,error,*999)
2384 & geometric_interpolated_point_metrics,fibre_interpolated_point,dzdnu,err,error,*999)
2386 jxxi=geometric_interpolated_point_metrics%JACOBIAN
2390 & materials_interpolated_point,darcy_dependent_interpolated_point, &
2391 & independent_interpolated_point,cauchytensor,jznu,dzdnu,element_number,gauss_idx,err,error,*999)
2395 & number_of_xi,dfdz,err,error,*999)
2399 DO component_idx=1,dependent_number_of_components
2400 dependent_component_interpolation_type=dependent_field%VARIABLES(var1)%COMPONENTS(component_idx)%INTERPOLATION_TYPE
2401 IF(dependent_component_interpolation_type==field_node_based_interpolation)
THEN 2402 dependent_basis=>dependent_field%VARIABLES(var1)%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY% &
2403 & elements%ELEMENTS(element_number)%BASIS
2404 number_of_field_component_interpolation_parameters=dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
2405 DO parameter_idx=1,number_of_field_component_interpolation_parameters
2406 element_dof_idx=element_dof_idx+1
2407 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)= &
2408 & nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)+ &
2409 & gauss_weight*jxxi*jznu*(cauchytensor(component_idx,1)*dfdz(parameter_idx,1,component_idx)+ &
2410 & cauchytensor(component_idx,2)*dfdz(parameter_idx,2,component_idx)+ &
2411 & cauchytensor(component_idx,3)*dfdz(parameter_idx,3,component_idx))
2413 ELSEIF(dependent_component_interpolation_type==field_element_based_interpolation)
THEN 2415 CALL flagerror(
"Finite elasticity with element based interpolation is not implemented.",err,error,*999)
2421 IF(decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)%BOUNDARY_ELEMENT.AND. &
2422 & total_number_of_surface_pressure_conditions>0)
THEN 2426 IF(
ASSOCIATED(rhs_vector))
THEN 2427 IF(
ASSOCIATED(source_field))
THEN 2428 IF(
ASSOCIATED(materials_field%VARIABLE_TYPE_MAP(field_v_variable_type)%PTR))
THEN 2429 density_interpolation_parameters=>equations%INTERPOLATION%MATERIALS_INTERP_PARAMETERS(field_v_variable_type)%PTR
2430 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
2431 & density_interpolation_parameters,err,error,*999)
2432 density_interpolated_point=>equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_v_variable_type)%PTR
2433 IF(darcy_density)
THEN 2434 darcy_materials_interpolation_parameters=>equations%INTERPOLATION%MATERIALS_INTERP_PARAMETERS( &
2435 & field_u1_variable_type)%PTR
2436 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
2437 & darcy_materials_interpolation_parameters,err,error,*999)
2438 darcy_materials_interpolated_point=>equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u1_variable_type)%PTR
2440 IF(rhs_vector%UPDATE_VECTOR)
THEN 2442 CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
2446 DO gauss_idx=1,dependent_number_of_gauss_points
2447 gauss_weight=dependent_quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)
2451 & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
2453 & density_interpolated_point,err,error,*999)
2454 IF(darcy_density)
THEN 2456 & darcy_materials_interpolated_point,err,error,*999)
2463 density=density_interpolated_point%VALUES(1,1)*(1.0_dp-darcy_materials_interpolated_point%VALUES(8,1)) + &
2464 & darcy_materials_interpolated_point%VALUES(7,1)*(jznu-1.0_dp+darcy_materials_interpolated_point%VALUES(8,1))
2466 density=density_interpolated_point%VALUES(1,1)
2468 CALL field_interpolated_point_metrics_calculate(dependent_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
2469 & dependent_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)
2471 DO component_idx=1,number_of_dimensions
2472 dependent_basis=>dependent_field%VARIABLES(var1)%COMPONENTS(component_idx)%DOMAIN%TOPOLOGY% &
2473 & elements%ELEMENTS(element_number)%BASIS
2474 DO parameter_idx=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
2475 element_dof_idx=element_dof_idx+1
2476 rhs_vector%ELEMENT_VECTOR%VECTOR(element_dof_idx)=rhs_vector%ELEMENT_VECTOR%VECTOR(element_dof_idx) + &
2478 & dependent_quadrature_scheme%GAUSS_BASIS_FNS(parameter_idx,
no_part_deriv,gauss_idx)*gauss_weight * &
2479 & equations%INTERPOLATION%DEPENDENT_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN
2487 CALL flagerror(
"RHS vector is not associated.",err,error,*999)
2491 IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling)
THEN 2492 CALL field_interpolationparametersscalefactorselementget(element_number,equations%INTERPOLATION% &
2493 & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
2495 DO mh=1,field_variable%NUMBER_OF_COMPONENTS
2497 dependent_component_interpolation_type=dependent_field%VARIABLES(field_var_type)%COMPONENTS(mh)%INTERPOLATION_TYPE
2498 IF(dependent_component_interpolation_type==field_node_based_interpolation)
THEN 2499 dependent_basis=>dependent_field%VARIABLES(field_var_type)%COMPONENTS(mh)%DOMAIN%TOPOLOGY% &
2500 & elements%ELEMENTS(element_number)%BASIS
2501 DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
2503 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)=nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(mhs)* &
2504 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
2505 IF(
ASSOCIATED(rhs_vector))
THEN 2506 IF(
ASSOCIATED(source_field))
THEN 2507 IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
2508 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
2517 CALL flagerror(
"Equations set equations is not associated.",err,error,*999)
2520 CALL flagerror(
"Equations set is not associated.",err,error,*999)
2525 IF(element_number == 1)
THEN 2527 field_variable=>dependent_field%VARIABLES(var1)
2528 DO mh=1,field_variable%NUMBER_OF_COMPONENTS
2529 SELECT CASE(field_variable%COMPONENTS(mh)%INTERPOLATION_TYPE)
2530 CASE(field_node_based_interpolation)
2531 mesh_component_1 = field_variable%COMPONENTS(mh)%MESH_COMPONENT_NUMBER
2532 dependent_basis_1 => dependent_field%DECOMPOSITION%DOMAIN(mesh_component_1)%ptr% &
2533 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
2534 ndofs = ndofs + dependent_basis_1%NUMBER_OF_ELEMENT_PARAMETERS
2536 CASE(field_element_based_interpolation)
2542 &
" is not valid for a finite elasticity equation.",err,error,*999)
2548 & element_number,err,error,*999)
2550 & nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(:), &
2551 &
'(4(X,E13.6))',
'4(4(X,E13.6))',err,error,*999)
2555 exits(
"FiniteElasticity_FiniteElementResidualEvaluate")
2557 999
errors(
"FiniteElasticity_FiniteElementResidualEvaluate",err,error)
2558 exits(
"FiniteElasticity_FiniteElementResidualEvaluate")
2572 INTEGER(INTG),
INTENT(OUT) :: ERR
2578 enters(
"FiniteElasticity_FiniteElementPreResidualEvaluate",err,error,*999)
2580 IF(
ASSOCIATED(equations_set))
THEN 2581 IF(.NOT.
ALLOCATED(equations_set%SPECIFICATION))
THEN 2582 CALL flagerror(
"Equations set specification is not allocated.",err,error,*999)
2583 ELSE IF(
SIZE(equations_set%SPECIFICATION,1)/=3)
THEN 2584 CALL flagerror(
"Equations set specification must have three entries for a finite elasticity type equations set.", &
2587 SELECT CASE(equations_set%SPECIFICATION(3))
2590 dependent_field=>equations_set%EQUATIONS%INTERPOLATION%DEPENDENT_FIELD
2592 & field_u1_variable_type,err,error,*999)
2620 local_error=
"The third equations set specification of "// &
2622 &
" is not valid for a finite elasticity type of an elasticity equation set." 2623 CALL flagerror(local_error,err,error,*999)
2626 CALL flagerror(
"Equations set is not associated.",err,error,*999)
2629 exits(
"FiniteElasticity_FiniteElementPreResidualEvaluate")
2631 999
errors(
"FiniteElasticity_FiniteElementPreResidualEvaluate",err,error)
2632 exits(
"FiniteElasticity_FiniteElementPreResidualEvaluate")
2646 INTEGER(INTG),
INTENT(OUT) :: ERR
2651 enters(
"FiniteElasticity_FiniteElementPostResidualEvaluate",err,error,*999)
2653 IF(
ASSOCIATED(equations_set))
THEN 2654 IF(.NOT.
ALLOCATED(equations_set%SPECIFICATION))
THEN 2655 CALL flagerror(
"Equations set specification is not allocated.",err,error,*999)
2656 ELSE IF(
SIZE(equations_set%SPECIFICATION,1)/=3)
THEN 2657 CALL flagerror(
"Equations set specification must have three entries for a finite elasticity type equations set.", &
2660 SELECT CASE(equations_set%SPECIFICATION(3))
2689 local_error=
"The third equations set specification of "// &
2691 &
" is not valid for a finite elasticity type of an elasticity equation set." 2692 CALL flagerror(local_error,err,error,*999)
2695 CALL flagerror(
"Equations set is not associated.",err,error,*999)
2698 exits(
"FiniteElasticity_FiniteElementPostResidualEvaluate")
2700 999
errors(
"FiniteElasticity_FiniteElementPostResidualEvaluate",err,error)
2701 exits(
"FiniteElasticity_FiniteElementPostResidualEvaluate")
2715 INTEGER(INTG),
INTENT(IN) :: derivedType
2716 INTEGER(INTG),
INTENT(OUT) :: err
2722 enters(
"FiniteElasticityEquationsSet_DerivedVariableCalculate",err,error,*999)
2724 NULLIFY(derivedvariable)
2726 IF(
ASSOCIATED(equationsset))
THEN 2727 IF(.NOT.equationsset%EQUATIONS_SET_FINISHED)
THEN 2728 CALL flagerror(
"Equations set has not been finished.",err,error,*999)
2730 IF(
ASSOCIATED(equationsset%equations))
THEN 2732 SELECT CASE(derivedtype)
2735 & derivedvariable%field,derivedvariable%variable_type,err,error,*999)
2737 CALL flagerror(
"Not implemented.",err,error,*999)
2740 &
" is not valid for a finite elasticity equations set type.",err,error,*999)
2743 CALL flagerror(
"Equations set equations are not associated.",err,error,*999)
2747 CALL flagerror(
"Equations set is not associated.",err,error,*999)
2750 exits(
"FiniteElasticityEquationsSet_DerivedVariableCalculate")
2752 999
errors(
"FiniteElasticityEquationsSet_DerivedVariableCalculate",err,error)
2753 exits(
"FiniteElasticityEquationsSet_DerivedVariableCalculate")
2765 TYPE(
field_type),
POINTER,
INTENT(INOUT) :: strainField
2766 INTEGER(INTG),
INTENT(IN) :: strainFieldVariableType
2767 INTEGER(INTG),
INTENT(OUT) :: err
2772 TYPE(
field_type),
POINTER :: dependentField,geometricField,fibreField
2775 & fibreInterpolationParameters
2781 INTEGER(INTG) :: componentIdx,dependentNumberOfComponents,elementIdx,elementNumber,fieldVariableType,gaussIdx, &
2782 & meshComponentNumber,numberOfComponents,numberOfDimensions,numberOfGauss,numberOfTimes,numberOfXi,partIdx, &
2783 & startIdx,finishIdx
2784 INTEGER(INTG) :: var1
2785 INTEGER(INTG) :: var2
2786 REAL(DP) :: dZdNu(3,3),Fg(3,3),Fe(3,3),J,Jg,Je,C(3,3),f(3,3),E(3,3)
2787 REAL(SP) :: elementUserElapsed,elementSystemElapsed,systemElapsed,systemTime1(1),systemTime2(1),systemTime3(1),systemTime4(1), &
2788 & userElapsed,userTime1(1),userTime2(1),userTime3(1),userTime4(1)
2790 enters(
"FiniteElasticity_StrainCalculate",err,error,*999)
2792 IF(
ASSOCIATED(equationsset))
THEN 2793 equations=>equationsset%equations
2794 IF(
ASSOCIATED(equations))
THEN 2798 IF(
ASSOCIATED(strainfield))
THEN 2799 CALL field_variabletypecheck(strainfield,strainfieldvariabletype,err,error,*999)
2802 numberofcomponents=6
2804 numberofcomponents=3
2806 numberofcomponents=1
2809 &
" is invalid.",err,error,*999)
2811 CALL field_numberofcomponentscheck(strainfield,strainfieldvariabletype,6,err,error,*999)
2812 DO componentidx=1,numberofcomponents
2813 CALL field_componentinterpolationcheck(strainfield,strainfieldvariabletype,componentidx, &
2814 & field_gauss_point_based_interpolation,err,error,*999)
2817 CALL flagerror(
"Strain field is not associated.",err,error,*999)
2823 var1=equations%EQUATIONS_MAPPING%NONLINEAR_MAPPING%RESIDUAL_VARIABLES(1)%PTR%VARIABLE_NUMBER
2824 var2=equations%EQUATIONS_MAPPING%RHS_MAPPING%RHS_VARIABLE%VARIABLE_NUMBER
2826 geometricfield=>equations%interpolation%GEOMETRIC_FIELD
2827 dependentfield=>equations%interpolation%DEPENDENT_FIELD
2828 fibrefield=>equations%interpolation%FIBRE_FIELD
2829 dependentnumberofcomponents=dependentfield%variables(var1)%NUMBER_OF_COMPONENTS
2831 decomposition=>dependentfield%decomposition
2832 meshcomponentnumber=decomposition%MESH_COMPONENT_NUMBER
2835 fieldvariabletype=equations%EQUATIONS_MAPPING%NONLINEAR_MAPPING%RESIDUAL_VARIABLES(1)%PTR%VARIABLE_TYPE
2836 geometricinterpolationparameters=>equations%interpolation%GEOMETRIC_INTERP_PARAMETERS(field_u_variable_type)%ptr
2837 geometricinterpolatedpoint=>equations%interpolation%GEOMETRIC_INTERP_POINT(field_u_variable_type)%ptr
2838 geometricinterpolatedpointmetrics=>equations%interpolation%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%ptr
2839 dependentinterpolationparameters=>equations%interpolation%DEPENDENT_INTERP_PARAMETERS(fieldvariabletype)%ptr
2840 dependentinterpolatedpoint=>equations%interpolation%DEPENDENT_INTERP_POINT(fieldvariabletype)%ptr
2841 dependentinterpolatedpointmetrics=>equations%interpolation%DEPENDENT_INTERP_POINT_METRICS(fieldvariabletype)%ptr
2842 IF(
ASSOCIATED(fibrefield))
THEN 2843 fibreinterpolationparameters=>equations%interpolation%FIBRE_INTERP_PARAMETERS(field_u_variable_type)%ptr
2844 fibreinterpolatedpoint=>equations%interpolation%FIBRE_INTERP_POINT(field_u_variable_type)%ptr
2846 NULLIFY(fibreinterpolationparameters)
2847 NULLIFY(fibreinterpolatedpoint)
2850 elementsmapping=>dependentfield%decomposition%domain(meshcomponentnumber)%ptr%mappings%elements
2863 startidx=elementsmapping%BOUNDARY_START
2864 finishidx=elementsmapping%GHOST_FINISH
2866 startidx=elementsmapping%INTERNAL_START
2867 finishidx=elementsmapping%INTERNAL_FINISH
2871 DO elementidx=startidx,finishidx
2873 numberoftimes=numberoftimes+1
2874 elementnumber=elementsmapping%DOMAIN_LIST(elementidx)
2880 dependentbasis=>decomposition%domain(meshcomponentnumber)%ptr%topology%elements%elements(elementnumber)%basis
2882 numberofgauss=dependentquadraturescheme%NUMBER_OF_GAUSS
2884 numberofxi=dependentbasis%NUMBER_OF_XI
2886 CALL field_interpolationparameterselementget(field_values_set_type,elementnumber,geometricinterpolationparameters, &
2888 CALL field_interpolationparameterselementget(field_values_set_type,elementnumber,dependentinterpolationparameters, &
2890 IF(
ASSOCIATED(fibrefield))
THEN 2891 CALL field_interpolationparameterselementget(field_values_set_type,elementnumber,fibreinterpolationparameters, &
2896 DO gaussidx=1,numberofgauss
2905 CALL field_interpolatedpointmetricscalculate(numberofxi,dependentinterpolatedpointmetrics,err,error,*999)
2908 CALL field_interpolatedpointmetricscalculate(numberofxi,geometricinterpolatedpointmetrics,err,error,*999)
2909 IF(
ASSOCIATED(fibrefield))
THEN 2916 & geometricinterpolatedpointmetrics,fibreinterpolatedpoint,dzdnu,err,error,*999)
2919 & dzdnu,fg,fe,jg,je,err,error,*999)
2928 CALL field_parametersetupdatelocalgausspoint(strainfield,strainfieldvariabletype,field_values_set_type, &
2929 & gaussidx,elementnumber,1,c(1,1),err,error,*999)
2930 CALL field_parametersetupdatelocalgausspoint(strainfield,strainfieldvariabletype,field_values_set_type, &
2931 & gaussidx,elementnumber,2,c(1,2),err,error,*999)
2932 CALL field_parametersetupdatelocalgausspoint(strainfield,strainfieldvariabletype,field_values_set_type, &
2933 & gaussidx,elementnumber,3,c(1,3),err,error,*999)
2934 CALL field_parametersetupdatelocalgausspoint(strainfield,strainfieldvariabletype,field_values_set_type, &
2935 & gaussidx,elementnumber,4,c(2,2),err,error,*999)
2936 CALL field_parametersetupdatelocalgausspoint(strainfield,strainfieldvariabletype,field_values_set_type, &
2937 & gaussidx,elementnumber,5,c(2,3),err,error,*999)
2938 CALL field_parametersetupdatelocalgausspoint(strainfield,strainfieldvariabletype,field_values_set_type, &
2939 & gaussidx,elementnumber,6,c(3,3),err,error,*999)
2943 CALL field_parametersetupdatelocalgausspoint(strainfield,strainfieldvariabletype,field_values_set_type, &
2944 & gaussidx,elementnumber,1,c(1,1),err,error,*999)
2945 CALL field_parametersetupdatelocalgausspoint(strainfield,strainfieldvariabletype,field_values_set_type, &
2946 & gaussidx,elementnumber,2,c(1,2),err,error,*999)
2947 CALL field_parametersetupdatelocalgausspoint(strainfield,strainfieldvariabletype,field_values_set_type, &
2948 & gaussidx,elementnumber,3,c(2,2),err,error,*999)
2951 CALL field_parametersetupdatelocalgausspoint(strainfield,strainfieldvariabletype,field_values_set_type, &
2952 & gaussidx,elementnumber,1,c(1,1),err,error,*999)
2956 CALL flagerror(localerror,err,error,*999)
2965 userelapsed=usertime2(1)-usertime1(1)
2966 systemelapsed=systemtime2(1)-systemtime1(1)
2967 elementuserelapsed=elementuserelapsed+userelapsed
2968 elementsystemelapsed=elementsystemelapsed+systemelapsed
2971 & userelapsed,err,error,*999)
2973 & systemelapsed,err,error,*999)
2976 & userelapsed,err,error,*999)
2978 & systemelapsed,err,error,*999)
2979 IF(numberoftimes>0)
THEN 2981 & elementuserelapsed/numberoftimes,err,error,*999)
2983 & elementsystemelapsed/numberoftimes,err,error,*999)
2994 CALL field_parametersetupdatestart(strainfield,strainfieldvariabletype,field_values_set_type,err,error,*999)
2997 CALL field_parametersetupdatefinish(strainfield,strainfieldvariabletype,field_values_set_type,err,error,*999)
3002 userelapsed=usertime4(1)-usertime3(1)
3003 systemelapsed=systemtime4(1)-systemtime3(1)
3014 CALL flagerror(
"Equations set equations is not associated.",err,error,*999)
3017 CALL flagerror(
"Equations set is not associated.",err,error,*999)
3020 exits(
"FiniteElasticity_StrainCalculate")
3022 999 errorsexits(
"FiniteElasticity_StrainCalculate",err,error)
3035 INTEGER(INTG),
INTENT(IN) :: tensorEvaluateType
3036 INTEGER(INTG),
INTENT(IN) :: userElementNumber
3037 REAL(DP),
INTENT(IN) :: xi(:)
3038 REAL(DP),
INTENT(OUT) :: values(3,3)
3039 INTEGER(INTG),
INTENT(OUT) :: err
3045 & fibreInterpolatedPoint,dependentInterpolatedPoint,materialsInterpolatedPoint, &
3046 & independentInterpolatedPoint,darcyInterpolatedPoint
3048 & dependentInterpolatedPointMetrics
3054 LOGICAL :: userElementExists,ghostElement
3055 INTEGER(INTG) :: dependentVarType,meshComponentNumber
3056 INTEGER(INTG) :: numberOfDimensions,numberOfXi
3057 INTEGER(INTG) :: localElementNumber,i,nh,mh
3058 REAL(DP) :: dZdNu(3,3),dZdNuT(3,3),AZL(3,3),E(3,3),cauchyStressTensor(3,3),cauchyStressVoigt(6),Jznu
3060 enters(
"FiniteElasticity_TensorInterpolateXi",err,error,*999)
3063 NULLIFY(dependentfield)
3064 NULLIFY(geometricinterpolatedpoint)
3065 NULLIFY(fibreinterpolatedpoint)
3066 NULLIFY(dependentinterpolatedpoint)
3067 NULLIFY(materialsinterpolatedpoint)
3068 NULLIFY(independentinterpolatedpoint)
3069 NULLIFY(darcyinterpolatedpoint)
3070 NULLIFY(decomposition)
3071 NULLIFY(decompositiontopology)
3072 NULLIFY(domaintopology)
3073 NULLIFY(elementbasis)
3075 IF(.NOT.
ASSOCIATED(equationsset))
THEN 3076 CALL flagerror(
"Equations set is not associated.",err,error,*999)
3078 equations=>equationsset%equations
3079 IF(.NOT.
ASSOCIATED(equations))
THEN 3080 CALL flagerror(
"Equations set equations is not associated.",err,error,*999)
3083 nonlinearmapping=>equations%equations_mapping%nonlinear_mapping
3084 IF(.NOT.
ASSOCIATED(equations))
THEN 3085 CALL flagerror(
"Equations nonlinear mapping is not associated.",err,error,*999)
3087 dependentvartype=nonlinearmapping%residual_variables(1)%ptr%variable_type
3089 IF(.NOT.
ASSOCIATED(equations%interpolation))
THEN 3090 CALL flagerror(
"Equations interpolation is not associated.",err,error,*999)
3092 dependentfield=>equations%interpolation%dependent_field
3093 IF(.NOT.
ASSOCIATED(dependentfield))
THEN 3094 CALL flagerror(
"Equations dependent field is not associated.",err,error,*999)
3096 decomposition=>dependentfield%decomposition
3097 IF(.NOT.
ASSOCIATED(decomposition))
THEN 3098 CALL flagerror(
"Dependent field decomposition is not associated.",err,error,*999)
3100 CALL decomposition_mesh_component_number_get(decomposition,meshcomponentnumber,err,error,*999)
3101 decompositiontopology=>decomposition%topology
3102 domaintopology=>decomposition%domain(meshcomponentnumber)%ptr%topology
3103 CALL decomposition_topology_element_check_exists(decompositiontopology,userelementnumber, &
3104 & userelementexists,localelementnumber,ghostelement,err,error,*999)
3105 IF(.NOT.userelementexists)
THEN 3106 CALL flagerror(
"The specified user element number of "// &
3108 &
" does not exist in the decomposition for the dependent field.",err,error,*999)
3110 CALL domaintopology_elementbasisget( &
3111 & domaintopology,userelementnumber,elementbasis,err,error,*999)
3114 CALL field_interpolation_parameters_element_get(field_values_set_type,localelementnumber, &
3115 & equations%interpolation%geometric_interp_parameters(field_u_variable_type)%ptr,err,error,*999)
3116 IF(
ASSOCIATED(equations%interpolation%fibre_interp_parameters))
THEN 3117 CALL field_interpolation_parameters_element_get(field_values_set_type,localelementnumber, &
3118 & equations%interpolation%fibre_interp_parameters(field_u_variable_type)%ptr,err,error,*999)
3120 CALL field_interpolation_parameters_element_get(field_values_set_type,localelementnumber, &
3121 & equations%interpolation%dependent_interp_parameters(dependentvartype)%ptr,err,error,*999)
3124 geometricinterpolatedpoint=>equations%interpolation%geometric_interp_point(field_u_variable_type)%ptr
3125 IF(
ASSOCIATED(equations%interpolation%fibre_interp_point))
THEN 3126 fibreinterpolatedpoint=>equations%interpolation%fibre_interp_point(field_u_variable_type)%ptr
3128 dependentinterpolatedpoint=>equations%interpolation%dependent_interp_point(dependentvartype)%ptr
3131 geometricinterpolatedpointmetrics=>equations%interpolation% &
3132 & geometric_interp_point_metrics(field_u_variable_type)%ptr
3133 dependentinterpolatedpointmetrics=>equations%interpolation% &
3134 & dependent_interp_point_metrics(dependentvartype)%ptr
3137 CALL field_interpolate_xi(
first_part_deriv,xi,dependentinterpolatedpoint,err,error,*999)
3138 CALL field_interpolate_xi(
first_part_deriv,xi,geometricinterpolatedpoint,err,error,*999)
3139 IF(
ASSOCIATED(fibreinterpolatedpoint))
THEN 3140 CALL field_interpolate_xi(
first_part_deriv,xi,fibreinterpolatedpoint,err,error,*999)
3144 CALL field_interpolated_point_metrics_calculate( &
3145 & elementbasis%number_of_xi,geometricinterpolatedpointmetrics,err,error,*999)
3146 CALL field_interpolated_point_metrics_calculate( &
3147 & elementbasis%number_of_xi,dependentinterpolatedpointmetrics,err,error,*999)
3150 numberofdimensions=equationsset%region%coordinate_system%number_of_dimensions
3151 numberofxi=elementbasis%number_of_xi
3153 & geometricinterpolatedpointmetrics,fibreinterpolatedpoint,dzdnu,err,error,*999)
3165 e(i,i)=e(i,i)-0.5_dp
3171 CALL field_interpolation_parameters_element_get(field_values_set_type,localelementnumber, &
3172 & equations%interpolation%materials_interp_parameters(field_u_variable_type)%ptr,err,error,*999)
3173 IF(
ASSOCIATED(equations%interpolation%independent_interp_parameters))
THEN 3174 CALL field_interpolation_parameters_element_get(field_values_set_type,localelementnumber, &
3175 & equations%interpolation%independent_interp_parameters(field_u_variable_type)%ptr,err,error,*999)
3179 materialsinterpolatedpoint=>equations%interpolation%materials_interp_point(field_u_variable_type)%ptr
3180 IF(
ASSOCIATED(equations%interpolation%independent_interp_point))
THEN 3181 independentinterpolatedpoint=>equations%interpolation%independent_interp_point(dependentvartype)%ptr
3185 CALL field_interpolate_xi(
no_part_deriv,xi,materialsinterpolatedpoint,err,error,*999)
3186 IF(
ASSOCIATED(independentinterpolatedpoint))
THEN 3187 CALL field_interpolate_xi(
first_part_deriv,xi,independentinterpolatedpoint,err,error,*999)
3190 SELECT CASE(equationsset%specification(3))
3193 jznu=dependentinterpolatedpointmetrics%JACOBIAN/geometricinterpolatedpointmetrics%JACOBIAN
3201 & materialsinterpolatedpoint,cauchystressvoigt,dzdnu,jznu,localelementnumber,0,err,error,*999)
3211 & materialsinterpolatedpoint,darcyinterpolatedpoint, &
3212 & independentinterpolatedpoint,cauchystresstensor,jznu,dzdnu,localelementnumber,0,err,error,*999)
3214 CALL flagerror(
"Not implemented ",err,error,*999)
3218 SELECT CASE(tensorevaluatetype)
3226 values=cauchystresstensor
3228 CALL flagerror(
"Not implemented.",err,error,*999)
3231 &
"for finite elasticity equation sets",err,error,*999)
3234 exits(
"FiniteElasticity_TensorInterpolateXi")
3236 999 errorsexits(
"FiniteElasticity_TensorInterpolateXi",err,error)
3252 INTEGER(INTG),
INTENT(IN) :: ELEMENT_NUMBER
3253 INTEGER(INTG),
INTENT(OUT) :: ERR
3275 INTEGER(INTG) :: FACE_NUMBER,xiDirection(3),orientation
3276 INTEGER(INTG) :: FIELD_VAR_U_TYPE,FIELD_VAR_DELUDELN_TYPE,MESH_COMPONENT_NUMBER
3277 INTEGER(INTG) :: oh,mh,ms,mhs,nh,ns,nhs,ng,naf
3278 INTEGER(INTG) :: NUMBER_OF_DIMENSIONS,NUMBER_OF_LOCAL_FACES
3279 INTEGER(INTG) :: SUM_ELEMENT_PARAMETERS
3280 INTEGER(INTG) :: ELEMENT_BASE_DOF_INDEX(3),NUMBER_OF_FACE_PARAMETERS(3)
3281 INTEGER(INTG),
PARAMETER :: OFF_DIAG_COMP(3)=[0,1,3],off_diag_dep_var1(3)=[1,1,2],off_diag_dep_var2(3)=[2,3,3]
3282 REAL(DP) :: PRESSURE_GAUSS,GW_PRESSURE
3283 REAL(DP) :: NORMAL(3),GW_PRESSURE_W(2),TEMP3, TEMP4
3284 REAL(DP) :: TEMPVEC1(2),TEMPVEC2(2),TEMPVEC3(3),TEMPVEC4(3),TEMPVEC5(3)
3285 LOGICAL :: NONZERO_PRESSURE
3287 enters(
"FiniteElasticity_SurfacePressureJacobianEvaluate",err,error,*999)
3289 NULLIFY(dependent_basis)
3290 NULLIFY(decomposition)
3292 NULLIFY(equations,equations_mapping,equations_matrices,nonlinear_mapping,nonlinear_matrices,jacobian_matrix)
3293 NULLIFY(dependent_interpolation_parameters,pressure_interpolation_parameters)
3294 NULLIFY(dependent_interp_point,dependent_interp_point_metrics,pressure_interp_point)
3295 NULLIFY(dependent_field)
3296 NULLIFY(field_variable)
3297 NULLIFY(dependent_quadrature_scheme)
3299 number_of_dimensions=equations_set%REGION%COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS
3301 equations=>equations_set%EQUATIONS
3302 equations_matrices=>equations%EQUATIONS_MATRICES
3303 nonlinear_matrices=>equations_matrices%NONLINEAR_MATRICES
3304 jacobian_matrix=>nonlinear_matrices%JACOBIANS(1)%PTR
3306 dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
3307 decomposition=>dependent_field%DECOMPOSITION
3308 mesh_component_number=decomposition%MESH_COMPONENT_NUMBER
3309 element=>decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)
3310 number_of_local_faces=dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR% &
3311 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS%NUMBER_OF_LOCAL_FACES
3313 field_variable=>equations%EQUATIONS_MAPPING%NONLINEAR_MAPPING%RESIDUAL_VARIABLES(1)%PTR
3314 field_var_u_type=equations%EQUATIONS_MAPPING%NONLINEAR_MAPPING%RESIDUAL_VARIABLES(1)%PTR%VARIABLE_TYPE
3315 field_var_deludeln_type=equations%EQUATIONS_MAPPING%RHS_MAPPING%RHS_VARIABLE_TYPE
3318 DO naf=1,number_of_local_faces
3319 face_number=element%ELEMENT_FACES(naf)
3320 face=>decomposition%TOPOLOGY%FACES%FACES(face_number)
3323 IF(face%BOUNDARY_FACE)
THEN 3324 xidirection(3)=abs(face%XI_DIRECTION)
3326 pressure_interpolation_parameters=>equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_deludeln_type)%PTR
3327 CALL field_interpolation_parameters_face_get(field_pressure_values_set_type,face_number, &
3328 & pressure_interpolation_parameters,err,error,*999,field_geometric_components_type)
3329 pressure_interp_point=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_deludeln_type)%PTR
3332 nonzero_pressure=any(abs(pressure_interpolation_parameters%PARAMETERS(:,xidirection(3)))>
zero_tolerance)
3335 IF(nonzero_pressure)
THEN 3336 mesh_component_number=decomposition%MESH_COMPONENT_NUMBER
3337 dependent_basis=>decomposition%DOMAIN(mesh_component_number)%PTR%TOPOLOGY%FACES%FACES(face_number)%BASIS
3340 dependent_interpolation_parameters=>equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_u_type)%PTR
3341 CALL field_interpolation_parameters_face_get(field_values_set_type,face_number, &
3342 & dependent_interpolation_parameters,err,error,*999,field_geometric_components_type)
3343 dependent_interp_point=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_u_type)%PTR
3344 dependent_interp_point_metrics=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT_METRICS(field_var_u_type)%PTR
3346 sum_element_parameters=0
3348 DO nh=1,number_of_dimensions
3349 mesh_component_number=field_variable%COMPONENTS(nh)%MESH_COMPONENT_NUMBER
3350 dependent_basis=>decomposition%DOMAIN(mesh_component_number)%PTR%TOPOLOGY%FACES%FACES(face_number)%BASIS
3351 bases(nh)%PTR=>decomposition%DOMAIN(mesh_component_number)%PTR% &
3352 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3354 number_of_face_parameters(nh)=dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
3355 element_base_dof_index(nh)=sum_element_parameters
3356 sum_element_parameters=sum_element_parameters+bases(nh)%PTR%NUMBER_OF_ELEMENT_PARAMETERS
3364 DO ng=1,dependent_quadrature_scheme%NUMBER_OF_GAUSS
3366 & pressure_interp_point,err,error,*999,field_geometric_components_type)
3368 & dependent_interp_point,err,error,*999,field_geometric_components_type)
3370 & dependent_interp_point_metrics,err,error,*999)
3372 CALL cross_product(dependent_interp_point_metrics%DX_DXI(:,1), &
3373 & dependent_interp_point_metrics%DX_DXI(:,2),normal,err,error,*999)
3374 pressure_gauss=pressure_interp_point%VALUES(xidirection(3),
no_part_deriv)*orientation
3375 gw_pressure=dependent_quadrature_scheme%GAUSS_WEIGHTS(ng)*pressure_gauss
3377 DO oh=1,off_diag_comp(number_of_dimensions)
3378 nh=off_diag_dep_var1(oh)
3379 mh=off_diag_dep_var2(oh)
3380 gw_pressure_w(1:2)=(normal(mh)*dependent_interp_point_metrics%DXI_DX(1:2,nh)- &
3381 & dependent_interp_point_metrics%DXI_DX(1:2,mh)*normal(nh))*gw_pressure
3382 DO ns=1,number_of_face_parameters(nh)
3384 nhs=element_base_dof_index(nh)+ &
3385 & bases(nh)%PTR%ELEMENT_PARAMETERS_IN_LOCAL_FACE(ns,naf)
3386 tempvec1(1:2)=gw_pressure_w(1:2)*quadrature_schemes(nh)%PTR% &
3388 DO ms=1,number_of_face_parameters(mh)
3389 mhs=element_base_dof_index(mh)+ &
3390 & bases(mh)%PTR%ELEMENT_PARAMETERS_IN_LOCAL_FACE(ms,naf)
3391 tempvec2=quadrature_schemes(mh)%PTR%GAUSS_BASIS_FNS(ms,
no_part_deriv,ng)
3392 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(mhs,nhs)+ &
3393 & dot_product(tempvec1,tempvec2)* &
3394 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_u_type)%PTR%SCALE_FACTORS(ms,mh)* &
3395 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_u_type)%PTR%SCALE_FACTORS(ns,nh)
3400 DO oh=1,off_diag_comp(number_of_dimensions)
3401 nh=off_diag_dep_var1(oh)
3402 mh=off_diag_dep_var2(oh)
3403 gw_pressure_w(1:2)=(normal(nh)*dependent_interp_point_metrics%DXI_DX(1:2,mh)- &
3404 & dependent_interp_point_metrics%DXI_DX(1:2,nh)*normal(mh))*gw_pressure
3405 DO ms=1,number_of_face_parameters(mh)
3407 mhs=element_base_dof_index(mh)+ &
3408 & bases(mh)%PTR%ELEMENT_PARAMETERS_IN_LOCAL_FACE(ms,naf)
3409 tempvec1(1:2)=gw_pressure_w(1:2)*quadrature_schemes(mh)%PTR% &
3411 DO ns=1,number_of_face_parameters(nh)
3412 nhs=element_base_dof_index(nh)+ &
3413 & bases(nh)%PTR%ELEMENT_PARAMETERS_IN_LOCAL_FACE(ns,naf)
3414 tempvec2=quadrature_schemes(nh)%PTR%GAUSS_BASIS_FNS(ns,
no_part_deriv,ng)
3415 jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(nhs,mhs)=jacobian_matrix%ELEMENT_JACOBIAN%MATRIX(nhs,mhs)+ &
3416 & dot_product(tempvec1,tempvec2)* &
3417 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_u_type)%PTR%SCALE_FACTORS(ms,mh)* &
3418 & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_u_type)%PTR%SCALE_FACTORS(ns,nh)
3427 exits(
"FiniteElasticity_SurfacePressureJacobianEvaluate")
3429 999
errors(
"FiniteElasticity_SurfacePressureJacobianEvaluate",err,error)
3430 exits(
"FiniteElasticity_SurfacePressureJacobianEvaluate")
3443 INTEGER(INTG),
INTENT(IN) :: ELEMENT_NUMBER
3444 INTEGER(INTG),
INTENT(IN) :: var1
3445 INTEGER(INTG),
INTENT(IN) :: var2
3446 INTEGER(INTG),
INTENT(OUT) :: ERR
3449 TYPE(
basis_type),
POINTER :: DEPENDENT_FACE_BASIS,COMPONENT_FACE_BASIS,COMPONENT_BASIS
3463 INTEGER(INTG) :: FIELD_VAR_U_TYPE,FIELD_VAR_DUDN_TYPE,MESH_COMPONENT_NUMBER
3464 INTEGER(INTG) :: element_face_idx,face_number,gauss_idx
3465 INTEGER(INTG) :: component_idx,element_base_dof_idx,element_dof_idx,parameter_idx,face_parameter_idx
3466 INTEGER(INTG) :: NUMBER_OF_DIMENSIONS,NUMBER_OF_LOCAL_FACES
3467 INTEGER(INTG) :: xiDirection(3),orientation
3468 REAL(DP) :: PRESSURE_GAUSS,GW_PRESSURE,GW_PRESSURE_NORMAL_COMPONENT
3469 REAL(DP) :: NORMAL(3)
3470 LOGICAL :: NONZERO_PRESSURE
3472 enters(
"FiniteElasticity_SurfacePressureResidualEvaluate",err,error,*999)
3474 NULLIFY(dependent_face_basis,component_face_basis,component_basis)
3475 NULLIFY(decomposition)
3476 NULLIFY(decomp_element)
3477 NULLIFY(decomp_face)
3479 NULLIFY(equations,nonlinear_matrices)
3480 NULLIFY(dependent_field,field_variable)
3481 NULLIFY(face_dependent_interpolation_parameters)
3482 NULLIFY(face_dependent_interpolated_point,face_dependent_interpolated_point_metrics)
3483 NULLIFY(face_pressure_interpolation_parameters,face_pressure_interpolated_point)
3484 NULLIFY(component_face_quadrature_scheme,face_quadrature_scheme)
3486 number_of_dimensions=equations_set%REGION%COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS
3489 equations=>equations_set%EQUATIONS
3490 nonlinear_matrices=>equations%EQUATIONS_MATRICES%NONLINEAR_MATRICES
3491 dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
3492 decomposition=>dependent_field%DECOMPOSITION
3493 mesh_component_number=decomposition%MESH_COMPONENT_NUMBER
3494 decomp_element=>decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(element_number)
3497 field_variable=>equations%EQUATIONS_MAPPING%NONLINEAR_MAPPING%RESIDUAL_VARIABLES(1)%PTR
3498 field_var_u_type=field_variable%VARIABLE_TYPE
3499 field_var_dudn_type=equations%EQUATIONS_MAPPING%RHS_MAPPING%RHS_VARIABLE_TYPE
3500 number_of_local_faces=decomposition%DOMAIN(mesh_component_number)%PTR%TOPOLOGY%ELEMENTS% &
3501 & elements(element_number)%BASIS%NUMBER_OF_LOCAL_FACES
3504 DO element_face_idx=1,number_of_local_faces
3505 face_number=decomp_element%ELEMENT_FACES(element_face_idx)
3506 decomp_face=>decomposition%TOPOLOGY%FACES%FACES(face_number)
3509 IF(decomp_face%BOUNDARY_FACE)
THEN 3510 xidirection(3)=abs(decomp_face%XI_DIRECTION)
3512 face_pressure_interpolation_parameters=>equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_dudn_type)%PTR
3513 CALL field_interpolation_parameters_face_get(field_pressure_values_set_type,face_number, &
3514 & face_pressure_interpolation_parameters,err,error,*999,field_geometric_components_type)
3515 face_pressure_interpolated_point=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT(var2)%PTR
3518 nonzero_pressure=any(abs(face_pressure_interpolation_parameters%PARAMETERS(:,xidirection(3)))>
zero_tolerance)
3521 IF(nonzero_pressure)
THEN 3523 dependent_face_basis=>decomposition%DOMAIN(mesh_component_number)%PTR%TOPOLOGY%FACES%FACES(face_number)%BASIS
3526 face_dependent_interpolation_parameters=>equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_u_type)%PTR
3527 CALL field_interpolation_parameters_face_get(field_values_set_type,face_number, &
3528 & face_dependent_interpolation_parameters,err,error,*999,field_geometric_components_type)
3529 face_dependent_interpolated_point=>equations%INTERPOLATION%DEPENDENT_INTERP_POINT(field_var_u_type)%PTR
3530 face_dependent_interpolated_point_metrics=>equations%INTERPOLATION% &
3531 & dependent_interp_point_metrics(field_var_u_type)%PTR
3540 DO gauss_idx=1,face_quadrature_scheme%NUMBER_OF_GAUSS
3543 & face_pressure_interpolated_point,err,error,*999,field_geometric_components_type)
3545 & face_dependent_interpolated_point,err,error,*999)
3547 & face_dependent_interpolated_point_metrics,err,error,*999)
3549 CALL cross_product(face_dependent_interpolated_point_metrics%DX_DXI(:,1), &
3550 & face_dependent_interpolated_point_metrics%DX_DXI(:,2),normal,err,error,*999)
3551 pressure_gauss=face_pressure_interpolated_point%VALUES(xidirection(3),
no_part_deriv)*orientation
3552 gw_pressure=face_quadrature_scheme%GAUSS_WEIGHTS(gauss_idx)*pressure_gauss
3553 element_base_dof_idx=0
3555 DO component_idx=1,number_of_dimensions
3556 mesh_component_number=field_variable%COMPONENTS(component_idx)%MESH_COMPONENT_NUMBER
3557 component_basis=>dependent_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR% &
3558 & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
3559 component_face_basis=>decomposition%DOMAIN(mesh_component_number)%PTR%TOPOLOGY%FACES%FACES(face_number)%BASIS
3560 component_face_quadrature_scheme=>component_face_basis% &
3562 gw_pressure_normal_component=gw_pressure*normal(component_idx)
3563 DO face_parameter_idx=1,component_face_basis%NUMBER_OF_ELEMENT_PARAMETERS
3564 parameter_idx=component_basis%ELEMENT_PARAMETERS_IN_LOCAL_FACE(face_parameter_idx,element_face_idx)
3565 element_dof_idx=element_base_dof_idx+parameter_idx
3566 nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)= &
3567 & nonlinear_matrices%ELEMENT_RESIDUAL%VECTOR(element_dof_idx)+ &
3568 & gw_pressure_normal_component * &
3569 & component_face_quadrature_scheme%GAUSS_BASIS_FNS(face_parameter_idx,
no_part_deriv,gauss_idx)
3572 element_base_dof_idx=element_base_dof_idx+component_basis%NUMBER_OF_ELEMENT_PARAMETERS
3579 exits(
"FiniteElasticity_SurfacePressureResidualEvaluate")
3581 999
errors(
"FiniteElasticity_SurfacePressureResidualEvaluate",err,error)
3582 exits(
"FiniteElasticity_SurfacePressureResidualEvaluate")
3593 & fibreinterpolatedpoint,dzdnu,err,error,*)
3598 REAL(DP),
INTENT(OUT) :: dZdNu(3,3)
3599 INTEGER(INTG),
INTENT(OUT) :: err
3602 INTEGER(INTG) :: numberOfXDimensions,numberOfXiDimensions,numberOfZDimensions
3603 REAL(DP) :: dNuDXi(3,3),dXidNu(3,3), dNudX(3,3),dXdNu(3,3)
3605 enters(
"FiniteElasticity_GaussDeformationGradientTensor",err,error,*999)
3607 IF(
ASSOCIATED(dependentinterppointmetrics))
THEN 3608 IF(
ASSOCIATED(geometricinterppointmetrics))
THEN 3609 numberofxdimensions=geometricinterppointmetrics%NUMBER_OF_X_DIMENSIONS
3610 numberofxidimensions=geometricinterppointmetrics%NUMBER_OF_XI_DIMENSIONS
3611 numberofzdimensions=dependentinterppointmetrics%NUMBER_OF_X_DIMENSIONS
3614 & dnudxi(1:numberofxdimensions,1:numberofxidimensions), &
3615 & dxidnu(1:numberofxidimensions,1:numberofxdimensions),err,error,*999)
3617 CALL matrixproduct(dependentinterppointmetrics%DX_DXI(1:numberofzdimensions,1:numberofxidimensions), &
3618 & dxidnu(1:numberofxidimensions,1:numberofxdimensions),dzdnu(1:numberofzdimensions,1:numberofxdimensions), &
3621 IF(numberofzdimensions == 2)
THEN 3622 dzdnu(:,3) = [0.0_dp,0.0_dp,1.0_dp]
3623 dzdnu(3,1:2) = 0.0_dp
3633 &
'(" dZdNu',
'(",I1,",:)',
' :",3(X,E13.6))',
'(15X,3(X,E13.6))',err,error,*999)
3637 CALL flagerror(
"Geometric interpolated point metrics is not associated.",err,error,*999)
3640 CALL flagerror(
"Dependent interpolated point metrics is not associated.",err,error,*999)
3643 exits(
"FiniteElasticity_GaussDeformationGradientTensor")
3645 999
errors(
"FiniteElasticity_GaussDeformationGradientTensor",err,error)
3646 exits(
"FiniteElasticity_GaussDeformationGradientTensor")
3657 & materials_interpolated_point,darcy_dependent_interpolated_point, &
3658 & independent_interpolated_point,cauchy_tensor,jznu,dzdnu,
element_number,gauss_point_number,err,error,*)
3665 REAL(DP),
INTENT(OUT) :: CAUCHY_TENSOR(:,:)
3666 REAL(DP),
INTENT(OUT) :: Jznu
3667 REAL(DP),
INTENT(IN) :: DZDNU(3,3)
3668 INTEGER(INTG),
INTENT(IN) :: ELEMENT_NUMBER,GAUSS_POINT_NUMBER
3669 INTEGER(INTG),
INTENT(OUT) :: ERR
3672 INTEGER(INTG) :: EQUATIONS_SET_SUBTYPE
3673 INTEGER(INTG) :: i,j,k,PRESSURE_COMPONENT,component_idx,dof_idx
3674 REAL(DP) :: activation
3675 REAL(DP) :: AZL(3,3),AZU(3,3),DZDNUT(3,3),PIOLA_TENSOR(3,3),E(3,3),P,IDENTITY(3,3),AZLT(3,3),AZUT(3,3)
3676 REAL(DP) :: AZL_SQUARED(3,3)
3677 REAL(DP) :: I1,I2,I3
3678 REAL(DP) :: ACTIVE_STRESS_11,ACTIVE_STRESS_22,ACTIVE_STRESS_33
3679 REAL(DP) :: TEMP(3,3),TEMPTERM
3682 REAL(DP),
DIMENSION (:),
POINTER :: C
3683 REAL(DP) :: a, B(3,3), Q
3684 REAL(DP) :: ffact,dfdJfact
3685 INTEGER(INTG) :: DARCY_MASS_INCREASE_ENTRY
3686 REAL(DP) ::
VALUE,VAL1,VAL2
3687 REAL(DP) :: WV_PRIME,TOL,TOL1,UP,LOW
3688 REAL(DP) :: F_e(3,3),F_a(3,3),F_a_inv(3,3),F_a_T(3,3),C_a(3,3),C_a_inv(3,3),lambda_a,C_e(3,3),F_e_T(3,3)
3689 REAL(DP) :: REFERENCE_VOLUME,XB_STIFFNESS,XB_DISTORTION,V_MAX
3690 REAL(DP) :: SARCO_LENGTH,FREE_ENERGY,FREE_ENERGY_0,XB_ENERGY_PER_VOLUME,SLOPE,lambda_f,A_1,A_2,x_1,x_2
3691 REAL(DP) :: MAX_XB_NUMBER_PER_VOLUME,ENERGY_PER_XB,FORCE_LENGTH,I_1e,EVALUES(3),EVECTOR_1(3),EVECTOR_2(3),EVECTOR_3(3)
3692 REAL(DP) :: EMATRIX_1(3,3),EMATRIX_2(3,3),EMATRIX_3(3,3),TEMP1(3,3),TEMP2(3,3),TEMP3(3,3),N1(3,3),N2(3,3),N3(3,3)
3693 REAL(DP),
DIMENSION(5) :: PAR
3694 INTEGER(INTG) :: LWORK,node1,node2
3695 INTEGER(INTG),
PARAMETER :: LWMAX=1000
3696 REAL(DP) :: WORK(lwmax),RIGHT_NODE(3),LEFT_NODE(3),delta_t,dist1,dist2,velo
3697 TYPE(
field_type),
POINTER :: DEPENDENT_FIELD,INDEPENDENT_FIELD
3698 REAL(DP) :: ISOMETRIC_FORCE_AT_FULL_ACT,LENGTH_HALF_SARCO
3699 REAL(DP) :: TITIN_VALUE,TITIN_VALUE_CROSS_FIBRE,TITIN_UNBOUND,TITIN_BOUND
3700 REAL(DP) :: TITIN_UNBOUND_CROSS_FIBRE,TITIN_BOUND_CROSS_FIBRE
3702 enters(
"FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR",err,error,*999)
3704 NULLIFY(field_variable)
3706 IF(.NOT.
ALLOCATED(equations_set%SPECIFICATION))
THEN 3707 CALL flagerror(
"Equations set specification is not allocated.",err,error,*999)
3708 ELSE IF(
SIZE(equations_set%SPECIFICATION,1)/=3)
THEN 3709 CALL flagerror(
"Equations set specification must have three entries for a finite elasticity type equations set.", &
3712 equations_set_subtype=equations_set%SPECIFICATION(3)
3713 c=>materials_interpolated_point%VALUES(:,1)
3725 pressure_component=dependent_interpolated_point%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
3726 p=dependent_interpolated_point%VALUES(pressure_component,1)
3728 CALL invert(azl,azu,i3,err,error,*999)
3732 e(i,i)=e(i,i)-0.5_dp
3736 & 3,3,e,
write_string_matrix_name_and_indices,
'(" E',
'(",I1,",:)',
' :",3(X,E13.6))', &
3737 &
'(17X,3(X,E13.6))',err,error,*999)
3741 identity(i,i)=1.0_dp
3744 SELECT CASE(equations_set_subtype)
3749 wv_prime = c(3)*(jznu - 1.0_dp)
3751 i1 = azl(1,1) + azl(2,2) + azl(3,3)
3753 i2 = 0.5_dp * (
i1**2 - azl_squared(1,1) - azl_squared(2,2) - azl_squared(3,3))
3755 piola_tensor=2.0_dp*jznu**(-2.0_dp/3.0_dp)*((c(1)+c(2)*
i1)*identity-c(2)*azl &
3756 & -(c(1)*
i1+2.0_dp*c(2)*i2-1.5_dp*wv_prime*jznu**(5.0_dp/3.0_dp))/3.0_dp*azu)
3763 i1 = azl(1,1) + azl(2,2) + azl(3,3)
3765 i2 = 0.5_dp * (
i1**2 - azl_squared(1,1) - azl_squared(2,2) - azl_squared(3,3))
3784 piola_tensor=2.0_dp*jznu**(-2.0_dp/3.0_dp)*((c(1)+c(2)*
i1)*identity-c(2)*azl &
3785 & -(c(1)*
i1+2.0_dp*c(2)*i2-1.5_dp*p*jznu**(5.0_dp/3.0_dp))/3.0_dp*azu)
3789 dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
3790 node1=dependent_field%decomposition%domain(1)%ptr%topology%elements%elements(element_number)%element_nodes(13)
3791 node2=dependent_field%decomposition%domain(1)%ptr%topology%elements%elements(element_number)%element_nodes(15)
3793 NULLIFY(field_variable)
3795 CALL field_variable_get(dependent_field,field_v_variable_type,field_variable,err,error,*999)
3796 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
3797 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,left_node(1), &
3799 dof_idx=field_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
3800 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,left_node(2), &
3802 dof_idx=field_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
3803 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,left_node(3), &
3806 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
3807 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,right_node(1), &
3809 dof_idx=field_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
3810 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,right_node(2), &
3812 dof_idx=field_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
3813 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,right_node(3), &
3816 dist1=sqrt((right_node(1)-left_node(1))*(right_node(1)-left_node(1))+ &
3817 & (right_node(2)-left_node(2))*(right_node(2)-left_node(2))+ &
3818 & (right_node(3)-left_node(3))*(right_node(3)-left_node(3)))
3820 NULLIFY(field_variable)
3822 CALL field_variable_get(dependent_field,field_u_variable_type,field_variable,err,error,*999)
3823 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
3824 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,left_node(1), &
3826 dof_idx=field_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
3827 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,left_node(2), &
3829 dof_idx=field_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
3830 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,left_node(3), &
3833 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
3834 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,right_node(1), &
3836 dof_idx=field_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
3837 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,right_node(2), &
3839 dof_idx=field_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
3840 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,right_node(3), &
3843 dist2=sqrt((right_node(1)-left_node(1))*(right_node(1)-left_node(1))+ &
3844 & (right_node(2)-left_node(2))*(right_node(2)-left_node(2))+ &
3845 & (right_node(3)-left_node(3))*(right_node(3)-left_node(3)))
3848 velo=(dist2-dist1)/delta_t
3863 xb_distortion=8.0e-9_dp*(1+velo/v_max)
3865 xb_stiffness=2.2e-3_dp
3867 reference_volume=1.4965e+06_dp
3868 max_xb_number_per_volume=120.0_dp*2.0_dp/reference_volume
3869 energy_per_xb=0.5_dp*xb_stiffness*xb_distortion**2
3871 sarco_length=dzdnu(1,1)
3874 IF(sarco_length.LE.0.635_dp)
THEN 3876 ELSE IF(sarco_length.LE.0.835_dp)
THEN 3877 force_length=4.2_dp*(sarco_length-0.635_dp)
3878 ELSE IF(sarco_length.LE.1.0_dp)
THEN 3879 force_length=0.84_dp+0.9697_dp*(sarco_length-0.835_dp)
3880 ELSE IF(sarco_length.LE.1.125_dp)
THEN 3882 ELSE IF(sarco_length.LE.1.825_dp)
THEN 3883 force_length=1.0_dp-1.4286_dp*(sarco_length-1.125_dp)
3889 xb_energy_per_volume=max_xb_number_per_volume*force_length*c(8)*energy_per_xb*10.0_dp**23
3897 f_a_inv(1,1)=1.0_dp/lambda_a
3912 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,-1,err)
3913 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
3914 lwork=min(lwmax,int(work(1)))
3915 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,lwork,err)
3916 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
3923 ematrix_1(i,j)=evector_1(i)*evector_1(j)
3924 ematrix_2(i,j)=evector_2(i)*evector_2(j)
3925 ematrix_3(i,j)=evector_3(i)*evector_3(j)
3936 free_energy_0=0.0_dp
3938 free_energy_0=free_energy_0+c(i)/c(i+3)*( &
3939 & evalues(1)**(c(i+3)/2.0_dp)+ &
3940 & evalues(2)**(c(i+3)/2.0_dp)+ &
3941 & evalues(3)**(c(i+3)/2.0_dp)-3.0_dp)
3943 free_energy_0=c(7)*free_energy_0
3945 free_energy=free_energy_0
3947 VALUE=xb_energy_per_volume-(free_energy-free_energy_0)
3955 DO WHILE (abs(
VALUE).GE.tol)
3957 IF (abs(
VALUE).GE.tol1)
THEN 3958 lambda_a=up-(up-low)/2.0_dp
3961 f_a_inv(1,1)=1.0_dp/lambda_a
3972 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,-1,err)
3973 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
3974 lwork=min(lwmax,int(work(1)))
3975 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,lwork,err)
3976 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
3983 ematrix_1(i,j)=evector_1(i)*evector_1(j)
3984 ematrix_2(i,j)=evector_2(i)*evector_2(j)
3985 ematrix_3(i,j)=evector_3(i)*evector_3(j)
3998 free_energy=free_energy+c(i)/c(i+3)*( &
3999 & evalues(1)**(c(i+3)/2.0_dp)+ &
4000 & evalues(2)**(c(i+3)/2.0_dp)+ &
4001 & evalues(3)**(c(i+3)/2.0_dp)-3.0_dp)
4003 free_energy=c(7)*free_energy
4005 VALUE=xb_energy_per_volume-(free_energy-free_energy_0)
4007 IF (
VALUE .GE. 0.0_dp)
THEN 4024 & c(i)*evalues(1)**(c(i+3)/2.0_dp-1.0_dp)*temp1+ &
4025 & c(i)*evalues(2)**(c(i+3)/2.0_dp-1.0_dp)*temp2+ &
4026 & c(i)*evalues(3)**(c(i+3)/2.0_dp-1.0_dp)*temp3
4028 slope=temp(1,1)*c(7)
4029 lambda_a=lambda_a-
VALUE/slope
4036 f_a_inv(1,1)=1.0_dp/lambda_a
4047 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,-1,err)
4048 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
4049 lwork=min(lwmax,int(work(1)))
4050 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,lwork,err)
4051 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
4058 ematrix_1(i,j)=evector_1(i)*evector_1(j)
4059 ematrix_2(i,j)=evector_2(i)*evector_2(j)
4060 ematrix_3(i,j)=evector_3(i)*evector_3(j)
4073 free_energy=free_energy+c(i)/c(i+3)*( &
4074 & evalues(1)**(c(i+3)/2.0_dp)+ &
4075 & evalues(2)**(c(i+3)/2.0_dp)+ &
4076 & evalues(3)**(c(i+3)/2.0_dp)-3.0_dp)
4078 free_energy=c(7)*free_energy
4080 VALUE=xb_energy_per_volume-(free_energy-free_energy_0)
4129 piola_tensor=piola_tensor+ &
4130 & c(i)*evalues(1)**(c(i+3)/2.0_dp-1.0_dp)*n1+ &
4131 & c(i)*evalues(2)**(c(i+3)/2.0_dp-1.0_dp)*n2+ &
4132 & c(i)*evalues(3)**(c(i+3)/2.0_dp-1.0_dp)*n3
4134 piola_tensor=piola_tensor*c(7)+2.0_dp*p*azu
4139 dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
4140 node1=dependent_field%decomposition%domain(1)%ptr%topology%elements%elements(element_number)%element_nodes(13)
4141 node2=dependent_field%decomposition%domain(1)%ptr%topology%elements%elements(element_number)%element_nodes(15)
4143 NULLIFY(field_variable)
4145 CALL field_variable_get(dependent_field,field_v_variable_type,field_variable,err,error,*999)
4146 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
4147 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,left_node(1), &
4149 dof_idx=field_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
4150 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,left_node(2), &
4152 dof_idx=field_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
4153 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,left_node(3), &
4156 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
4157 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,right_node(1), &
4159 dof_idx=field_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
4160 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,right_node(2), &
4162 dof_idx=field_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
4163 CALL field_parameter_set_get_local_dof(dependent_field,field_v_variable_type,field_values_set_type,dof_idx,right_node(3), &
4166 dist1=sqrt((right_node(1)-left_node(1))*(right_node(1)-left_node(1))+ &
4167 & (right_node(2)-left_node(2))*(right_node(2)-left_node(2))+ &
4168 & (right_node(3)-left_node(3))*(right_node(3)-left_node(3)))
4170 NULLIFY(field_variable)
4172 CALL field_variable_get(dependent_field,field_u_variable_type,field_variable,err,error,*999)
4173 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
4174 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,left_node(1), &
4176 dof_idx=field_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
4177 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,left_node(2), &
4179 dof_idx=field_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node1)%DERIVATIVES(1)%VERSIONS(1)
4180 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,left_node(3), &
4183 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
4184 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,right_node(1), &
4186 dof_idx=field_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
4187 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,right_node(2), &
4189 dof_idx=field_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP%NODES(node2)%DERIVATIVES(1)%VERSIONS(1)
4190 CALL field_parameter_set_get_local_dof(dependent_field,field_u_variable_type,field_values_set_type,dof_idx,right_node(3), &
4193 dist2=sqrt((right_node(1)-left_node(1))*(right_node(1)-left_node(1))+ &
4194 & (right_node(2)-left_node(2))*(right_node(2)-left_node(2))+ &
4195 & (right_node(3)-left_node(3))*(right_node(3)-left_node(3)))
4198 velo=(dist2-dist1)/delta_t
4204 CALL field_parameter_set_update_gauss_point(dependent_field,field_u1_variable_type,field_values_set_type,gauss_point_number, &
4205 & element_number,2,velo,err,error,*999)
4209 NULLIFY(independent_field)
4210 independent_field=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
4211 NULLIFY(field_variable)
4212 CALL field_variable_get(independent_field,field_u_variable_type,field_variable,err,error,*999)
4214 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4216 CALL field_parameter_set_get_local_dof(independent_field,field_u_variable_type,field_values_set_type,dof_idx,a_1, &
4218 dof_idx=field_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4220 CALL field_parameter_set_get_local_dof(independent_field,field_u_variable_type,field_values_set_type,dof_idx,a_2, &
4222 dof_idx=field_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4224 CALL field_parameter_set_get_local_dof(independent_field,field_u_variable_type,field_values_set_type,dof_idx,x_1, &
4226 dof_idx=field_variable%COMPONENTS(4)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4228 CALL field_parameter_set_get_local_dof(independent_field,field_u_variable_type,field_values_set_type,dof_idx,x_2, &
4232 sarco_length=dzdnu(1,1)
4234 IF(sarco_length.LE.0.635_dp)
THEN 4236 ELSE IF(sarco_length.LE.0.835_dp)
THEN 4237 force_length=4.2_dp*(sarco_length-0.635_dp)
4238 ELSE IF(sarco_length.LE.1.0_dp)
THEN 4239 force_length=0.84_dp+0.9697_dp*(sarco_length-0.835_dp)
4240 ELSE IF(sarco_length.LE.1.125_dp)
THEN 4242 ELSE IF(sarco_length.LE.1.825_dp)
THEN 4243 force_length=1.0_dp-1.4286_dp*(sarco_length-1.125_dp)
4248 reference_volume=1.4965e+06_dp
4249 max_xb_number_per_volume=120.0_dp*2.0_dp/reference_volume
4250 energy_per_xb=0.5_dp*x_2**2*c(8)
4253 xb_energy_per_volume=max_xb_number_per_volume*force_length*energy_per_xb*a_2*10.0_dp**23
4259 f_a_inv(1,1)=1.0_dp/lambda_a
4269 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,-1,err)
4270 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
4271 lwork=min(lwmax,int(work(1)))
4272 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,lwork,err)
4273 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
4280 ematrix_1(i,j)=evector_1(i)*evector_1(j)
4281 ematrix_2(i,j)=evector_2(i)*evector_2(j)
4282 ematrix_3(i,j)=evector_3(i)*evector_3(j)
4293 free_energy_0=0.0_dp
4295 free_energy_0=free_energy_0+c(i)/c(i+3)*( &
4296 & evalues(1)**(c(i+3)/2.0_dp)+ &
4297 & evalues(2)**(c(i+3)/2.0_dp)+ &
4298 & evalues(3)**(c(i+3)/2.0_dp)-3.0_dp)
4300 free_energy_0=c(7)*free_energy_0
4302 free_energy=free_energy_0
4304 VALUE=xb_energy_per_volume-(free_energy-free_energy_0)
4315 DO WHILE (abs(
VALUE).GE.tol)
4318 IF (abs(
VALUE).GE.tol1)
THEN 4319 lambda_a=up-(up-low)/2.0_dp
4322 IF(lambda_a<tol)
THEN 4323 CALL flagwarning(
"lambda_a is close to zero",err,error,*999)
4327 lambda_a=lambda_a+tol
4329 f_a_inv(1,1)=1.0_dp/lambda_a
4337 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,-1,err)
4338 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
4339 lwork=min(lwmax,int(work(1)))
4340 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,lwork,err)
4341 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
4348 ematrix_1(i,j)=evector_1(i)*evector_1(j)
4349 ematrix_2(i,j)=evector_2(i)*evector_2(j)
4350 ematrix_3(i,j)=evector_3(i)*evector_3(j)
4363 free_energy=free_energy+c(i)/c(i+3)*( &
4364 & evalues(1)**(c(i+3)/2.0_dp)+ &
4365 & evalues(2)**(c(i+3)/2.0_dp)+ &
4366 & evalues(3)**(c(i+3)/2.0_dp)-3.0_dp)
4368 free_energy=c(7)*free_energy
4370 VALUE=xb_energy_per_volume-(free_energy-free_energy_0)
4372 IF (
VALUE.GE.0)
THEN 4390 & c(i)*evalues(1)**(c(i+3)/2.0_dp-1.0_dp)*temp1+ &
4391 & c(i)*evalues(2)**(c(i+3)/2.0_dp-1.0_dp)*temp2+ &
4392 & c(i)*evalues(3)**(c(i+3)/2.0_dp-1.0_dp)*temp3
4394 slope=temp(1,1)*c(7)
4395 lambda_a=lambda_a-
VALUE/slope
4402 f_a_inv(1,1)=1.0_dp/lambda_a
4410 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,-1,err)
4411 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
4412 lwork=min(lwmax,int(work(1)))
4413 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,lwork,err)
4414 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
4421 ematrix_1(i,j)=evector_1(i)*evector_1(j)
4422 ematrix_2(i,j)=evector_2(i)*evector_2(j)
4423 ematrix_3(i,j)=evector_3(i)*evector_3(j)
4436 free_energy=free_energy+c(i)/c(i+3)*( &
4437 & evalues(1)**(c(i+3)/2.0_dp)+ &
4438 & evalues(2)**(c(i+3)/2.0_dp)+ &
4439 & evalues(3)**(c(i+3)/2.0_dp)-3.0_dp)
4441 free_energy=c(7)*free_energy
4443 VALUE=xb_energy_per_volume-(free_energy-free_energy_0)
4449 piola_tensor=piola_tensor+ &
4450 & c(i)*evalues(1)**(c(i+3)/2.0_dp-1.0_dp)*n1+ &
4451 & c(i)*evalues(2)**(c(i+3)/2.0_dp-1.0_dp)*n2+ &
4452 & c(i)*evalues(3)**(c(i+3)/2.0_dp-1.0_dp)*n3
4454 piola_tensor=piola_tensor*c(7)+2.0_dp*p*azu
4457 lambda_f=sqrt(azl(1,1))
4458 CALL field_parameter_set_update_gauss_point(dependent_field,field_u1_variable_type,field_values_set_type,gauss_point_number, &
4459 & element_number,1,lambda_f,err,error,*999)
4464 NULLIFY(independent_field)
4465 independent_field=>equations_set%INDEPENDENT%INDEPENDENT_FIELD
4466 NULLIFY(field_variable)
4467 CALL field_variable_get(independent_field,field_u_variable_type,field_variable,err,error,*999)
4469 dof_idx=field_variable%COMPONENTS(5)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4471 CALL field_parameter_set_get_local_dof(independent_field,field_u_variable_type,field_values_set_type,dof_idx,lambda_a, &
4475 f_a_inv(1,1)=1.0_dp/lambda_a
4485 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,-1,err)
4486 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
4487 lwork=min(lwmax,int(work(1)))
4488 CALL dsyev(
'V',
'U',3,c_e,3,evalues,work,lwork,err)
4489 IF(err.NE.0)
CALL flagerror(
"Error in Eigenvalue computation",err,error,*999)
4496 ematrix_1(i,j)=evector_1(i)*evector_1(j)
4497 ematrix_2(i,j)=evector_2(i)*evector_2(j)
4498 ematrix_3(i,j)=evector_3(i)*evector_3(j)
4511 piola_tensor=piola_tensor+ &
4512 & c(i)*evalues(1)**(c(i+3)/2.0_dp-1.0_dp)*n1+ &
4513 & c(i)*evalues(2)**(c(i+3)/2.0_dp-1.0_dp)*n2+ &
4514 & c(i)*evalues(3)**(c(i+3)/2.0_dp-1.0_dp)*n3
4516 piola_tensor=piola_tensor*c(7)+2.0_dp*p*azu
4531 piola_tensor(1,3)=2.0_dp*(c(2)*(-azl(3,1)))+p*azu(1,3)
4532 piola_tensor(2,3)=2.0_dp*(c(2)*(-azl(3,2)))+p*azu(2,3)
4533 piola_tensor(3,1)=piola_tensor(1,3)
4534 piola_tensor(3,2)=piola_tensor(2,3)
4535 piola_tensor(3,3)=2.0_dp*(c(1)+c(2)*(azl(1,1)+azl(2,2)))+p*azu(3,3)
4539 azl(3,3) = 1.0_dp / ((azl(1,1) * azl(2,2)) - (azl(1,2) * azl(2,1)))
4541 p = -1.0_dp*((c(1) + c(2) * (azl(1,1) + azl(2,2))) * azl(3,3))
4543 piola_tensor(:,3) = 0.0_dp
4544 piola_tensor(3,:) = 0.0_dp
4546 piola_tensor(1,1)=2.0_dp*(c(1)+c(2)*(azl(2,2)+azl(3,3)))+p*azu(1,1)
4547 piola_tensor(1,2)=2.0_dp*( c(2)*(-azl(2,1)))+p*azu(1,2)
4548 piola_tensor(2,1)=piola_tensor(1,2)
4549 piola_tensor(2,2)=2.0_dp*(c(1)+c(2)*(azl(3,3)+azl(1,1)))+p*azu(2,2)
4552 SELECT CASE(equations_set_subtype)
4558 CALL field_parametersetgetlocalgausspoint(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
4559 & field_u_variable_type,field_values_set_type,gauss_point_number,element_number,1,active_stress_11, &
4562 CALL field_parametersetgetlocalgausspoint(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
4563 & field_u_variable_type,field_values_set_type,gauss_point_number,element_number,2,active_stress_22, &
4566 CALL field_parametersetgetlocalgausspoint(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
4567 & field_u_variable_type,field_values_set_type,gauss_point_number,element_number,3,active_stress_33, &
4570 piola_tensor(1,1)=piola_tensor(1,1)+active_stress_11
4571 piola_tensor(2,2)=piola_tensor(2,2)+active_stress_22
4572 piola_tensor(3,3)=piola_tensor(3,3)+active_stress_33
4576 piola_tensor(1,1)=piola_tensor(1,1)+independent_interpolated_point%VALUES(1,
no_part_deriv)
4580 IF(azl(1,1) > 1.0_dp)
THEN 4581 piola_tensor(1,1)=piola_tensor(1,1)+c(3)/azl(1,1)*(azl(1,1)**(c(4)/2.0_dp)-1.0_dp)
4584 CALL field_parametersetgetlocalgausspoint(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
4585 & field_u_variable_type,field_values_set_type,gauss_point_number,element_number,1,
VALUE, &
4588 VALUE=
VALUE/sqrt(azl(1,1))*c(5)
4597 piola_tensor(1,1)=piola_tensor(1,1)+
VALUE 4601 IF(azl(1,1) > 1.0_dp)
THEN 4602 piola_tensor(1,1)=piola_tensor(1,1)+c(3)/azl(1,1)*(azl(1,1)**(c(4)/2.0_dp)-1.0_dp)
4605 CALL field_variable_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,field_variable,err,error,*999)
4606 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4608 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
4609 & field_values_set_type,dof_idx,
VALUE,err,error,*999)
4611 IF(
VALUE.LT.0.0_dp)
VALUE=0.0_dp
4614 VALUE=
VALUE/sqrt(azl(1,1))*c(5)
4616 piola_tensor(1,1)=piola_tensor(1,1)+
VALUE 4619 dof_idx=field_variable%COMPONENTS(2)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4621 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
4622 & field_values_set_type,dof_idx,titin_unbound,err,error,*999)
4624 dof_idx=field_variable%COMPONENTS(3)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4626 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
4627 & field_values_set_type,dof_idx,titin_bound,err,error,*999)
4629 dof_idx=field_variable%COMPONENTS(6)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4631 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
4632 & field_values_set_type,dof_idx,activation,err,error,*999)
4634 IF(activation.GT.1.0_dp) activation=1.0_dp
4635 IF(activation.LT.0.0_dp) activation=0.0_dp
4638 activation=c(6)*activation
4641 titin_value=activation*titin_bound+(1.0_dp-activation)*titin_unbound
4644 titin_value=titin_value/sqrt(azl(1,1))*c(5)
4646 piola_tensor(1,1)=piola_tensor(1,1)+titin_value
4649 dof_idx=field_variable%COMPONENTS(4)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4651 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
4652 & field_values_set_type,dof_idx,titin_unbound_cross_fibre,err,error,*999)
4654 dof_idx=field_variable%COMPONENTS(5)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4656 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
4657 & field_values_set_type,dof_idx,titin_bound_cross_fibre,err,error,*999)
4660 titin_value_cross_fibre=activation*titin_bound_cross_fibre+(1.0_dp-activation)*titin_unbound_cross_fibre
4662 titin_value_cross_fibre=titin_value_cross_fibre*c(5)
4664 piola_tensor(2,2)=piola_tensor(2,2)+titin_value_cross_fibre
4665 piola_tensor(3,3)=piola_tensor(3,3)+titin_value_cross_fibre
4669 IF(azl(1,1) > 1.0_dp)
THEN 4672 piola_tensor(1,1)=piola_tensor(1,1)+0.355439810963035_dp/azl(1,1)*(azl(1,1)**(12.660539325481963_dp/2.0_dp)-1.0_dp)
4675 IF(azl(2,2) > 1.0_dp)
THEN 4676 piola_tensor(2,2)=piola_tensor(2,2)+5316.372204148964_dp/azl(2,2)*(azl(2,2)**(0.014991843974911_dp/2.0_dp)-1.0_dp)
4678 IF(azl(3,3) > 1.0_dp)
THEN 4679 piola_tensor(3,3)=piola_tensor(3,3)+5316.372204148964_dp/azl(3,3)*(azl(3,3)**(0.014991843974911_dp/2.0_dp)-1.0_dp)
4683 CALL field_variable_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,field_variable,err,error,*999)
4684 dof_idx=field_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP%GAUSS_POINTS(gauss_point_number, &
4686 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
4687 & field_values_set_type,dof_idx,
VALUE,err,error,*999)
4692 VALUE=
VALUE/sqrt(azl(1,1))*c(5)
4699 val2=100.0_dp*(sqrt(azl(1,1))-1)
4700 VALUE=
VALUE+val2/sqrt(azl(1,1))
4701 piola_tensor(1,1)=piola_tensor(1,1)+
VALUE 4715 IF(azl(1,1) > 1.0_dp)
THEN 4716 piola_tensor(1,1)=piola_tensor(1,1)+c(3)/azl(1,1)*(azl(1,1)**(c(4)/2.0_dp)-1.0_dp)
4724 IF(azl(1,1) > 1.0_dp)
THEN 4725 piola_tensor(1,1)=piola_tensor(1,1)+c(3)/azl(1,1)*(azl(1,1)**(c(4)/2.0_dp)-1.0_dp)
4728 if((sqrt(azl(1,1))>0.72_dp).AND.(sqrt(azl(1,1))<1.68_dp))
then 4730 VALUE=(-25.0_dp/4.0_dp*azl(1,1)/1.2_dp/1.2_dp + 25.0_dp/2.0_dp*sqrt(azl(1,1))/1.2_dp - 5.25_dp)
4731 VALUE=
VALUE*(1.0_dp/sqrt(azl(1,1)))*20.0_dp*c(5)
4732 piola_tensor(1,1)=piola_tensor(1,1)+
VALUE 4747 IF(azl(1,1) > 1.0_dp)
THEN 4748 piola_tensor(1,1)=piola_tensor(1,1)+c(3)/azl(1,1)*(azl(1,1)**(c(4)/2.0_dp)-1.0_dp)
4750 IF(azl(2,2) > 1.0_dp)
THEN 4751 piola_tensor(2,2)=piola_tensor(2,2)+c(5)/azl(2,2)*(azl(2,2)**(c(6)/2.0_dp)-1.0_dp)
4753 IF(azl(3,3) > 1.0_dp)
THEN 4754 piola_tensor(3,3)=piola_tensor(3,3)+c(7)/azl(3,3)*(azl(3,3)**(c(8)/2.0_dp)-1.0_dp)
4773 IF(azl(1,1) > 1.0_dp)
THEN 4774 piola_tensor(1,1)=piola_tensor(1,1)+c(3)/azl(1,1)*(azl(1,1)**(c(4)/2.0_dp)-1.0_dp)
4776 IF(azl(2,2) > 1.0_dp)
THEN 4777 piola_tensor(2,2)=piola_tensor(2,2)+c(5)/azl(2,2)*(azl(2,2)**(c(6)/2.0_dp)-1.0_dp)
4779 IF(azl(3,3) > 1.0_dp)
THEN 4780 piola_tensor(3,3)=piola_tensor(3,3)+c(7)/azl(3,3)*(azl(3,3)**(c(8)/2.0_dp)-1.0_dp)
4783 val1=sqrt(azl(1,1))/c(9)
4784 IF((val1>0.7_dp).AND.(val1<1.3_dp))
THEN 4786 VALUE=(-11.1111_dp*val1*val1+22.2222_dp*val1-10.1111_dp)
4788 VALUE=
VALUE*c(10)*c(11)/sqrt(azl(1,1))
4794 val2=c(11)*c(12)*(sqrt(azl(1,1))-1)
4795 VALUE=
VALUE+val2/sqrt(azl(1,1))
4796 piola_tensor(1,1)=piola_tensor(1,1)+
VALUE 4824 val1=c(1)*c(10)+c(5)*(1.0_dp-c(10))
4825 val2=c(2)*c(10)+c(6)*(1.0_dp-c(10))
4828 piola_tensor(1,1)=2.0_dp*(val1+val2*(azl(2,2)+azl(3,3))+p*azu(1,1))
4829 piola_tensor(1,2)=2.0_dp*( val2*(-azl(2,1)) +p*azu(1,2))
4830 piola_tensor(1,3)=2.0_dp*( val2*(-azl(3,1)) +p*azu(1,3))
4831 piola_tensor(2,1)=piola_tensor(1,2)
4832 piola_tensor(2,2)=2.0_dp*(val1+val2*(azl(3,3)+azl(1,1))+p*azu(2,2))
4833 piola_tensor(2,3)=2.0_dp*( val2*(-azl(3,2)) +p*azu(2,3))
4834 piola_tensor(3,1)=piola_tensor(1,3)
4835 piola_tensor(3,2)=piola_tensor(2,3)
4836 piola_tensor(3,3)=2.0_dp*(val1+val2*(azl(1,1)+azl(2,2))+p*azu(3,3))
4839 IF(azl(1,1) > 1.0_dp)
THEN 4840 val1=c(3)/azl(1,1)*(azl(1,1)**(c(4)/2.0_dp)-1.0_dp)
4841 val2=c(7)/azl(1,1)*(azl(1,1)**(c(8)/2.0_dp)-1.0_dp)
4842 piola_tensor(1,1)=piola_tensor(1,1)+(val1*c(10)+val2*(1.0_dp-c(10)))
4846 IF((sqrt(azl(1,1))>0.84_dp).AND.(sqrt(azl(1,1))<1.96_dp))
THEN 4847 VALUE=(-25.0_dp/4.0_dp*azl(1,1)/1.4_dp/1.4_dp + 25.0_dp/2.0_dp*sqrt(azl(1,1))/1.4_dp - 5.25_dp)
4848 VALUE=
VALUE*(1.0_dp/sqrt(azl(1,1)))*c(9)*c(10)*c(11)
4849 piola_tensor(1,1)=piola_tensor(1,1)+
VALUE 4856 piola_tensor=c(1)*c(2)*exp(c(2)*(azl(1,1)+azl(2,2)+azl(3,3)-3.0_dp))*identity+2.0_dp*p*azu
4865 p=darcy_dependent_interpolated_point%VALUES(1,
no_part_deriv)
4867 i1=azl(1,1)+azl(2,2)+azl(3,3)
4868 temp=matmul(azl,azl)
4869 i2=0.5_dp*(
i1**2.0_dp-temp(1,1)-temp(2,2)-temp(3,3))
4873 piola_tensor=2.0_dp*c(1)*jznu**(-2.0_dp/3.0_dp)*(identity-(1.0_dp/3.0_dp)*
i1*azu)
4874 piola_tensor=piola_tensor+2.0_dp*c(2)*jznu**(-4.0_dp/3.0_dp)*(
i1*identity-azlt-(2.0_dp/3.0_dp)*i2*azu)
4875 piola_tensor=piola_tensor+(c(3)-c(4)*c(5)**2)*(jznu-1.0_dp)*azu
4876 piola_tensor=piola_tensor-c(5)*(p-c(6))*jznu*azu
4877 piola_tensor=piola_tensor+0.5_dp*((p-c(6))**2/c(4))*(dfdjfact/(ffact**2))*jznu*azu
4899 i1=azl(1,1)+azl(2,2)+azl(3,3)
4900 temp=matmul(azl,azl)
4901 i2=0.5_dp*(
i1**2.0_dp-temp(1,1)-temp(2,2)-temp(3,3))
4904 tempterm=2.0_dp*c(4)*c(1)*exp(c(2)*(
i1 - 3.0_dp) + c(3)*(i2 - 3.0_dp)) / (i3**(c(2)+2.0_dp*c(3)))
4905 piola_tensor=c(2)*tempterm*identity + c(3)*tempterm*(
i1*identity-azlt) - (c(2)+2.0_dp*c(3))*tempterm*azut
4906 piola_tensor=piola_tensor - darcy_dependent_interpolated_point%VALUES(1,
no_part_deriv)*jznu*azu
4931 i1=azl(1,1)+azl(2,2)+azl(3,3)
4932 temp=matmul(azl,azl)
4933 i2=0.5_dp*(
i1**2.0_dp-temp(1,1)-temp(2,2)-temp(3,3))
4936 tempterm=2.0_dp*c(4)*c(1)*exp(c(2)*(
i1 - 3.0_dp) + c(3)*(i2 - 3.0_dp)) / (i3**(c(2)+2.0_dp*c(3)))
4937 piola_tensor=c(2)*tempterm*identity + c(3)*tempterm*(
i1*identity-azlt) - (c(2)+2.0_dp*c(3))*tempterm*azut
4938 piola_tensor=piola_tensor - darcy_dependent_interpolated_point%VALUES(1,
no_part_deriv)*jznu*azu
4940 IF((sqrt(azl(1,1))>0.72_dp).AND.(sqrt(azl(1,1))<1.68_dp))
THEN 4941 VALUE=(-25.0_dp/4.0_dp*azl(1,1)/1.2_dp/1.2_dp + 25.0_dp/2.0_dp*sqrt(azl(1,1))/1.2_dp - 5.25_dp)
4946 piola_tensor(1,1) = piola_tensor(1,1) + 1.0_dp/sqrt(azl(1,1))*c(5)*c(6)*
VALUE 4952 piola_tensor(1,3)=(2.0_dp*c(2)*e(1,3))+(2.0_dp*p*azu(1,3))
4953 piola_tensor(2,3)=(2.0_dp*c(2)*e(2,3))+(2.0_dp*p*azu(2,3))
4954 piola_tensor(3,1)=piola_tensor(1,3)
4955 piola_tensor(3,2)=piola_tensor(2,3)
4956 piola_tensor(3,3)=c(1)*(e(1,1)+e(2,2)+e(3,3))+(2.0_dp*e(3,3)*c(2)+(2.0_dp*p*azu(3,3)))
4958 piola_tensor(1,1)=c(1)*(e(1,1)+e(2,2)+e(3,3))+(2.0_dp*e(1,1)*c(2)+(2.0_dp*p*azu(1,1)))
4959 piola_tensor(1,2)=(2.0_dp*c(2)*e(1,2))+(2.0_dp*p*azu(1,2))
4960 piola_tensor(2,1)=piola_tensor(1,2)
4961 piola_tensor(2,2)=c(1)*(e(1,1)+e(2,2)+e(3,3))+(2.0_dp*e(2,2)*c(2)+(2.0_dp*p*azu(2,2)))
4963 CALL field_parametersetgetlocalgausspoint(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
4964 & field_u_variable_type,field_values_set_type,gauss_point_number,element_number,1,active_stress_11, &
4967 CALL field_parametersetgetlocalgausspoint(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
4968 & field_u_variable_type,field_values_set_type,gauss_point_number,element_number,2,active_stress_22, &
4971 CALL field_parametersetgetlocalgausspoint(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
4972 & field_u_variable_type, field_values_set_type,gauss_point_number,element_number,3,active_stress_33, &
4975 piola_tensor(1,1)=piola_tensor(1,1)+active_stress_11
4976 piola_tensor(2,2)=piola_tensor(2,2)+active_stress_22
4977 piola_tensor(3,3)=piola_tensor(3,3)+active_stress_33
4985 tempterm=c(1)*exp(2.0*c(2)*(e(1,1)+e(2,2)+e(3,3))+c(3)*e(1,1)**2+c(4)*(e(2,2)**2+e(3,3)**2+2.0_dp*e(2,3)**2)+ &
4986 & c(5)*2.0_dp*(e(1,2)**2+e(1,3)**2))
4987 piola_tensor(1,1)=(c(2)+c(3)*e(1,1))*tempterm+2.0_dp*p*azu(1,1)
4988 piola_tensor(1,2)=c(5)*e(1,2)*tempterm+2.0_dp*p*azu(1,2)
4989 piola_tensor(1,3)=c(5)*e(1,3)*tempterm+2.0_dp*p*azu(1,3)
4990 piola_tensor(2,1)=piola_tensor(1,2)
4991 piola_tensor(2,2)=(c(2)+c(4)*e(2,2))*tempterm+2.0_dp*p*azu(2,2)
4992 piola_tensor(2,3)=c(4)*e(2,3)*tempterm+2.0_dp*p*azu(2,3)
4993 piola_tensor(3,1)=piola_tensor(1,3)
4994 piola_tensor(3,2)=piola_tensor(2,3)
4995 piola_tensor(3,3)=(c(2)+c(4)*e(3,3))*tempterm+2.0_dp*p*azu(3,3)
5000 q=c(2)*e(1,1)**2 + c(3)*(e(2,2)**2+e(3,3)**2+2.0_dp*e(2,3)**2) + 2.0_dp*c(4)*(e(1,2)**2+e(1,3)**2)
5001 tempterm=0.5_dp*c(1)*exp(q)
5002 piola_tensor(1,1) = 2.0_dp*c(2) * e(1,1)
5003 piola_tensor(2,2) = 2.0_dp*c(3) * e(2,2)
5004 piola_tensor(3,3) = 2.0_dp*c(3) * e(3,3)
5005 piola_tensor(1,2) = 2.0_dp*c(4) * e(1,2)
5006 piola_tensor(2,1) = piola_tensor(1,2)
5007 piola_tensor(1,3) = 2.0_dp*c(4) * e(1,3)
5008 piola_tensor(3,1) = piola_tensor(1,3)
5009 piola_tensor(3,2) = 2.0_dp*c(3) * e(2,3)
5010 piola_tensor(2,3) = piola_tensor(3,2)
5011 piola_tensor = piola_tensor * tempterm
5032 piola_tensor = piola_tensor + p*azu
5037 CALL field_variable_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,field_variable,err,error,*999)
5038 DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
5039 dof_idx=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP% &
5040 & gauss_points(gauss_point_number,element_number)
5041 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
5042 & field_values_set_type,dof_idx,
VALUE,err,error,*999)
5043 piola_tensor(component_idx,component_idx)=piola_tensor(component_idx,component_idx)+
VALUE 5049 i1=azl(1,1)+azl(2,2)+azl(3,3)
5050 piola_tensor(1,1)=c(1)*c(2)*exp(c(2)*(
i1-3))+ &
5051 & c(3)*2.0_dp*(sqrt(azl(1,1))-1)*c(4)*exp(c(4)*(sqrt(azl(1,1))-1)**2)/(2*sqrt(azl(1,1)))+p*azu(1,1)
5052 piola_tensor(2,2)=c(1)*c(2)*exp(c(2)*(
i1-3))+p*azu(2,2)
5053 piola_tensor(3,3)=c(1)*c(2)*exp(c(2)*(
i1-3))+p*azu(3,3)
5054 piola_tensor(1,2)=p*azu(1,2)
5055 piola_tensor(1,3)=p*azu(1,3)
5056 piola_tensor(2,3)=p*azu(2,3)
5057 piola_tensor(2,1)=piola_tensor(1,2)
5058 piola_tensor(3,1)=piola_tensor(1,3)
5059 piola_tensor(3,2)=piola_tensor(2,3)
5060 piola_tensor=piola_tensor*2.0_dp
5066 a = materials_interpolated_point%VALUES(1,1)
5067 b(1,1) = materials_interpolated_point%VALUES(1+1,1)
5068 b(1,2) = materials_interpolated_point%VALUES(1+2,1)
5069 b(1,3) = materials_interpolated_point%VALUES(1+3,1)
5071 b(2,2) = materials_interpolated_point%VALUES(1+4,1)
5072 b(2,3) = materials_interpolated_point%VALUES(1+5,1)
5075 b(3,3) = materials_interpolated_point%VALUES(1+6,1)
5080 e(i,j) = 0.5_dp * (azl(i,j)-1);
5082 e(i,j) = 0.5_dp * azl(i,j);
5084 q = q + b(i,j) * e(i,j) * e(i,j)
5090 piola_tensor(i,j)=a*b(i,j)*e(i,j)*q + p*azu(i,j);
5096 & equations_set%EQUATIONS%INTERPOLATION%MATERIALS_FIELD, piola_tensor(1,1),e(1,1), &
5097 & element_number,gauss_point_number,err,error,*999)
5104 c(1)=materials_interpolated_point%VALUES(1,1)
5105 c(2)=materials_interpolated_point%VALUES(2,1)
5107 piola_tensor(1,1)=c(1)+c(2)*(azl(2,2)+azl(3,3))
5108 piola_tensor(1,2)=c(2)*(-azl(2,1))
5109 piola_tensor(1,3)=c(2)*(-azl(3,1))
5110 piola_tensor(2,1)=piola_tensor(1,2)
5111 piola_tensor(2,2)=c(1)+c(2)*(azl(3,3)+azl(1,1))
5112 piola_tensor(2,3)=c(2)*(-azl(3,2))
5113 piola_tensor(3,1)=piola_tensor(1,3)
5114 piola_tensor(3,2)=piola_tensor(2,3)
5115 piola_tensor(3,3)=c(1)+c(2)*(azl(1,1)+azl(2,2))
5116 piola_tensor=piola_tensor*2.0_dp
5122 & 3,3,azl,
write_string_matrix_name_and_indices,
'(" AZL',
'(",I1,",:)',
' :",3(X,E13.6))', &
5123 &
'(17X,3(X,E13.6))',err,error,*999)
5128 CALL field_parametersetgetlocalgausspoint(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
5129 & field_u_variable_type,field_values_set_type,gauss_point_number,element_number,1,active_stress_11, &
5132 CALL field_parametersetgetlocalgausspoint(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
5133 & field_u_variable_type,field_values_set_type,gauss_point_number,element_number,2,active_stress_22, &
5136 CALL field_parametersetgetlocalgausspoint(equations_set%INDEPENDENT%INDEPENDENT_FIELD, &
5137 & field_u_variable_type,field_values_set_type,gauss_point_number,element_number,3,active_stress_33, &
5140 piola_tensor(1,1)=piola_tensor(1,1)+active_stress_11
5141 piola_tensor(2,2)=piola_tensor(2,2)+active_stress_22
5142 piola_tensor(3,3)=piola_tensor(3,3)+active_stress_33
5146 c(3)=materials_interpolated_point%VALUES(3,1)
5147 piola_tensor=piola_tensor+2.0_dp*c(3)*(i3-sqrt(i3))*azu
5151 SELECT CASE (equations_set_subtype)
5153 c(3)=materials_interpolated_point%VALUES(3,1)
5159 piola_tensor=piola_tensor+c(3)*(sqrt(i3)-1.0_dp)*azu
5160 darcy_mass_increase_entry = 5
5176 c(1)=materials_interpolated_point%VALUES(1,1)
5177 c(2)=materials_interpolated_point%VALUES(2,1)
5178 c(3)=materials_interpolated_point%VALUES(3,1)
5181 tempterm=jznu**(-2.0_dp/3.0_dp)
5185 i1=azl(1,1)+azl(2,2)+azl(3,3)
5186 piola_tensor=c(1)* (temp-1.0_dp/3.0_dp*
i1*tempterm*azu)
5189 temp=matmul(azl,azl)
5190 i2=0.5_dp*(
i1**2.0_dp-(temp(1,1)+temp(2,2)+temp(3,3)))
5191 tempterm=jznu**(-4.0_dp/3.0_dp)
5193 temp(1,1)=azl(2,2)+azl(3,3)
5195 temp(1,2)=-1.0_dp*azl(1,2)
5197 temp(1,3)=-1.0_dp*azl(1,3)
5199 temp(2,2)=azl(1,1)+azl(3,3)
5201 temp(2,3)=-1.0_dp*azl(2,3)
5204 temp(3,3)=azl(1,1)+azl(2,2)
5205 piola_tensor=piola_tensor+c(2)* (tempterm*temp-2.0_dp/3.0_dp*i2*tempterm*azu)
5208 piola_tensor=piola_tensor+(2.0_dp*c(3)*(jznu-1.0_dp)+p)*jznu*azu
5211 piola_tensor=2.0_dp*piola_tensor
5214 darcy_mass_increase_entry = 4
5245 c(1)=materials_interpolated_point%VALUES(1,1)
5246 c(2)=materials_interpolated_point%VALUES(2,1)
5247 c(3)=materials_interpolated_point%VALUES(3,1)
5248 c(4)=materials_interpolated_point%VALUES(4,1)
5249 c(5)=materials_interpolated_point%VALUES(5,1)
5250 c(6)=materials_interpolated_point%VALUES(6,1)
5251 c(7)=materials_interpolated_point%VALUES(7,1)
5252 c(8)=materials_interpolated_point%VALUES(8,1)
5253 i1=azl(1,1)+azl(2,2)+azl(3,3)
5254 tempterm=c(1)*exp(c(2)*(
i1-3.0_dp))
5255 piola_tensor(1,1)=-p*azu(1,1)+tempterm
5256 IF(azl(1,1)>1.0_dp)
THEN 5257 piola_tensor(1,1)=piola_tensor(1,1)+2.0_dp*c(3)*(azl(1,1)-1.0_dp)*exp(c(5)*(azl(1,1)-1.0_dp)**2.0_dp)
5259 piola_tensor(1,2)=-p*azu(1,2)+c(7)*azl(1,2)*exp(c(8)*azl(1,2)**2.0_dp)
5260 piola_tensor(1,3)=-p*azu(1,3)
5261 piola_tensor(2,1)=piola_tensor(1,2)
5262 piola_tensor(2,2)=-p*azu(2,2)+tempterm
5263 IF(azl(2,2)>1.0_dp)
THEN 5264 piola_tensor(2,2)=piola_tensor(2,2)+2.0_dp*c(4)*(azl(2,2)-1.0_dp)*exp(c(6)*(azl(2,2)-1.0_dp)**2.0_dp)
5266 piola_tensor(2,3)=-p*azu(2,3)
5267 piola_tensor(3,1)=piola_tensor(1,3)
5268 piola_tensor(3,2)=piola_tensor(2,3)
5269 piola_tensor(3,3)=-p*azu(3,3)+tempterm
5275 CALL field_variable_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,field_variable,err,error,*999)
5276 DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
5277 dof_idx=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP% &
5278 & gauss_points(gauss_point_number,element_number)
5279 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
5280 & field_values_set_type,dof_idx,
VALUE,err,error,*999)
5281 piola_tensor(component_idx,component_idx)=piola_tensor(component_idx,component_idx)+
VALUE 5286 local_error=
"The third equations set specification of "//
trim(
number_to_vstring(equations_set_subtype,
"*",err,error))// &
5287 &
" is not valid for a finite elasticity type of an elasticity equation set." 5288 CALL flagerror(local_error,err,error,*999)
5294 cauchy_tensor=cauchy_tensor/jznu
5299 & 3,3,piola_tensor,
write_string_matrix_name_and_indices,
'(" PIOLA_TENSOR',
'(",I1,",:)',
' :",3(X,E13.6))', &
5300 &
'(17X,3(X,E13.6))',err,error,*999)
5302 & 3,3,cauchy_tensor,
write_string_matrix_name_and_indices,
'(" CAUCHY_TENSOR',
'(",I1,",:)',
' :",3(X,E13.6))', &
5303 &
'(17X,3(X,E13.6))',err,error,*999)
5307 exits(
"FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR")
5309 999 errorsexits(
"FINITE_ELASTICITY_GAUSS_CAUCHY_TENSOR",err,error)
5319 & dependentfield, deformationgradienttensor,growthtensor,elasticdeformationgradienttensor,jg,je,err,error,*)
5323 INTEGER(INTG),
INTENT(IN) :: numberOfDimensions
5324 INTEGER(INTG),
INTENT(IN) :: gaussPointNumber
5325 INTEGER(INTG),
INTENT(IN) :: elementNumber
5327 REAL(DP),
INTENT(IN) :: deformationGradientTensor(3,3)
5328 REAL(DP),
INTENT(OUT) :: growthTensor(3,3)
5329 REAL(DP),
INTENT(OUT) :: elasticDeformationGradientTensor(3,3)
5330 REAL(DP),
INTENT(OUT) :: Jg
5331 REAL(DP),
INTENT(OUT) :: Je
5332 INTEGER(INTG),
INTENT(OUT) :: err
5335 REAL(DP) :: growthTensorInverse(3,3), growthTensorInverseTranspose(3,3)
5337 enters(
"FiniteElasticity_GaussGrowthTensor",err,error,*999)
5339 IF(
ASSOCIATED(equationsset))
THEN 5342 CALL field_parametersetgetlocalgausspoint(dependentfield,field_u3_variable_type,field_values_set_type, &
5343 & gausspointnumber,elementnumber,1,growthtensor(1,1),err,error,*999)
5344 IF(numberofdimensions>1)
THEN 5345 CALL field_parametersetgetlocalgausspoint(dependentfield,field_u3_variable_type,field_values_set_type, &
5346 & gausspointnumber,elementnumber,2,growthtensor(2,2),err,error,*999)
5347 IF(numberofdimensions>2)
THEN 5348 CALL field_parametersetgetlocalgausspoint(dependentfield,field_u3_variable_type,field_values_set_type, &
5349 & gausspointnumber,elementnumber,3,growthtensor(3,3),err,error,*999)
5353 CALL invert(growthtensor,growthtensorinverse,jg,err,error,*999)
5355 CALL matrixproduct(deformationgradienttensor,growthtensorinverse,elasticdeformationgradienttensor,err,error,*999)
5358 elasticdeformationgradienttensor=deformationgradienttensor
5360 je=
determinant(elasticdeformationgradienttensor,err,error)
5363 CALL flagerror(
"Equations set is not associated.",err,error,*999)
5370 CALL writestringmatrix(
diagnostic_output_type,1,1,3,1,1,3,3,3,deformationgradienttensor, &
5371 &
write_string_matrix_name_and_indices,
'(" F',
'(",I1,",:)',
' :",3(X,E13.6))',
'(13X,3(X,E13.6))',err,error,*999)
5375 CALL writestringmatrix(
diagnostic_output_type,1,1,3,1,1,3,3,3,elasticdeformationgradienttensor, &
5376 &
write_string_matrix_name_and_indices,
'(" Fe',
'(",I1,",:)',
' :",3(X,E13.6))',
'(13X,3(X,E13.6))',err,error,*999)
5379 CALL writestringmatrix(
diagnostic_output_type,1,1,3,1,1,3,3,3,growthtensor, &
5380 &
write_string_matrix_name_and_indices,
'(" Fg',
'(",I1,",:)',
' :",3(X,E13.6))',
'(13X,3(X,E13.6))',err,error,*999)
5384 exits(
"FiniteElasticity_GaussGrowthTensor")
5386 999 errorsexits(
"FiniteElasticity_GaussGrowthTensor",err,error)
5397 & fingerdeformationtensor, jacobian,greenstraintensor,err,error,*)
5400 REAL(DP),
INTENT(IN) :: deformationGradientTensor(3,3)
5401 REAL(DP),
INTENT(OUT) :: rightCauchyDeformationTensor(3,3)
5402 REAL(DP),
INTENT(OUT) :: fingerDeformationTensor(3,3)
5403 REAL(DP),
INTENT(OUT) :: Jacobian
5404 REAL(DP),
INTENT(OUT) :: greenStrainTensor(3,3)
5405 INTEGER(INTG),
INTENT(OUT) :: err
5411 enters(
"FiniteElasticity_StrainTensor",err,error,*999)
5413 CALL matrixtransposeproduct(deformationgradienttensor,deformationgradienttensor,rightcauchydeformationtensor,err,error,*999)
5414 CALL invert(rightcauchydeformationtensor,fingerdeformationtensor,i3,err,error,*999)
5415 jacobian=
determinant(deformationgradienttensor,err,error)
5417 greenstraintensor=0.5_dp*rightcauchydeformationtensor
5419 greenstraintensor(i,i)=greenstraintensor(i,i)-0.5_dp
5428 &
' :",3(X,E13.6))',
'(12X,3(X,E13.6))',err,error,*999)
5432 &
' :",3(X,E13.6))',
'(12X,3(X,E13.6))',err,error,*999)
5437 &
' :",3(X,E13.6))',
'(12X,3(X,E13.6))',err,error,*999)
5440 exits(
"FiniteElasticity_StrainTensor")
5442 999 errorsexits(
"FiniteElasticity_StrainTensor",err,error)
5453 & materials_interpolated_point,stress_tensor,dzdnu,jznu,
element_number,gauss_point_number,err,error,*)
5458 REAL(DP),
INTENT(OUT) :: STRESS_TENSOR(:)
5459 REAL(DP),
INTENT(IN) :: DZDNU(3,3)
5460 REAL(DP),
INTENT(IN) :: Jznu
5461 INTEGER(INTG),
INTENT(IN) :: ELEMENT_NUMBER,GAUSS_POINT_NUMBER
5462 INTEGER(INTG),
INTENT(OUT) :: ERR
5465 INTEGER(INTG) :: PRESSURE_COMPONENT,component_idx,dof_idx
5468 REAL(DP) :: TEMPTERM1,TEMPTERM2,
VALUE 5469 REAL(DP) :: ONETHIRD_TRACE
5472 REAL(DP) :: MOD_DZDNU(3,3),MOD_DZDNUT(3,3),AZL(3,3)
5473 REAL(DP) :: B(6),E(6),DQ_DE(6)
5474 REAL(DP),
POINTER :: C(:)
5476 enters(
"FINITE_ELASTICITY_GAUSS_STRESS_TENSOR",err,error,*999)
5478 NULLIFY(field_variable,c)
5483 mod_dzdnu=dzdnu*jznu**(-1.0_dp/3.0_dp)
5488 SELECT CASE(equations_set%specification(3))
5490 pressure_component=dependent_interpolated_point%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
5491 p=dependent_interpolated_point%VALUES(pressure_component,
no_part_deriv)
5496 i1=azl(1,1)+azl(2,2)+azl(3,3)
5497 tempterm1=-2.0_dp*c(2)
5498 tempterm2=2.0_dp*(c(1)+
i1*c(2))
5499 stress_tensor(1)=tempterm1*azl(1,1)+tempterm2
5500 stress_tensor(2)=tempterm1*azl(2,2)+tempterm2
5501 stress_tensor(3)=tempterm1*azl(3,3)+tempterm2
5502 stress_tensor(4)=tempterm1*azl(2,1)
5503 stress_tensor(5)=tempterm1*azl(3,1)
5504 stress_tensor(6)=tempterm1*azl(3,2)
5510 CALL field_variable_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,field_variable,err,error,*999)
5511 DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
5512 dof_idx=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP% &
5513 & gauss_points(gauss_point_number,element_number)
5514 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
5515 & field_values_set_type,dof_idx,
VALUE,err,error,*999)
5516 stress_tensor(component_idx)=stress_tensor(component_idx)+
VALUE 5523 onethird_trace=sum(stress_tensor(1:3))/3.0_dp
5524 stress_tensor(1:3)=stress_tensor(1:3)-onethird_trace+p
5527 pressure_component=dependent_interpolated_point%INTERPOLATION_PARAMETERS%FIELD_VARIABLE%NUMBER_OF_COMPONENTS
5528 p=dependent_interpolated_point%VALUES(pressure_component,
no_part_deriv)
5529 b=[2.0_dp*c(2),2.0_dp*c(3),2.0_dp*c(3),c(4),c(4),c(3)]
5530 e=[0.5_dp*(azl(1,1)-1.0_dp),0.5_dp*(azl(2,2)-1.0_dp),0.5_dp*(azl(3,3)-1.0_dp),azl(2,1),azl(3,1),azl(3,2)]
5532 tempterm1=0.5_dp*c(1)*exp(0.5_dp*dot_product(e,dq_de))
5534 stress_tensor=tempterm1*dq_de
5539 CALL field_variable_get(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type,field_variable,err,error,*999)
5540 DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS
5541 dof_idx=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP%GAUSS_POINT_PARAM2DOF_MAP% &
5542 & gauss_points(gauss_point_number,element_number)
5543 CALL field_parameter_set_get_local_dof(equations_set%INDEPENDENT%INDEPENDENT_FIELD,field_u_variable_type, &
5544 & field_values_set_type,dof_idx,
VALUE,err,error,*999)
5545 stress_tensor(component_idx)=stress_tensor(component_idx)+
VALUE 5551 onethird_trace=sum(stress_tensor(1:3))/3.0_dp
5552 stress_tensor(1:3)=stress_tensor(1:3)-onethird_trace+p
5554 local_error=
"The third equations set specification of "// &
5556 &
" is not valid for a finite elasticity type of an elasticity equation set." 5557 CALL flagerror(local_error,err,error,*999)
5560 exits(
"FINITE_ELASTICITY_GAUSS_STRESS_TENSOR")
5562 999 errorsexits(
"FINITE_ELASTICITY_GAUSS_STRESS_TENSOR",err,error)
5575 TYPE(
field_type),
POINTER,
INTENT(IN) :: INDEPENDENT_FIELD, MATERIALS_FIELD
5576 REAL(DP),
INTENT(INOUT) :: PIOLA_FF
5577 REAL(DP),
INTENT(IN) :: E_FF
5578 INTEGER(INTG),
INTENT(IN) :: ELEMENT_NUMBER,GAUSS_POINT_NUMBER
5579 INTEGER(INTG),
INTENT(OUT) :: ERR
5583 REAL(DP) :: S, LAMBDA, ISO_TA, TA, ACTIVTIME, TIME, DT
5584 REAL(DP),
DIMENSION(1:4) :: QL
5586 REAL(DP),
PARAMETER :: PERIOD = 1000
5587 REAL(DP),
PARAMETER,
DIMENSION(28) :: TIMES = [ 0, 20, 30, 40, 60, 80, 100, 120, 150, 160, 170, 175, 180, 190, 200,&