From 94843becdafeab9d743d5496556a5b749c15cc3c Mon Sep 17 00:00:00 2001 From: Abisheik Panneerselvam Date: Wed, 14 Jul 2021 19:33:54 +0200 Subject: [PATCH 01/26] writing ip displacements (algorithm contain bugs) --- PRIVATE | 2 +- src/mesh/mesh_mech_FEM.f90 | 43 ++++++++++++++++++++++++++++++++++---- 2 files changed, 40 insertions(+), 5 deletions(-) diff --git a/PRIVATE b/PRIVATE index 1ca9a0e9f..918eed1d9 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 1ca9a0e9f2b333d3244b1ab44480380b3b22bebe +Subproject commit 918eed1d9f67eed75b4a4c66945c8c3053fa10fb diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 6fa2f668b..43b997b70 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -670,7 +670,7 @@ end subroutine FEM_mechanical_converged !-------------------------------------------------------------------------------------------------- -!> @brief Calculate current coordinates (FEM nodal coordinates only at the moment) +!> @brief Calculate current coordinates (both nodal and ip coordinates) !-------------------------------------------------------------------------------------------------- subroutine FEM_mechanical_updateCoords() @@ -678,28 +678,63 @@ subroutine FEM_mechanical_updateCoords() nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes) real(pReal), pointer, dimension(:,:) :: & nodeCoords !< nodal coordinates (3,Nnodes) + real(pReal), pointer, dimension(:,:,:) :: & + ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems) DM :: dm_local Vec :: x_local PetscErrorCode :: ierr - PetscInt :: dimPlex, pStart, pEnd, p, s, e + PetscInt :: dimPlex, pStart, pEnd, p, s, e,& + cellStart, cellEnd, c, qPt, comp, nc, nqpoints,qOffset,nOffset,n,Nnodes PetscSection :: section + PetscFE :: mechFE + PetscDS :: mechDS + PetscQuadrature :: mechQuad + PetscReal, dimension(:), pointer :: qPoints, qWeights, basisField, basisFieldDer + PetscScalar, dimension(:), pointer :: x_scal call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr) + call DMGetDS(dm_local,mechQuad,ierr); CHKERRQ(ierr) call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetDimension(dm_local,dimPlex,ierr); CHKERRQ(ierr) + + ! write cell vertex displacements call DMPlexGetDepthStratum(dm_local,0,pStart,pEnd,ierr); CHKERRQ(ierr) allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal) call VecGetArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr) - do p=pStart, pEnd-1 call DMPlexGetPointLocal(dm_local, p, s, e, ierr); CHKERRQ(ierr) nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e) end do - call discretization_setNodeCoords(nodeCoords) call VecRestoreArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr) + + ! write ip displacements + call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) + call PetscDSGetTabulation(mechQuad,0,basisField,basisFieldDer,ierr); CHKERRQ(ierr) + allocate(ipCoords(0:2,0:nQuadrature-1,0:mesh_NcpElems-1),source=0.0_pReal) + ipCoords=0.0_pReal + do c=cellStart,cellEnd-1 + qOffset=0 + call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) !< get nodal coordinates of each element + Nnodes=size(x_scal)/dimplex !< no. of nodes for each element + do qPt=0,nQuadrature-1 + qOffset=qPt*(cellDof*dimPlex) + do comp=0,dimPlex-1 !< loop over components + nOffset=0 + do n=0,Nnodes-1 + ipCoords(comp,qPt,c)=ipCoords(comp,qPt,c)+sum(basisField(qOffset*dimplex+1:qOffset*dimplex+2)*& + x_scal(nOffset+1:nOffset+2)) + qOffset=qOffset+dimplex + nOffset=nOffset+dimPlex + enddo + qOffset=qPt*(cellDof*dimPlex)+1 + enddo + enddo + call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) + end do + call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature])) call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) end subroutine FEM_mechanical_updateCoords From 4d6a9a51ed8be27ef3b2d81411369a26f5dd5d89 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 15 Jul 2021 17:30:00 +0200 Subject: [PATCH 02/26] correct logic --- PRIVATE | 2 +- src/mesh/mesh_mech_FEM.f90 | 28 +++++++++++++++++----------- 2 files changed, 18 insertions(+), 12 deletions(-) diff --git a/PRIVATE b/PRIVATE index 918eed1d9..8fec909d1 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 918eed1d9f67eed75b4a4c66945c8c3053fa10fb +Subproject commit 8fec909d1931b092b223b0560dd30c3339c6e5a7 diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 43b997b70..c63f92a2c 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -681,10 +681,17 @@ subroutine FEM_mechanical_updateCoords() real(pReal), pointer, dimension(:,:,:) :: & ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems) + integer :: & + qPt, & + comp, & + nc, & + qOffset, & + nOffset + DM :: dm_local Vec :: x_local PetscErrorCode :: ierr - PetscInt :: dimPlex, pStart, pEnd, p, s, e,& + PetscInt :: dimPlex, pStart, pEnd, p, s, e, q, & cellStart, cellEnd, c, qPt, comp, nc, nqpoints,qOffset,nOffset,n,Nnodes PetscSection :: section PetscFE :: mechFE @@ -713,23 +720,22 @@ subroutine FEM_mechanical_updateCoords() ! write ip displacements call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call PetscDSGetTabulation(mechQuad,0,basisField,basisFieldDer,ierr); CHKERRQ(ierr) - allocate(ipCoords(0:2,0:nQuadrature-1,0:mesh_NcpElems-1),source=0.0_pReal) - ipCoords=0.0_pReal + allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pReal) do c=cellStart,cellEnd-1 qOffset=0 call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) !< get nodal coordinates of each element - Nnodes=size(x_scal)/dimplex !< no. of nodes for each element do qPt=0,nQuadrature-1 - qOffset=qPt*(cellDof*dimPlex) + qOffset= qPt * (size(basisField)/nQuadrature) do comp=0,dimPlex-1 !< loop over components nOffset=0 - do n=0,Nnodes-1 - ipCoords(comp,qPt,c)=ipCoords(comp,qPt,c)+sum(basisField(qOffset*dimplex+1:qOffset*dimplex+2)*& - x_scal(nOffset+1:nOffset+2)) - qOffset=qOffset+dimplex - nOffset=nOffset+dimPlex + q = comp + do n=0,nBasis-1 + ipCoords(comp+1,qPt+1,c+1)=ipCoords(comp+1,qPt+1,c+1)+& + sum(basisField(qOffset+(q*dimPlex)+1:qOffset+(q*dimPlex)+dimPlex)*& + x_scal(nOffset+1:nOffset+dimPlex)) + q = q+dimPlex + nOffset = nOffset+dimPlex enddo - qOffset=qPt*(cellDof*dimPlex)+1 enddo enddo call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) From f05ac877eb9eb557dddb8254049c5f1f5f644c03 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 15 Jul 2021 18:53:51 +0200 Subject: [PATCH 03/26] ip displacement added to reference file --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 8fec909d1..feed4bff5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 8fec909d1931b092b223b0560dd30c3339c6e5a7 +Subproject commit feed4bff505c7d556393d35ac9938379eda6d2e4 From 97203ff551e468277f9980a1b7f0e4217edabc66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 18 Jul 2021 23:30:30 +0200 Subject: [PATCH 04/26] document and keep definitions together --- src/mesh/FEM_quadrature.f90 | 7 +++++++ src/prec.f90 | 5 ----- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index ed0bbef9a..3af535c2f 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -20,6 +20,10 @@ module FEM_quadrature -1.0_pReal, 1.0_pReal, -1.0_pReal, & -1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4]) + type :: group_float !< variable length datatype used for storage of state + real(pReal), dimension(:), pointer :: p + end type group_float + integer, dimension(2:3,maxOrder), public, protected :: & FEM_nQuadrature !< number of quadrature points for a given spatial dimension(2-3) and interpolation order(1-maxOrder) type(group_float), dimension(2:3,maxOrder), public, protected :: & @@ -39,6 +43,9 @@ subroutine FEM_quadrature_init print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6) + print*, 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009' + print*, 'https://www.jstor.org/stable/43693493' + !-------------------------------------------------------------------------------------------------- ! 2D linear FEM_nQuadrature(2,1) = 1 diff --git a/src/prec.f90 b/src/prec.f90 index 7613cfe46..733286a8a 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -27,11 +27,6 @@ module prec real(pReal), parameter :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation) - - type :: group_float !< variable length datatype used for storage of state - real(pReal), dimension(:), pointer :: p - end type group_float - type :: tState integer :: & sizeState = 0, & !< size of state From 28179963a292877022aab8161569038e0e33313b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 19 Jul 2021 13:04:16 +0200 Subject: [PATCH 05/26] avoid failing tests --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index feed4bff5..806d65296 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit feed4bff505c7d556393d35ac9938379eda6d2e4 +Subproject commit 806d6529674b21931bb42957492b3f063032febc From a654cd4fb1caa30bff7178986eeb43acf376acf4 Mon Sep 17 00:00:00 2001 From: Abisheik Panneerselvam Date: Mon, 19 Jul 2021 17:07:43 +0200 Subject: [PATCH 06/26] To avoid fortran double variable truncation for quadrature points --- PRIVATE | 2 +- src/mesh/FEM_quadrature.f90 | 25 ++++++++++++++----------- src/mesh/mesh_mech_FEM.f90 | 8 +++----- 3 files changed, 18 insertions(+), 17 deletions(-) diff --git a/PRIVATE b/PRIVATE index 806d65296..2bafc2dbc 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 806d6529674b21931bb42957492b3f063032febc +Subproject commit 2bafc2dbc9b3b6682ef7c77c193dcdfe55d6b9c2 diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 3af535c2f..5aa7ff947 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -71,12 +71,12 @@ subroutine FEM_quadrature_init FEM_nQuadrature(2,3) = 6 allocate(FEM_quadrature_weights(2,3)%p(6)) - FEM_quadrature_weights(2,3)%p(1:3) = 0.22338158967801146570_pReal - FEM_quadrature_weights(2,3)%p(4:6) = 0.10995174365532186764_pReal + FEM_quadrature_weights(2,3)%p(1:3) = 0.2233815896780115_pReal + FEM_quadrature_weights(2,3)%p(4:6) = 0.1099517436553218_pReal allocate(FEM_quadrature_points (2,3)%p(12)) - FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.44594849091596488632_pReal]) - FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.091576213509770743460_pReal]) + FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.4459484909159649_pReal]) + FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.0915762135097707_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quartic @@ -97,6 +97,7 @@ subroutine FEM_quadrature_init FEM_nQuadrature(2,5) = 16 allocate(FEM_quadrature_weights(2,5)%p(16)) + FEM_quadrature_weights(2,5)%p(1 ) = 0.14431560767779_pReal FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728_pReal FEM_quadrature_weights(2,5)%p(5:7) = 0.10321737053472_pReal @@ -104,6 +105,7 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443_pReal allocate(FEM_quadrature_points (2,5)%p(32)) + FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.33333333333333_pReal]) FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.45929258829272_pReal]) FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.17056930775176_pReal]) @@ -128,21 +130,22 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal allocate(FEM_quadrature_points (3,2)%p(12)) - FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.13819660112501051518_pReal]) + + FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.1381966011250105_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D cubic FEM_nQuadrature(3,3) = 14 allocate(FEM_quadrature_weights(3,3)%p(14)) - FEM_quadrature_weights(3,3)%p(5:8) = 0.11268792571801585080_pReal - FEM_quadrature_weights(3,3)%p(1:4) = 0.073493043116361949544_pReal - FEM_quadrature_weights(3,3)%p(9:14) = 0.042546020777081466438_pReal + FEM_quadrature_weights(3,3)%p(5:8) = 0.11268792571801574_pReal + FEM_quadrature_weights(3,3)%p(1:4) = 0.0734930431163620_pReal + FEM_quadrature_weights(3,3)%p(9:14) = 0.0425460207770815_pReal allocate(FEM_quadrature_points (3,3)%p(42)) - FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.092735250310891226402_pReal]) - FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.31088591926330060980_pReal]) - FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.045503704125649649492_pReal]) + FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.0927352503108912_pReal]) + FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.3108859192633006_pReal]) + FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.0455037041256497_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D quartic diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 1bbb3d178..f63ec347f 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -690,13 +690,11 @@ subroutine FEM_mechanical_updateCoords() DM :: dm_local Vec :: x_local PetscErrorCode :: ierr - PetscInt :: dimPlex, pStart, pEnd, p, s, e, q, & - cellStart, cellEnd, c, nqpoints,n + PetscInt :: pStart, pEnd, p, s, e, q, & + cellStart, cellEnd, c, n PetscSection :: section - PetscFE :: mechFE - PetscDS :: mechDS PetscQuadrature :: mechQuad - PetscReal, dimension(:), pointer :: qPoints, qWeights, basisField, basisFieldDer + PetscReal, dimension(:), pointer :: basisField, basisFieldDer PetscScalar, dimension(:), pointer :: x_scal call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr) From 86d918367fe0d4b0ea95b8b3a6b8d3487a4d2bf8 Mon Sep 17 00:00:00 2001 From: Abisheik Panneerselvam Date: Tue, 20 Jul 2021 15:34:13 +0200 Subject: [PATCH 07/26] Round quadrature points --- src/mesh/FEM_quadrature.f90 | 49 ++++++++++++++++++++++--------------- 1 file changed, 29 insertions(+), 20 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 5aa7ff947..32516f4e1 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -71,10 +71,12 @@ subroutine FEM_quadrature_init FEM_nQuadrature(2,3) = 6 allocate(FEM_quadrature_weights(2,3)%p(6)) + FEM_quadrature_weights(2,3)%p(1:3) = 0.2233815896780115_pReal - FEM_quadrature_weights(2,3)%p(4:6) = 0.1099517436553218_pReal + FEM_quadrature_weights(2,3)%p(4:6) = 0.1099517436553219_pReal allocate(FEM_quadrature_points (2,3)%p(12)) + FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.4459484909159649_pReal]) FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.0915762135097707_pReal]) @@ -83,34 +85,35 @@ subroutine FEM_quadrature_init FEM_nQuadrature(2,4) = 12 allocate(FEM_quadrature_weights(2,4)%p(12)) - FEM_quadrature_weights(2,4)%p(1:3) = 0.11678627572638_pReal - FEM_quadrature_weights(2,4)%p(4:6) = 0.05084490637021_pReal - FEM_quadrature_weights(2,4)%p(7:12) = 0.08285107561837_pReal + + FEM_quadrature_weights(2,4)%p(1:3) = 0.1167862757263800_pReal + FEM_quadrature_weights(2,4)%p(4:6) = 0.0508449063702100_pReal + FEM_quadrature_weights(2,4)%p(7:12) = 0.0828510756183700_pReal allocate(FEM_quadrature_points (2,4)%p(24)) - FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.24928674517091_pReal]) - FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.06308901449150_pReal]) - FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal]) + + FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.2492867451709100_pReal]) + FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.0630890144915000_pReal]) + FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.3103524510337800_pReal, 0.6365024991214000_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quintic FEM_nQuadrature(2,5) = 16 allocate(FEM_quadrature_weights(2,5)%p(16)) - - FEM_quadrature_weights(2,5)%p(1 ) = 0.14431560767779_pReal - FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728_pReal - FEM_quadrature_weights(2,5)%p(5:7) = 0.10321737053472_pReal - FEM_quadrature_weights(2,5)%p(8:10) = 0.03245849762320_pReal - FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443_pReal + FEM_quadrature_weights(2,5)%p(1 ) = 0.1443156076777900_pReal + FEM_quadrature_weights(2,5)%p(2:4) = 0.0950916342672800_pReal + FEM_quadrature_weights(2,5)%p(5:7) = 0.1032173705347200_pReal + FEM_quadrature_weights(2,5)%p(8:10) = 0.0324584976232000_pReal + FEM_quadrature_weights(2,5)%p(11:16)= 0.0272303141744300_pReal allocate(FEM_quadrature_points (2,5)%p(32)) - FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.33333333333333_pReal]) - FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.45929258829272_pReal]) - FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.17056930775176_pReal]) - FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.05054722831703_pReal]) - FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal]) + FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.3333333333333300_pReal]) + FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.4592925882927200_pReal]) + FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.1705693077517600_pReal]) + FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.0505472283170300_pReal]) + FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.2631128296346400_pReal, 0.7284923929554000_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D linear @@ -138,11 +141,13 @@ subroutine FEM_quadrature_init FEM_nQuadrature(3,3) = 14 allocate(FEM_quadrature_weights(3,3)%p(14)) - FEM_quadrature_weights(3,3)%p(5:8) = 0.11268792571801574_pReal + + FEM_quadrature_weights(3,3)%p(5:8) = 0.1126879257180159_pReal FEM_quadrature_weights(3,3)%p(1:4) = 0.0734930431163620_pReal FEM_quadrature_weights(3,3)%p(9:14) = 0.0425460207770815_pReal allocate(FEM_quadrature_points (3,3)%p(42)) + FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.0927352503108912_pReal]) FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.3108859192633006_pReal]) FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.0455037041256497_pReal]) @@ -152,6 +157,7 @@ subroutine FEM_quadrature_init FEM_nQuadrature(3,4) = 35 allocate(FEM_quadrature_weights(3,4)%p(35)) + FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal @@ -159,6 +165,7 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal allocate(FEM_quadrature_points (3,4)%p(105)) + FEM_quadrature_points (3,4)%p(1:12) = permutationStar31([0.0267367755543735_pReal]) FEM_quadrature_points (3,4)%p(13:48) = permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]) FEM_quadrature_points (3,4)%p(49:66) = permutationStar22([0.4547545999844830_pReal]) @@ -170,6 +177,7 @@ subroutine FEM_quadrature_init FEM_nQuadrature(3,5) = 56 allocate(FEM_quadrature_weights(3,5)%p(56)) + FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal @@ -178,6 +186,7 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal allocate(FEM_quadrature_points (3,5)%p(168)) + FEM_quadrature_points (3,5)%p(1:12) = permutationStar31([0.0149520651530592_pReal]) FEM_quadrature_points (3,5)%p(13:48) = permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]) FEM_quadrature_points (3,5)%p(49:84) = permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]) @@ -372,4 +381,4 @@ pure function permutationStar1111(point) result(qPt) end function permutationStar1111 -end module FEM_quadrature +end module FEM_quadrature \ No newline at end of file From 9427885b43b8b6ac66385de5da6c1e09d24a8005 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 21 Jul 2021 13:29:59 +0200 Subject: [PATCH 08/26] updated private --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 2bafc2dbc..4d696af2b 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2bafc2dbc9b3b6682ef7c77c193dcdfe55d6b9c2 +Subproject commit 4d696af2bf78bb6d707fd51e958843b42b1279ff From f25dee4c3a7f1ddc73a577072fe7f2a87b78a50f Mon Sep 17 00:00:00 2001 From: Abisheik Panneerselvam Date: Wed, 21 Jul 2021 14:34:18 +0200 Subject: [PATCH 09/26] quad points updated and rounded based on literature values --- PRIVATE | 2 +- src/mesh/FEM_quadrature.f90 | 48 +++++++++++++++++-------------------- 2 files changed, 23 insertions(+), 27 deletions(-) diff --git a/PRIVATE b/PRIVATE index 4d696af2b..567067daf 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 4d696af2bf78bb6d707fd51e958843b42b1279ff +Subproject commit 567067daf05866c8e022fc0af6f441840f143b81 diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 32516f4e1..e310cae06 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -70,50 +70,46 @@ subroutine FEM_quadrature_init ! 2D cubic FEM_nQuadrature(2,3) = 6 - allocate(FEM_quadrature_weights(2,3)%p(6)) - + allocate(FEM_quadrature_weights(2,3)%p(6)) FEM_quadrature_weights(2,3)%p(1:3) = 0.2233815896780115_pReal FEM_quadrature_weights(2,3)%p(4:6) = 0.1099517436553219_pReal allocate(FEM_quadrature_points (2,3)%p(12)) - FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.4459484909159649_pReal]) - FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.0915762135097707_pReal]) + FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.09157621350977074_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quartic FEM_nQuadrature(2,4) = 12 - allocate(FEM_quadrature_weights(2,4)%p(12)) - - FEM_quadrature_weights(2,4)%p(1:3) = 0.1167862757263800_pReal - FEM_quadrature_weights(2,4)%p(4:6) = 0.0508449063702100_pReal - FEM_quadrature_weights(2,4)%p(7:12) = 0.0828510756183700_pReal + allocate(FEM_quadrature_weights(2,4)%p(12)) + FEM_quadrature_weights(2,4)%p(1:3) = 0.1167862757263794_pReal + FEM_quadrature_weights(2,4)%p(4:6) = 0.0508449063702068_pReal + FEM_quadrature_weights(2,4)%p(7:12) = 0.08285107561837358_pReal - allocate(FEM_quadrature_points (2,4)%p(24)) - - FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.2492867451709100_pReal]) - FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.0630890144915000_pReal]) - FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.3103524510337800_pReal, 0.6365024991214000_pReal]) + allocate(FEM_quadrature_points (2,4)%p(24)) + FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.2492867451709104_pReal]) + FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.06308901449150223_pReal]) + FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.3103524510337844_pReal, 0.05314504984481695_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quintic FEM_nQuadrature(2,5) = 16 allocate(FEM_quadrature_weights(2,5)%p(16)) - FEM_quadrature_weights(2,5)%p(1 ) = 0.1443156076777900_pReal - FEM_quadrature_weights(2,5)%p(2:4) = 0.0950916342672800_pReal - FEM_quadrature_weights(2,5)%p(5:7) = 0.1032173705347200_pReal - FEM_quadrature_weights(2,5)%p(8:10) = 0.0324584976232000_pReal - FEM_quadrature_weights(2,5)%p(11:16)= 0.0272303141744300_pReal + FEM_quadrature_weights(2,5)%p(1 ) = 0.1443156076777871_pReal + FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728462_pReal + FEM_quadrature_weights(2,5)%p(5:7) = 0.1032173705347183_pReal + FEM_quadrature_weights(2,5)%p(8:10) = 0.03245849762319808_pReal + FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443499_pReal allocate(FEM_quadrature_points (2,5)%p(32)) - FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.3333333333333300_pReal]) - FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.4592925882927200_pReal]) - FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.1705693077517600_pReal]) - FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.0505472283170300_pReal]) - FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.2631128296346400_pReal, 0.7284923929554000_pReal]) + FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.3333333333333333_pReal]) + FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.4592925882927231_pReal]) + FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.1705693077517602_pReal]) + FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.0505472283170310_pReal]) + FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.2631128296346381_pReal, 0.008394777409957605_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D linear @@ -148,9 +144,9 @@ subroutine FEM_quadrature_init allocate(FEM_quadrature_points (3,3)%p(42)) - FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.0927352503108912_pReal]) + FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.09273525031089123_pReal]) FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.3108859192633006_pReal]) - FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.0455037041256497_pReal]) + FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.04550370412564965_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D quartic From b1ba64e6afbc7ba9e2b013c2750ccb7ef07d3d87 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 17:51:58 +0200 Subject: [PATCH 10/26] understanding floating point precision double can use 16 or 17 significant figures - if the leading number is high (above 5 or so), only 16 figures are significant. - a trailing 1 is never significant --- src/mesh/FEM_quadrature.f90 | 70 ++++++++++++++++--------------------- 1 file changed, 31 insertions(+), 39 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index e310cae06..f99bc58b3 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -51,65 +51,64 @@ subroutine FEM_quadrature_init FEM_nQuadrature(2,1) = 1 allocate(FEM_quadrature_weights(2,1)%p(1)) - FEM_quadrature_weights(2,1)%p(1) = 1.0_pReal + FEM_quadrature_weights(2,1)%p(1) = 1._pReal allocate(FEM_quadrature_points (2,1)%p(2)) - FEM_quadrature_points (2,1)%p(1:2) = permutationStar3([1.0_pReal/3.0_pReal]) + FEM_quadrature_points (2,1)%p(1:2) = permutationStar3([1._pReal/3._pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quadratic FEM_nQuadrature(2,2) = 3 allocate(FEM_quadrature_weights(2,2)%p(3)) - FEM_quadrature_weights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal + FEM_quadrature_weights(2,2)%p(1:3) = 1._pReal/3._pReal allocate(FEM_quadrature_points (2,2)%p(6)) - FEM_quadrature_points (2,2)%p(1:6) = permutationStar21([1.0_pReal/6.0_pReal]) + FEM_quadrature_points (2,2)%p(1:6) = permutationStar21([1._pReal/6._pReal]) !-------------------------------------------------------------------------------------------------- ! 2D cubic FEM_nQuadrature(2,3) = 6 allocate(FEM_quadrature_weights(2,3)%p(6)) - FEM_quadrature_weights(2,3)%p(1:3) = 0.2233815896780115_pReal - FEM_quadrature_weights(2,3)%p(4:6) = 0.1099517436553219_pReal + FEM_quadrature_weights(2,3)%p(1:3) = 2.2338158967801147e-1_pReal + FEM_quadrature_weights(2,3)%p(4:6) = 1.0995174365532187e-1_pReal allocate(FEM_quadrature_points (2,3)%p(12)) - FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.4459484909159649_pReal]) - FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.09157621350977074_pReal]) + FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([4.4594849091596489e-1_pReal]) + FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([9.157621350977074e-2_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quartic FEM_nQuadrature(2,4) = 12 allocate(FEM_quadrature_weights(2,4)%p(12)) - FEM_quadrature_weights(2,4)%p(1:3) = 0.1167862757263794_pReal - FEM_quadrature_weights(2,4)%p(4:6) = 0.0508449063702068_pReal - FEM_quadrature_weights(2,4)%p(7:12) = 0.08285107561837358_pReal + FEM_quadrature_weights(2,4)%p(1:3) = 1.1678627572637937e-1_pReal + FEM_quadrature_weights(2,4)%p(4:6) = 5.0844906370206817e-2_pReal + FEM_quadrature_weights(2,4)%p(7:12) = 8.285107561837358e-2_pReal allocate(FEM_quadrature_points (2,4)%p(24)) - FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.2492867451709104_pReal]) - FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.06308901449150223_pReal]) - FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.3103524510337844_pReal, 0.05314504984481695_pReal]) + FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([2.4928674517091042e-1_pReal]) + FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([6.308901449150223e-2_pReal]) + FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([3.1035245103378440e-1_pReal, 5.3145049844816947e-2_pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quintic FEM_nQuadrature(2,5) = 16 allocate(FEM_quadrature_weights(2,5)%p(16)) - FEM_quadrature_weights(2,5)%p(1 ) = 0.1443156076777871_pReal - FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728462_pReal - FEM_quadrature_weights(2,5)%p(5:7) = 0.1032173705347183_pReal - FEM_quadrature_weights(2,5)%p(8:10) = 0.03245849762319808_pReal - FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443499_pReal + FEM_quadrature_weights(2,5)%p(1 ) = 1.4431560767778717e-1_pReal + FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal + FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal + FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal + FEM_quadrature_weights(2,5)%p(11:16)= 2.7230314174434994e-2_pReal allocate(FEM_quadrature_points (2,5)%p(32)) - - FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.3333333333333333_pReal]) - FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.4592925882927231_pReal]) - FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.1705693077517602_pReal]) - FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.0505472283170310_pReal]) - FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.2631128296346381_pReal, 0.008394777409957605_pReal]) + FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([1._pReal/3._pReal]) + FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([4.5929258829272316e-1_pReal]) + FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([1.705693077517602e-1_pReal]) + FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([5.0547228317030975e-2_pReal]) + FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([2.631128296346381e-1_pReal, 8.3947774099576053e-2_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D linear @@ -129,7 +128,6 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal allocate(FEM_quadrature_points (3,2)%p(12)) - FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.1381966011250105_pReal]) !-------------------------------------------------------------------------------------------------- @@ -137,23 +135,20 @@ subroutine FEM_quadrature_init FEM_nQuadrature(3,3) = 14 allocate(FEM_quadrature_weights(3,3)%p(14)) - - FEM_quadrature_weights(3,3)%p(5:8) = 0.1126879257180159_pReal - FEM_quadrature_weights(3,3)%p(1:4) = 0.0734930431163620_pReal - FEM_quadrature_weights(3,3)%p(9:14) = 0.0425460207770815_pReal + FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal + FEM_quadrature_weights(3,3)%p(1:4) = 7.3493043116361949e-2_pReal + FEM_quadrature_weights(3,3)%p(9:14) = 4.2546020777081467e-2_pReal allocate(FEM_quadrature_points (3,3)%p(42)) - - FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.09273525031089123_pReal]) - FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.3108859192633006_pReal]) - FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.04550370412564965_pReal]) + FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([9.273525031089123e-2_pReal]) + FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([3.108859192633006e-1_pReal]) + FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([4.5503704125649649e-2_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D quartic FEM_nQuadrature(3,4) = 35 allocate(FEM_quadrature_weights(3,4)%p(35)) - FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal @@ -161,7 +156,6 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal allocate(FEM_quadrature_points (3,4)%p(105)) - FEM_quadrature_points (3,4)%p(1:12) = permutationStar31([0.0267367755543735_pReal]) FEM_quadrature_points (3,4)%p(13:48) = permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]) FEM_quadrature_points (3,4)%p(49:66) = permutationStar22([0.4547545999844830_pReal]) @@ -173,7 +167,6 @@ subroutine FEM_quadrature_init FEM_nQuadrature(3,5) = 56 allocate(FEM_quadrature_weights(3,5)%p(56)) - FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal @@ -182,7 +175,6 @@ subroutine FEM_quadrature_init FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal allocate(FEM_quadrature_points (3,5)%p(168)) - FEM_quadrature_points (3,5)%p(1:12) = permutationStar31([0.0149520651530592_pReal]) FEM_quadrature_points (3,5)%p(13:48) = permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]) FEM_quadrature_points (3,5)%p(49:84) = permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]) @@ -377,4 +369,4 @@ pure function permutationStar1111(point) result(qPt) end function permutationStar1111 -end module FEM_quadrature \ No newline at end of file +end module FEM_quadrature From 000007ba596390ac52427e812805567897865cb5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 18:06:37 +0200 Subject: [PATCH 11/26] better test automatically --- src/mesh/FEM_quadrature.f90 | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index f99bc58b3..31d4acec0 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -39,7 +39,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief initializes FEM interpolation data !-------------------------------------------------------------------------------------------------- -subroutine FEM_quadrature_init +subroutine FEM_quadrature_init() print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6) @@ -182,6 +182,8 @@ subroutine FEM_quadrature_init FEM_quadrature_points (3,5)%p(121:156)= permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]) FEM_quadrature_points (3,5)%p(157:168)= permutationStar31([0.1344783347929940_pReal]) + call selfTest + end subroutine FEM_quadrature_init @@ -369,4 +371,21 @@ pure function permutationStar1111(point) result(qPt) end function permutationStar1111 + +!-------------------------------------------------------------------------------------------------- +!> @brief Check correctness of quadrature weights. +!-------------------------------------------------------------------------------------------------- +subroutine selfTest + + integer :: o, d + + + do d = lbound(FEM_quadrature_weights,1), ubound(FEM_quadrature_weights,1) + do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1) + if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) error stop 'quadrature weights' + enddo + enddo + +end subroutine selfTest + end module FEM_quadrature From c56979e2ad7d1b08accc75081ee98fc3e4840668 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 18:09:44 +0200 Subject: [PATCH 12/26] single source of truth --- src/mesh/FEM_quadrature.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 31d4acec0..2498e3947 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -50,7 +50,7 @@ subroutine FEM_quadrature_init() ! 2D linear FEM_nQuadrature(2,1) = 1 - allocate(FEM_quadrature_weights(2,1)%p(1)) + allocate(FEM_quadrature_weights(2,1)%p(FEM_nQuadrature(2,1))) FEM_quadrature_weights(2,1)%p(1) = 1._pReal allocate(FEM_quadrature_points (2,1)%p(2)) @@ -60,7 +60,7 @@ subroutine FEM_quadrature_init() ! 2D quadratic FEM_nQuadrature(2,2) = 3 - allocate(FEM_quadrature_weights(2,2)%p(3)) + allocate(FEM_quadrature_weights(2,2)%p(FEM_nQuadrature(2,2))) FEM_quadrature_weights(2,2)%p(1:3) = 1._pReal/3._pReal allocate(FEM_quadrature_points (2,2)%p(6)) @@ -70,7 +70,7 @@ subroutine FEM_quadrature_init() ! 2D cubic FEM_nQuadrature(2,3) = 6 - allocate(FEM_quadrature_weights(2,3)%p(6)) + allocate(FEM_quadrature_weights(2,3)%p(FEM_nQuadrature(2,3))) FEM_quadrature_weights(2,3)%p(1:3) = 2.2338158967801147e-1_pReal FEM_quadrature_weights(2,3)%p(4:6) = 1.0995174365532187e-1_pReal @@ -82,7 +82,7 @@ subroutine FEM_quadrature_init() ! 2D quartic FEM_nQuadrature(2,4) = 12 - allocate(FEM_quadrature_weights(2,4)%p(12)) + allocate(FEM_quadrature_weights(2,4)%p(FEM_nQuadrature(2,4))) FEM_quadrature_weights(2,4)%p(1:3) = 1.1678627572637937e-1_pReal FEM_quadrature_weights(2,4)%p(4:6) = 5.0844906370206817e-2_pReal FEM_quadrature_weights(2,4)%p(7:12) = 8.285107561837358e-2_pReal @@ -96,7 +96,7 @@ subroutine FEM_quadrature_init() ! 2D quintic FEM_nQuadrature(2,5) = 16 - allocate(FEM_quadrature_weights(2,5)%p(16)) + allocate(FEM_quadrature_weights(2,5)%p(FEM_nQuadrature(2,5))) FEM_quadrature_weights(2,5)%p(1 ) = 1.4431560767778717e-1_pReal FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal @@ -114,7 +114,7 @@ subroutine FEM_quadrature_init() ! 3D linear FEM_nQuadrature(3,1) = 1 - allocate(FEM_quadrature_weights(3,1)%p(1)) + allocate(FEM_quadrature_weights(3,1)%p(FEM_nQuadrature(3,1))) FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal allocate(FEM_quadrature_points (3,1)%p(3)) @@ -124,7 +124,7 @@ subroutine FEM_quadrature_init() ! 3D quadratic FEM_nQuadrature(3,2) = 4 - allocate(FEM_quadrature_weights(3,2)%p(4)) + allocate(FEM_quadrature_weights(3,2)%p(FEM_nQuadrature(3,2))) FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal allocate(FEM_quadrature_points (3,2)%p(12)) @@ -134,7 +134,7 @@ subroutine FEM_quadrature_init() ! 3D cubic FEM_nQuadrature(3,3) = 14 - allocate(FEM_quadrature_weights(3,3)%p(14)) + allocate(FEM_quadrature_weights(3,3)%p(FEM_nQuadrature(3,3))) FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal FEM_quadrature_weights(3,3)%p(1:4) = 7.3493043116361949e-2_pReal FEM_quadrature_weights(3,3)%p(9:14) = 4.2546020777081467e-2_pReal @@ -148,7 +148,7 @@ subroutine FEM_quadrature_init() ! 3D quartic FEM_nQuadrature(3,4) = 35 - allocate(FEM_quadrature_weights(3,4)%p(35)) + allocate(FEM_quadrature_weights(3,4)%p(FEM_nQuadrature(3,4))) FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal @@ -166,7 +166,7 @@ subroutine FEM_quadrature_init() ! 3D quintic FEM_nQuadrature(3,5) = 56 - allocate(FEM_quadrature_weights(3,5)%p(56)) + allocate(FEM_quadrature_weights(3,5)%p(FEM_nQuadrature(3,5))) FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal From 5a04a1d661a384c8469f9b22bc3f03e9b5a40f15 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 18:59:34 +0200 Subject: [PATCH 13/26] let the computer do the counting --- src/mesh/FEM_quadrature.f90 | 104 +++++++++++++++++------------------- 1 file changed, 50 insertions(+), 54 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 2498e3947..e76e9849f 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -9,7 +9,7 @@ module FEM_quadrature private integer, parameter :: & - maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary) + maxOrder = 5 !< maximum integration order real(pReal), dimension(2,3), parameter :: & triangle = reshape([-1.0_pReal, -1.0_pReal, & 1.0_pReal, -1.0_pReal, & @@ -21,7 +21,7 @@ module FEM_quadrature -1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4]) type :: group_float !< variable length datatype used for storage of state - real(pReal), dimension(:), pointer :: p + real(pReal), dimension(:), allocatable :: p end type group_float integer, dimension(2:3,maxOrder), public, protected :: & @@ -53,8 +53,7 @@ subroutine FEM_quadrature_init() allocate(FEM_quadrature_weights(2,1)%p(FEM_nQuadrature(2,1))) FEM_quadrature_weights(2,1)%p(1) = 1._pReal - allocate(FEM_quadrature_points (2,1)%p(2)) - FEM_quadrature_points (2,1)%p(1:2) = permutationStar3([1._pReal/3._pReal]) + FEM_quadrature_points (2,1)%p = permutationStar3([1._pReal/3._pReal]) !-------------------------------------------------------------------------------------------------- ! 2D quadratic @@ -63,8 +62,7 @@ subroutine FEM_quadrature_init() allocate(FEM_quadrature_weights(2,2)%p(FEM_nQuadrature(2,2))) FEM_quadrature_weights(2,2)%p(1:3) = 1._pReal/3._pReal - allocate(FEM_quadrature_points (2,2)%p(6)) - FEM_quadrature_points (2,2)%p(1:6) = permutationStar21([1._pReal/6._pReal]) + FEM_quadrature_points (2,2)%p = permutationStar21([1._pReal/6._pReal]) !-------------------------------------------------------------------------------------------------- ! 2D cubic @@ -74,9 +72,9 @@ subroutine FEM_quadrature_init() FEM_quadrature_weights(2,3)%p(1:3) = 2.2338158967801147e-1_pReal FEM_quadrature_weights(2,3)%p(4:6) = 1.0995174365532187e-1_pReal - allocate(FEM_quadrature_points (2,3)%p(12)) - FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([4.4594849091596489e-1_pReal]) - FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([9.157621350977074e-2_pReal]) + FEM_quadrature_points (2,3)%p = [ & + permutationStar21([4.4594849091596489e-1_pReal]), & + permutationStar21([9.157621350977074e-2_pReal]) ] !-------------------------------------------------------------------------------------------------- ! 2D quartic @@ -87,10 +85,10 @@ subroutine FEM_quadrature_init() FEM_quadrature_weights(2,4)%p(4:6) = 5.0844906370206817e-2_pReal FEM_quadrature_weights(2,4)%p(7:12) = 8.285107561837358e-2_pReal - allocate(FEM_quadrature_points (2,4)%p(24)) - FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([2.4928674517091042e-1_pReal]) - FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([6.308901449150223e-2_pReal]) - FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([3.1035245103378440e-1_pReal, 5.3145049844816947e-2_pReal]) + FEM_quadrature_points (2,4)%p = [ & + permutationStar21([2.4928674517091042e-1_pReal]), & + permutationStar21([6.308901449150223e-2_pReal]), & + permutationStar111([3.1035245103378440e-1_pReal, 5.3145049844816947e-2_pReal]) ] !-------------------------------------------------------------------------------------------------- ! 2D quintic @@ -103,22 +101,21 @@ subroutine FEM_quadrature_init() FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal FEM_quadrature_weights(2,5)%p(11:16)= 2.7230314174434994e-2_pReal - allocate(FEM_quadrature_points (2,5)%p(32)) - FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([1._pReal/3._pReal]) - FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([4.5929258829272316e-1_pReal]) - FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([1.705693077517602e-1_pReal]) - FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([5.0547228317030975e-2_pReal]) - FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([2.631128296346381e-1_pReal, 8.3947774099576053e-2_pReal]) + FEM_quadrature_points (2,5)%p = [ & + permutationStar3([1._pReal/3._pReal]), & + permutationStar21([4.5929258829272316e-1_pReal]), & + permutationStar21([1.705693077517602e-1_pReal]), & + permutationStar21([5.0547228317030975e-2_pReal]), & + permutationStar111([2.631128296346381e-1_pReal, 8.3947774099576053e-2_pReal]) ] !-------------------------------------------------------------------------------------------------- ! 3D linear FEM_nQuadrature(3,1) = 1 allocate(FEM_quadrature_weights(3,1)%p(FEM_nQuadrature(3,1))) - FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal + FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal - allocate(FEM_quadrature_points (3,1)%p(3)) - FEM_quadrature_points (3,1)%p(1:3)= permutationStar4([0.25_pReal]) + FEM_quadrature_points (3,1)%p = permutationStar4([0.25_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D quadratic @@ -127,60 +124,59 @@ subroutine FEM_quadrature_init() allocate(FEM_quadrature_weights(3,2)%p(FEM_nQuadrature(3,2))) FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal - allocate(FEM_quadrature_points (3,2)%p(12)) - FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.1381966011250105_pReal]) + FEM_quadrature_points (3,2)%p = permutationStar31([0.1381966011250105_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D cubic FEM_nQuadrature(3,3) = 14 allocate(FEM_quadrature_weights(3,3)%p(FEM_nQuadrature(3,3))) - FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal FEM_quadrature_weights(3,3)%p(1:4) = 7.3493043116361949e-2_pReal + FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal FEM_quadrature_weights(3,3)%p(9:14) = 4.2546020777081467e-2_pReal - allocate(FEM_quadrature_points (3,3)%p(42)) - FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([9.273525031089123e-2_pReal]) - FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([3.108859192633006e-1_pReal]) - FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([4.5503704125649649e-2_pReal]) + FEM_quadrature_points (3,3)%p = [ & + permutationStar31([9.273525031089123e-2_pReal]), & + permutationStar31([3.108859192633006e-1_pReal]), & + permutationStar22([4.5503704125649649e-2_pReal]) ] !-------------------------------------------------------------------------------------------------- ! 3D quartic FEM_nQuadrature(3,4) = 35 allocate(FEM_quadrature_weights(3,4)%p(FEM_nQuadrature(3,4))) - FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal - FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal - FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal - FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal - FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal + FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal + FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal + FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal + FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal + FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal - allocate(FEM_quadrature_points (3,4)%p(105)) - FEM_quadrature_points (3,4)%p(1:12) = permutationStar31([0.0267367755543735_pReal]) - FEM_quadrature_points (3,4)%p(13:48) = permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]) - FEM_quadrature_points (3,4)%p(49:66) = permutationStar22([0.4547545999844830_pReal]) - FEM_quadrature_points (3,4)%p(67:102) = permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal]) - FEM_quadrature_points (3,4)%p(103:105)= permutationStar4([0.25_pReal]) + FEM_quadrature_points (3,4)%p = [ & + permutationStar31([0.0267367755543735_pReal]), & + permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]), & + permutationStar22([0.4547545999844830_pReal]), & + permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal]), & + permutationStar4([0.25_pReal]) ] !-------------------------------------------------------------------------------------------------- ! 3D quintic FEM_nQuadrature(3,5) = 56 allocate(FEM_quadrature_weights(3,5)%p(FEM_nQuadrature(3,5))) - FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal - FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal - FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal - FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal - FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal - FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal + FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal + FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal + FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal + FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal + FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal + FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal - allocate(FEM_quadrature_points (3,5)%p(168)) - FEM_quadrature_points (3,5)%p(1:12) = permutationStar31([0.0149520651530592_pReal]) - FEM_quadrature_points (3,5)%p(13:48) = permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]) - FEM_quadrature_points (3,5)%p(49:84) = permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]) - FEM_quadrature_points (3,5)%p(85:120) = permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal]) - FEM_quadrature_points (3,5)%p(121:156)= permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]) - FEM_quadrature_points (3,5)%p(157:168)= permutationStar31([0.1344783347929940_pReal]) + FEM_quadrature_points (3,5)%p = [ & + permutationStar31([0.0149520651530592_pReal]), & + permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]), & + permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]), & + permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal]), & + permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]), & + permutationStar31([0.1344783347929940_pReal]) ] call selfTest @@ -218,7 +214,7 @@ pure function permutationStar21(point) result(qPt) temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)] temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)] - qPt = reshape(matmul(triangle, temp),[6]) + qPt = reshape(matmul(triangle, temp),[6]) end function permutationStar21 From 23077fd58c13f8f584a8e8ff30b6b5fa212ab1ce Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 19:42:53 +0200 Subject: [PATCH 14/26] testing coordinates I don't understand why this pattern exists (with w(2) = 3, w(3) = 2), but it exists --- src/mesh/FEM_quadrature.f90 | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index e76e9849f..53aadbe6f 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -373,12 +373,21 @@ end function permutationStar1111 !-------------------------------------------------------------------------------------------------- subroutine selfTest - integer :: o, d - + integer :: o, d, n + real(pReal), dimension(2:3), parameter :: w = [3.0_pReal,2.0_pReal] do d = lbound(FEM_quadrature_weights,1), ubound(FEM_quadrature_weights,1) do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1) - if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) error stop 'quadrature weights' + if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) & + error stop 'quadrature weights' + enddo + enddo + + do d = lbound(FEM_quadrature_points,1), ubound(FEM_quadrature_points,1) + do o = lbound(FEM_quadrature_points(d,:),1), ubound(FEM_quadrature_points(d,:),1) + n = size(FEM_quadrature_points(d,o)%p,1)/d + if (any(dNeq(sum(reshape(FEM_quadrature_points(d,o)%p,[d,n]),2),-real(n,pReal)/w(d),1.e-14_pReal))) & + error stop 'quadrature points' enddo enddo From b817c620a371e7dcb32b9172b756a5bd59e83879 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 20:35:42 +0200 Subject: [PATCH 15/26] let the computer count --- src/mesh/FEM_quadrature.f90 | 147 ++++++++++++++++++++---------------- 1 file changed, 82 insertions(+), 65 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 53aadbe6f..278bfb5db 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -51,7 +51,7 @@ subroutine FEM_quadrature_init() FEM_nQuadrature(2,1) = 1 allocate(FEM_quadrature_weights(2,1)%p(FEM_nQuadrature(2,1))) - FEM_quadrature_weights(2,1)%p(1) = 1._pReal + FEM_quadrature_weights(2,1)%p(1) = 1._pReal FEM_quadrature_points (2,1)%p = permutationStar3([1._pReal/3._pReal]) @@ -95,11 +95,11 @@ subroutine FEM_quadrature_init() FEM_nQuadrature(2,5) = 16 allocate(FEM_quadrature_weights(2,5)%p(FEM_nQuadrature(2,5))) - FEM_quadrature_weights(2,5)%p(1 ) = 1.4431560767778717e-1_pReal - FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal - FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal - FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal - FEM_quadrature_weights(2,5)%p(11:16)= 2.7230314174434994e-2_pReal + FEM_quadrature_weights(2,5)%p(1:1) = 1.4431560767778717e-1_pReal + FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal + FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal + FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal + FEM_quadrature_weights(2,5)%p(11:16) = 2.7230314174434994e-2_pReal FEM_quadrature_points (2,5)%p = [ & permutationStar3([1._pReal/3._pReal]), & @@ -193,7 +193,9 @@ pure function permutationStar3(point) result(qPt) real(pReal), dimension(3,1) :: temp - temp(:,1) = [point(1), point(1), point(1)] + + temp = reshape([ & + point(1), point(1), point(1)],shape(temp)) qPt = reshape(matmul(triangle, temp),[2]) @@ -210,9 +212,11 @@ pure function permutationStar21(point) result(qPt) real(pReal), dimension(3,3) :: temp - temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)] - temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)] - temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)] + + temp = reshape([ & + point(1), point(1), 1.0_pReal - 2.0_pReal*point(1), & + point(1), 1.0_pReal - 2.0_pReal*point(1), point(1), & + 1.0_pReal - 2.0_pReal*point(1), point(1), point(1)],shape(temp)) qPt = reshape(matmul(triangle, temp),[6]) @@ -229,12 +233,14 @@ pure function permutationStar111(point) result(qPt) real(pReal), dimension(3,6) :: temp - temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)] - temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)] - temp(:,3) = [point(2), point(1), 1.0_pReal - point(1) - point(2)] - temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)] - temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)] - temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)] + + temp = reshape([ & + point(1), point(2), 1.0_pReal - point(1) - point(2), & + point(1), 1.0_pReal - point(1) - point(2), point(2), & + point(2), point(1), 1.0_pReal - point(1) - point(2), & + point(2), 1.0_pReal - point(1) - point(2), point(1), & + 1.0_pReal - point(1) - point(2), point(2), point(1), & + 1.0_pReal - point(1) - point(2), point(1), point(2)],shape(temp)) qPt = reshape(matmul(triangle, temp),[12]) @@ -251,7 +257,9 @@ pure function permutationStar4(point) result(qPt) real(pReal), dimension(4,1) :: temp - temp(:,1) = [point(1), point(1), point(1), point(1)] + + temp = reshape([ & + point(1), point(1), point(1), point(1)],shape(temp)) qPt = reshape(matmul(tetrahedron, temp),[3]) @@ -268,10 +276,12 @@ pure function permutationStar31(point) result(qPt) real(pReal), dimension(4,4) :: temp - temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)] - temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)] - temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)] - temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)] + + temp = reshape([ & + point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), & + point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), & + point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), & + 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)],shape(temp)) qPt = reshape(matmul(tetrahedron, temp),[12]) @@ -281,19 +291,21 @@ end function permutationStar31 !-------------------------------------------------------------------------------------------------- !> @brief star 22 permutation of input !-------------------------------------------------------------------------------------------------- -pure function permutationStar22(point) result(qPt) +function permutationStar22(point) result(qPt) real(pReal), dimension(18) :: qPt real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(4,6) :: temp - temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)] - temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)] - temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)] - temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)] - temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)] - temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)] + + temp = reshape([ & + point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), & + point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), & + 0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1), & + 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1), & + 0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1), & + point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)],shape(temp)) qPt = reshape(matmul(tetrahedron, temp),[18]) @@ -310,18 +322,20 @@ pure function permutationStar211(point) result(qPt) real(pReal), dimension(4,12) :: temp - temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)] - temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)] - temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)] - temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)] - temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)] - temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)] - temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)] - temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)] - temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)] - temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)] - temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)] - temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)] + + temp = reshape([ & + point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), & + point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), & + point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), & + point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), & + point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), & + point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), & + point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), & + point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), & + point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), & + 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2), & + 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1), & + 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)],shape(temp)) qPt = reshape(matmul(tetrahedron, temp),[36]) @@ -338,30 +352,32 @@ pure function permutationStar1111(point) result(qPt) real(pReal), dimension(4,24) :: temp - temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)] - temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)] - temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)] - temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)] - temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)] - temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)] - temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)] - temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)] - temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)] - temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)] - temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)] - temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)] - temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)] - temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)] - temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)] - temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)] - temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)] - temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)] - temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)] + + temp = reshape([ & + point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), & + point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), & + point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), & + point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), & + point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), & + point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), & + point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), & + point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), & + point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), & + point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), & + point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), & + point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), & + point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), & + point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), & + point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), & + point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), & + point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), & + point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), & + 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3), & + 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2), & + 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3), & + 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1), & + 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2), & + 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)],shape(temp)) qPt = reshape(matmul(tetrahedron, temp),[72]) @@ -369,13 +385,14 @@ end function permutationStar1111 !-------------------------------------------------------------------------------------------------- -!> @brief Check correctness of quadrature weights. +!> @brief Check correctness of quadrature weights and points. !-------------------------------------------------------------------------------------------------- subroutine selfTest integer :: o, d, n real(pReal), dimension(2:3), parameter :: w = [3.0_pReal,2.0_pReal] + do d = lbound(FEM_quadrature_weights,1), ubound(FEM_quadrature_weights,1) do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1) if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) & From 58a79219666525ae8be6d86051e190da983f3573 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 22:45:04 +0200 Subject: [PATCH 16/26] no need for temp, pack instead of reshape for 1D array --- src/mesh/FEM_quadrature.f90 | 68 ++++++++++--------------------------- 1 file changed, 18 insertions(+), 50 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index 278bfb5db..bf1f1567f 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -20,12 +20,12 @@ module FEM_quadrature -1.0_pReal, 1.0_pReal, -1.0_pReal, & -1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4]) - type :: group_float !< variable length datatype used for storage of state + type :: group_float !< variable length datatype real(pReal), dimension(:), allocatable :: p end type group_float integer, dimension(2:3,maxOrder), public, protected :: & - FEM_nQuadrature !< number of quadrature points for a given spatial dimension(2-3) and interpolation order(1-maxOrder) + FEM_nQuadrature !< number of quadrature points for spatial dimension(2-3) and interpolation order (1-maxOrder) type(group_float), dimension(2:3,maxOrder), public, protected :: & FEM_quadrature_weights, & !< quadrature weights for each quadrature rule FEM_quadrature_points !< quadrature point coordinates (in simplical system) for each quadrature rule @@ -191,13 +191,9 @@ pure function permutationStar3(point) result(qPt) real(pReal), dimension(2) :: qPt real(pReal), dimension(1), intent(in) :: point - real(pReal), dimension(3,1) :: temp - - temp = reshape([ & - point(1), point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(triangle, temp),[2]) + qPt = pack(matmul(triangle,reshape([ & + point(1), point(1), point(1)],[3,1])),.true.) end function permutationStar3 @@ -210,15 +206,11 @@ pure function permutationStar21(point) result(qPt) real(pReal), dimension(6) :: qPt real(pReal), dimension(1), intent(in) :: point - real(pReal), dimension(3,3) :: temp - - temp = reshape([ & + qPt = pack(matmul(triangle,reshape([ & point(1), point(1), 1.0_pReal - 2.0_pReal*point(1), & point(1), 1.0_pReal - 2.0_pReal*point(1), point(1), & - 1.0_pReal - 2.0_pReal*point(1), point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(triangle, temp),[6]) + 1.0_pReal - 2.0_pReal*point(1), point(1), point(1)],[3,3])),.true.) end function permutationStar21 @@ -231,18 +223,14 @@ pure function permutationStar111(point) result(qPt) real(pReal), dimension(12) :: qPt real(pReal), dimension(2), intent(in) :: point - real(pReal), dimension(3,6) :: temp - - temp = reshape([ & + qPt = pack(matmul(triangle,reshape([ & point(1), point(2), 1.0_pReal - point(1) - point(2), & point(1), 1.0_pReal - point(1) - point(2), point(2), & point(2), point(1), 1.0_pReal - point(1) - point(2), & point(2), 1.0_pReal - point(1) - point(2), point(1), & 1.0_pReal - point(1) - point(2), point(2), point(1), & - 1.0_pReal - point(1) - point(2), point(1), point(2)],shape(temp)) - - qPt = reshape(matmul(triangle, temp),[12]) + 1.0_pReal - point(1) - point(2), point(1), point(2)],[3,6])),.true.) end function permutationStar111 @@ -255,13 +243,9 @@ pure function permutationStar4(point) result(qPt) real(pReal), dimension(3) :: qPt real(pReal), dimension(1), intent(in) :: point - real(pReal), dimension(4,1) :: temp - - temp = reshape([ & - point(1), point(1), point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(tetrahedron, temp),[3]) + qPt = pack(matmul(tetrahedron,reshape([ & + point(1), point(1), point(1), point(1)],[4,1])),.true.) end function permutationStar4 @@ -274,16 +258,12 @@ pure function permutationStar31(point) result(qPt) real(pReal), dimension(12) :: qPt real(pReal), dimension(1), intent(in) :: point - real(pReal), dimension(4,4) :: temp - - temp = reshape([ & + qPt = pack(matmul(tetrahedron,reshape([ & point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), & point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), & point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), & - 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(tetrahedron, temp),[12]) + 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)],[4,4])),.true.) end function permutationStar31 @@ -296,18 +276,14 @@ function permutationStar22(point) result(qPt) real(pReal), dimension(18) :: qPt real(pReal), dimension(1), intent(in) :: point - real(pReal), dimension(4,6) :: temp - - temp = reshape([ & + qPt = pack(matmul(tetrahedron,reshape([ & point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), & point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), & 0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1), & 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1), & 0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1), & - point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(tetrahedron, temp),[18]) + point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)],[4,6])),.true.) end function permutationStar22 @@ -320,10 +296,8 @@ pure function permutationStar211(point) result(qPt) real(pReal), dimension(36) :: qPt real(pReal), dimension(2), intent(in) :: point - real(pReal), dimension(4,12) :: temp - - temp = reshape([ & + qPt = pack(matmul(tetrahedron,reshape([ & point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), & point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), & point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), & @@ -335,9 +309,7 @@ pure function permutationStar211(point) result(qPt) point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), & 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2), & 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1), & - 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)],shape(temp)) - - qPt = reshape(matmul(tetrahedron, temp),[36]) + 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)],[4,12])),.true.) end function permutationStar211 @@ -350,10 +322,8 @@ pure function permutationStar1111(point) result(qPt) real(pReal), dimension(72) :: qPt real(pReal), dimension(3), intent(in) :: point - real(pReal), dimension(4,24) :: temp - - temp = reshape([ & + qPt = pack(matmul(tetrahedron,reshape([ & point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), & point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), & point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), & @@ -377,9 +347,7 @@ pure function permutationStar1111(point) result(qPt) 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3), & 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1), & 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2), & - 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)],shape(temp)) - - qPt = reshape(matmul(tetrahedron, temp),[72]) + 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)],[4,24])),.true.) end function permutationStar1111 From de7ef430958183843867ea2b42430f16919af16f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 23:43:55 +0200 Subject: [PATCH 17/26] documenting/silencing gfortran --- src/mesh/FEM_quadrature.f90 | 6 +++--- src/mesh/discretization_mesh.f90 | 6 ++---- src/quit.f90 | 5 ++++- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/mesh/FEM_quadrature.f90 b/src/mesh/FEM_quadrature.f90 index bf1f1567f..0c470ca6b 100644 --- a/src/mesh/FEM_quadrature.f90 +++ b/src/mesh/FEM_quadrature.f90 @@ -124,7 +124,7 @@ subroutine FEM_quadrature_init() allocate(FEM_quadrature_weights(3,2)%p(FEM_nQuadrature(3,2))) FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal - FEM_quadrature_points (3,2)%p = permutationStar31([0.1381966011250105_pReal]) + FEM_quadrature_points (3,2)%p = permutationStar31([1.3819660112501052e-1_pReal]) !-------------------------------------------------------------------------------------------------- ! 3D cubic @@ -141,7 +141,7 @@ subroutine FEM_quadrature_init() permutationStar22([4.5503704125649649e-2_pReal]) ] !-------------------------------------------------------------------------------------------------- -! 3D quartic +! 3D quartic (lower precision/unknown source) FEM_nQuadrature(3,4) = 35 allocate(FEM_quadrature_weights(3,4)%p(FEM_nQuadrature(3,4))) @@ -159,7 +159,7 @@ subroutine FEM_quadrature_init() permutationStar4([0.25_pReal]) ] !-------------------------------------------------------------------------------------------------- -! 3D quintic +! 3D quintic (lower precision/unknown source) FEM_nQuadrature(3,5) = 56 allocate(FEM_quadrature_weights(3,5)%p(FEM_nQuadrature(3,5))) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 49eaf0e18..422f4ca85 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -71,16 +71,14 @@ subroutine discretization_mesh_init(restart) logical, intent(in) :: restart - integer, allocatable, dimension(:) :: chunkPos integer :: dimPlex, & mesh_Nnodes, & !< total number of nodes in mesh - j, l, & + j, & debug_element, debug_ip PetscSF :: sf DM :: globalMesh PetscInt :: nFaceSets PetscInt, pointer, dimension(:) :: pFaceSets - character(len=pStringLen), dimension(:), allocatable :: fileContent IS :: faceSetIS PetscErrorCode :: ierr integer, dimension(:), allocatable :: & @@ -88,7 +86,7 @@ subroutine discretization_mesh_init(restart) class(tNode), pointer :: & num_mesh integer :: integrationOrder !< order of quadrature rule required - type(tvec) :: coords_node0 + type(tvec) :: coords_node0 print'(/,a)', ' <<<+- discretization_mesh init -+>>>' diff --git a/src/quit.f90 b/src/quit.f90 index 786241d85..6361ff783 100644 --- a/src/quit.f90 +++ b/src/quit.f90 @@ -6,7 +6,10 @@ !-------------------------------------------------------------------------------------------------- subroutine quit(stop_id) #include - use PetscSys + use PETScSys +#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) + use MPI_f08 +#endif use HDF5 implicit none From 9767866c1e8d944483fb40857458dda05cb55fab Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 07:24:36 +0200 Subject: [PATCH 18/26] cleaned+updated test --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 567067daf..96cf7af03 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 567067daf05866c8e022fc0af6f441840f143b81 +Subproject commit 96cf7af039b3153f91c037bcce62eb456af51f4b From 173a5f8e55e365a6f875096686953b6d8e1db717 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 07:48:01 +0200 Subject: [PATCH 19/26] less public/unneeded variables --- src/mesh/DAMASK_mesh.f90 | 1 + src/mesh/FEM_utilities.f90 | 17 ++++------------- src/mesh/discretization_mesh.f90 | 10 +++++----- 3 files changed, 10 insertions(+), 18 deletions(-) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 342dec5c0..2eb6a7d91 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -21,6 +21,7 @@ program DAMASK_mesh use mesh_mechanical_FEM implicit none + integer :: nActiveFields = 0 !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 5bb6d5f44..1c073dd69 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -23,14 +23,8 @@ module FEM_utilities implicit none private -!-------------------------------------------------------------------------------------------------- - logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill - integer, public, parameter :: maxFields = 6 - integer, public :: nActiveFields = 0 - -!-------------------------------------------------------------------------------------------------- -! grid related information information - real(pReal), public :: wgt !< weighting factor 1/Nelems + logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill + real(pReal), public, protected :: wgt !< weighting factor 1/Nelems !-------------------------------------------------------------------------------------------------- @@ -49,10 +43,6 @@ module FEM_utilities COMPONENT_MECH_Z_ID end enum -!-------------------------------------------------------------------------------------------------- -! variables controlling debugging - logical :: & - debugPETSc !< use some in debug defined options for more verbose PETSc solution !-------------------------------------------------------------------------------------------------- ! derived types @@ -109,8 +99,9 @@ subroutine FEM_utilities_init integer :: structOrder !< order of displacement shape functions character(len=*), parameter :: & PETSCDEBUG = ' -snes_view -snes_monitor ' - PetscErrorCode :: ierr + logical :: debugPETSc !< use some in debug defined options for more verbose PETSc solution + print'(/,a)', ' <<<+- FEM_utilities init -+>>>' diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 422f4ca85..e5f989484 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -40,6 +40,11 @@ module discretization_mesh mesh_maxNips !< max number of IPs in any CP element !!!! BEGIN DEPRECATED !!!!! + DM, public :: geomMesh + + PetscInt, dimension(:), allocatable, public, protected :: & + mesh_boundaries + real(pReal), dimension(:,:), allocatable :: & mesh_ipVolume, & !< volume associated with IP (initially!) mesh_node0 !< node x,y,z coordinates (initially!) @@ -50,11 +55,6 @@ module discretization_mesh real(pReal), dimension(:,:,:), allocatable :: & mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) - DM, public :: geomMesh - - PetscInt, dimension(:), allocatable, public, protected :: & - mesh_boundaries - public :: & discretization_mesh_init, & mesh_FEM_build_ipVolumes, & From eb834b635dee5e301dcea4d0d8fef29e01a38b79 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 07:58:34 +0200 Subject: [PATCH 20/26] unify notation with grid consistent style: Symbols, not descriptions. also removed untested/unused loginc functionality. Once load cases are written in YAML, we can introduce the scaling as in DAMASK_grid --- PRIVATE | 2 +- examples/mesh/tensionZ.load | 4 +-- src/mesh/DAMASK_mesh.f90 | 70 +++++++++++++++---------------------- src/mesh/FEM_utilities.f90 | 10 ------ 4 files changed, 32 insertions(+), 54 deletions(-) diff --git a/PRIVATE b/PRIVATE index 96cf7af03..d92b030e5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 96cf7af039b3153f91c037bcce62eb456af51f4b +Subproject commit d92b030e5777e718c77edc2e1e93abfa0981b024 diff --git a/examples/mesh/tensionZ.load b/examples/mesh/tensionZ.load index eb9a7c426..b3588daef 100644 --- a/examples/mesh/tensionZ.load +++ b/examples/mesh/tensionZ.load @@ -1,11 +1,11 @@ #initial elastic step -$Loadcase 1 time 0.0005 incs 1 frequency 5 +$Loadcase 1 t 0.0005 N 1 f_out 5 Face 1 X 0.01 Face 2 X 0.0 Face 2 Y 0.0 Face 2 Z 0.0 $EndLoadcase -$Loadcase 2 time 10.0 incs 200 frequency 5 +$Loadcase 2 t 10.0 N 200 f_out 5 Face 1 X 0.01 Face 2 X 0.0 Face 2 Y 0.0 diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 2eb6a7d91..d1568e87f 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -23,6 +23,15 @@ program DAMASK_mesh implicit none integer :: nActiveFields = 0 + type :: tLoadCase + real(pReal) :: time = 0.0_pReal !< length of increment + integer :: incs = 0, & !< number of increments + outputfrequency = 1 !< frequency of result writes + logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase + integer, allocatable, dimension(:) :: faceID + type(tFieldBC), allocatable, dimension(:) :: fieldBC + end type tLoadCase + !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file integer, allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing @@ -104,8 +113,8 @@ program DAMASK_mesh chunkPos = IO_stringPos(line) do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase - select case (IO_lc(IO_stringValue(line,chunkPos,i))) - case('$loadcase') + select case (IO_stringValue(line,chunkPos,i)) + case('$Loadcase') N_def = N_def + 1 end select enddo ! count all identifiers to allocate memory and do sanity check @@ -152,39 +161,36 @@ program DAMASK_mesh chunkPos = IO_stringPos(line) do i = 1, chunkPos(1) - select case (IO_lc(IO_stringValue(line,chunkPos,i))) + select case (IO_stringValue(line,chunkPos,i)) !-------------------------------------------------------------------------------------------------- ! loadcase information - case('$loadcase') + case('$Loadcase') currentLoadCase = IO_intValue(line,chunkPos,i+1) - case('face') + case('Face') currentFace = IO_intValue(line,chunkPos,i+1) currentFaceSet = -1 do faceSet = 1, mesh_Nboundaries if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet enddo if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') - case('t','time','delta') ! increment time + case('t') loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1) - case('n','incs','increments','steps') ! number of increments + case('N') loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1) - case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling) - loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1) - loadCases(currentLoadCase)%logscale = 1 - case('freq','frequency','outputfreq') ! frequency of result writings + case('f_out') loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1) - case('guessreset','dropguessing') + case('estimate_rate') loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory !-------------------------------------------------------------------------------------------------- ! boundary condition information - case('x','y','z') - select case(IO_lc(IO_stringValue(line,chunkPos,i))) - case('x') + case('X','Y','Z') + select case(IO_stringValue(line,chunkPos,i)) + case('X') ID = COMPONENT_MECH_X_ID - case('y') + case('Y') ID = COMPONENT_MECH_Y_ID - case('z') + case('Z') ID = COMPONENT_MECH_Z_ID end select @@ -223,8 +229,8 @@ program DAMASK_mesh do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) & print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), & - ' Component ', component, & - ' Value ', loadCases(currentLoadCase)%fieldBC(field)% & + ' Component ', component, & + ' Value ', loadCases(currentLoadCase)%fieldBC(field)% & componentBC(component)%Value(faceSet) enddo enddo @@ -267,32 +273,14 @@ program DAMASK_mesh !-------------------------------------------------------------------------------------------------- ! forwarding time timeIncOld = timeinc ! last timeinc that brought former inc to an end - if (loadCases(currentLoadCase)%logscale == 0) then ! linear scale - timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) - else - if (currentLoadCase == 1) then ! 1st load case of logarithmic scale - if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd - else ! not-1st inc of 1st load case of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal)) - endif - else ! not-1st load case of logarithmic scale - timeinc = time0 * & - ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/& - real(loadCases(currentLoadCase)%incs ,pReal))& - -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1 ,pReal)/& - real(loadCases(currentLoadCase)%incs ,pReal))) - endif - endif + timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step - - - stepFraction = 0 ! fraction scaled by stepFactor**cutLevel + stepFraction = 0 ! fraction scaled by stepFactor**cutLevel subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time - time = time + timeinc ! forward target time - stepFraction = stepFraction + 1 ! count step + time = time + timeinc ! forward target time + stepFraction = stepFraction + 1 ! count step !-------------------------------------------------------------------------------------------------- ! report begin of new step diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 1c073dd69..319997ad6 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -64,16 +64,6 @@ module FEM_utilities type(tComponentBC), allocatable :: componentBC(:) end type tFieldBC - type, public :: tLoadCase - real(pReal) :: time = 0.0_pReal !< length of increment - integer :: incs = 0, & !< number of increments - outputfrequency = 1, & !< frequency of result writes - logscale = 0 !< linear/logarithmic time inc flag - logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase - integer, allocatable, dimension(:) :: faceID - type(tFieldBC), allocatable, dimension(:) :: fieldBC - end type tLoadCase - public :: & FEM_utilities_init, & utilities_constitutiveResponse, & From 5b66db8a394ffe34b303039fe8cd17b370e60a30 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 10:12:43 +0200 Subject: [PATCH 21/26] only mechanics at the moment will be extended, but most likely differently --- src/mesh/DAMASK_mesh.f90 | 113 ++++++++++++------------------------- src/mesh/FEM_utilities.f90 | 8 +-- 2 files changed, 41 insertions(+), 80 deletions(-) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index d1568e87f..a8446d9fa 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -21,7 +21,6 @@ program DAMASK_mesh use mesh_mechanical_FEM implicit none - integer :: nActiveFields = 0 type :: tLoadCase real(pReal) :: time = 0.0_pReal !< length of increment @@ -78,7 +77,7 @@ program DAMASK_mesh type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tSolutionState), allocatable, dimension(:) :: solres - PetscInt :: faceSet, currentFaceSet, field, dimPlex + PetscInt :: faceSet, currentFaceSet, dimPlex PetscErrorCode :: ierr integer(kind(COMPONENT_UNDEFINED_ID)) :: ID external :: & @@ -101,8 +100,7 @@ program DAMASK_mesh ! reading basic information from load case file and allocate data structure containing load cases call DMGetDimension(geomMesh,dimPlex,ierr) !< dimension of mesh (2D or 3D) CHKERRA(ierr) - nActiveFields = 1 - allocate(solres(nActiveFields)) + allocate(solres(1)) !-------------------------------------------------------------------------------------------------- ! reading basic information from load case file and allocate data structure containing load cases @@ -124,32 +122,26 @@ program DAMASK_mesh allocate(loadCases(N_def)) do i = 1, size(loadCases) - allocate(loadCases(i)%fieldBC(nActiveFields)) - field = 1 - loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID + allocate(loadCases(i)%fieldBC(1)) + loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID enddo do i = 1, size(loadCases) - do field = 1, nActiveFields - select case (loadCases(i)%fieldBC(field)%ID) - case(FIELD_MECH_ID) - loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements - allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents)) - do component = 1, loadCases(i)%fieldBC(field)%nComponents - select case (component) - case (1) - loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_X_ID - case (2) - loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Y_ID - case (3) - loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Z_ID - end select - enddo + loadCases(i)%fieldBC(1)%nComponents = dimPlex !< X, Y (, Z) displacements + allocate(loadCases(i)%fieldBC(1)%componentBC(loadCases(i)%fieldBC(1)%nComponents)) + do component = 1, loadCases(i)%fieldBC(1)%nComponents + select case (component) + case (1) + loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_X_ID + case (2) + loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Y_ID + case (3) + loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID end select - do component = 1, loadCases(i)%fieldBC(field)%nComponents - allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal) - allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) - enddo + enddo + do component = 1, loadCases(i)%fieldBC(1)%nComponents + allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal) + allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) enddo enddo @@ -194,16 +186,12 @@ program DAMASK_mesh ID = COMPONENT_MECH_Z_ID end select - do field = 1, nActiveFields - if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then - do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents - if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == ID) then - loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = & - .true. - loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = & - IO_floatValue(line,chunkPos,i+1) - endif - enddo + do component = 1, loadcases(currentLoadCase)%fieldBC(1)%nComponents + if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%ID == ID) then + loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask (currentFaceSet) = & + .true. + loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = & + IO_floatValue(line,chunkPos,i+1) endif enddo end select @@ -219,21 +207,16 @@ program DAMASK_mesh print'(a,i0)', ' load case: ', currentLoadCase if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & print'(a)', ' drop guessing along trajectory' - do field = 1, nActiveFields - select case (loadCases(currentLoadCase)%fieldBC(field)%ID) - case(FIELD_MECH_ID) - print'(a)', ' Field '//trim(FIELD_MECH_label) + print'(a)', ' Field '//trim(FIELD_MECH_label) - end select - do faceSet = 1, mesh_Nboundaries - do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents - if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) & - print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), & - ' Component ', component, & - ' Value ', loadCases(currentLoadCase)%fieldBC(field)% & - componentBC(component)%Value(faceSet) - enddo - enddo + do faceSet = 1, mesh_Nboundaries + do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents + if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) & + print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), & + ' Component ', component, & + ' Value ', loadCases(currentLoadCase)%fieldBC(1)% & + componentBC(component)%Value(faceSet) + enddo enddo print'(a,f12.6)', ' time: ', loadCases(currentLoadCase)%time if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count @@ -247,12 +230,7 @@ program DAMASK_mesh !-------------------------------------------------------------------------------------------------- ! doing initialization depending on active solvers call FEM_Utilities_init - do field = 1, nActiveFields - select case (loadCases(1)%fieldBC(field)%ID) - case(FIELD_MECH_ID) - call FEM_mechanical_init(loadCases(1)%fieldBC(field)) - end select - enddo + call FEM_mechanical_init(loadCases(1)%fieldBC(1)) if (worldrank == 0) then open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') @@ -295,33 +273,16 @@ program DAMASK_mesh '-',stepFraction, '/', subStepFactor**cutBackLevel flush(IO_STDOUT) -!-------------------------------------------------------------------------------------------------- -! forward fields - do field = 1, nActiveFields - select case (loadCases(currentLoadCase)%fieldBC(field)%ID) - case(FIELD_MECH_ID) - call FEM_mechanical_forward (& - guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) - - end select - enddo + call FEM_mechanical_forward(guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) !-------------------------------------------------------------------------------------------------- ! solve fields stagIter = 0 stagIterate = .true. do while (stagIterate) - do field = 1, nActiveFields - select case (loadCases(currentLoadCase)%fieldBC(field)%ID) - case(FIELD_MECH_ID) - solres(field) = FEM_mechanical_solution (& - incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) + solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) + if(.not. solres(1)%converged) exit - end select - - if(.not. solres(field)%converged) exit ! no solution found - - enddo stagIter = stagIter + 1 stagIterate = stagIter < stagItMax & .and. all(solres(:)%converged) & diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 319997ad6..64bcc3896 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -53,15 +53,15 @@ module FEM_utilities end type tSolutionState type, public :: tComponentBC - integer(kind(COMPONENT_UNDEFINED_ID)) :: ID - real(pReal), allocatable, dimension(:) :: Value - logical, allocatable, dimension(:) :: Mask + integer(kind(COMPONENT_UNDEFINED_ID)) :: ID + real(pReal), allocatable, dimension(:) :: Value + logical, allocatable, dimension(:) :: Mask end type tComponentBC type, public :: tFieldBC integer(kind(FIELD_UNDEFINED_ID)) :: ID integer :: nComponents = 0 - type(tComponentBC), allocatable :: componentBC(:) + type(tComponentBC), allocatable, dimension(:) :: componentBC end type tFieldBC public :: & From aa03e5930632147cf984710d0fef6f41c68eb37d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 11:22:25 +0200 Subject: [PATCH 22/26] same statement, different syntax but silences gfortran runtime warning with automated LHS, compiling in DEBUG mode (tested for mesh solver) gives: home/m/DAMASK/src/phase_mechanical.f90:1010:43: runtime error: signed integer overflow: -9223372036854775808 - 1 cannot be represented in type 'integer(kind=8)' --- src/phase_mechanical.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 8b04a0009..3c14accb5 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1007,7 +1007,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) - subState0 = plasticState(ph)%State0(:,en) + allocate(subState0,source=plasticState(ph)%State0(:,en)) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) From 69cc0b528bf8a02695a59ed315cc316c852a07fb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 11:38:03 +0200 Subject: [PATCH 23/26] do calculation for initialization --- src/mesh/mesh_mech_FEM.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 1468e3fcc..6bcf44c3a 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -109,7 +109,7 @@ subroutine FEM_mechanical_init(fieldBC) character(len=*), parameter :: prefix = 'mechFE_' PetscErrorCode :: ierr - + real(pReal), dimension(3,3) :: devNull class(tNode), pointer :: & num_mesh @@ -258,6 +258,7 @@ subroutine FEM_mechanical_init(fieldBC) call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr) CHKERRQ(ierr) enddo + call utilities_constitutiveResponse(0.0_pReal,devNull,.true.) end subroutine FEM_mechanical_init @@ -397,7 +398,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(params%timeinc,P_av,ForwardData) + call utilities_constitutiveResponse(params%timeinc,P_av,ForwardData) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) ForwardData = .false. From 5ac592eb9e331de2988f97946237092882ad5795 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 12:31:10 +0200 Subject: [PATCH 24/26] initialize L_(i/p), L_(i,p)0 --- src/mesh/mesh_mech_FEM.f90 | 4 ++-- src/phase_mechanical.f90 | 28 ++++++++++++++-------------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 6bcf44c3a..d6d314a42 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -289,8 +289,8 @@ type(tSolutionState) function FEM_mechanical_solution( & params%timeinc = timeinc params%fieldBC = fieldBC - call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution) - call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? + call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution) + call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? terminallyIll = .false. if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 3c14accb5..a14102513 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -215,8 +215,8 @@ module subroutine mechanical_init(materials,phases) allocate(phase_mechanical_F0(phases%length)) allocate(phase_mechanical_Li(phases%length)) allocate(phase_mechanical_Li0(phases%length)) - allocate(phase_mechanical_Lp0(phases%length)) allocate(phase_mechanical_Lp(phases%length)) + allocate(phase_mechanical_Lp0(phases%length)) allocate(phase_mechanical_S(phases%length)) allocate(phase_mechanical_P(phases%length)) allocate(phase_mechanical_S0(phases%length)) @@ -224,20 +224,20 @@ module subroutine mechanical_init(materials,phases) do ph = 1, phases%length Nmembers = count(material_phaseID == ph) - allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_F(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers)) + allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers),source=0.0_pReal) + allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers),source=0.0_pReal) + allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers),source=0.0_pReal) + allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal) - allocate(phase_mechanical_F(ph)%data(3,3,Nmembers)) - allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers)) phase => phases%get(ph) mech => phase%get('mechanical') @@ -494,7 +494,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) enddo LpLoop call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & - S, Fi_new, ph,en) + S, Fi_new, ph,en) !* update current residuum and check for convergence of loop atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error @@ -1130,12 +1130,12 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF) en = material_phaseEntry(co,ce) call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & - phase_mechanical_Fe(ph)%data(1:3,1:3,en), & - phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en) + phase_mechanical_Fe(ph)%data(1:3,1:3,en), & + phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en) call phase_LiAndItsTangents(devNull,dLidS,dLidFi, & - phase_mechanical_S(ph)%data(1:3,1:3,en), & - phase_mechanical_Fi(ph)%data(1:3,1:3,en), & - ph,en) + phase_mechanical_S(ph)%data(1:3,1:3,en), & + phase_mechanical_Fi(ph)%data(1:3,1:3,en), & + ph,en) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)) invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,en)) From 11cdc63a2768013ec5b5e7e283431f65b6f5cab1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 22 Jul 2021 13:50:41 +0200 Subject: [PATCH 25/26] use reasonable load case, even for compile test --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index d92b030e5..f2e2d6c71 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit d92b030e5777e718c77edc2e1e93abfa0981b024 +Subproject commit f2e2d6c71ea798bfc63230a756b7cf9748599bec From 93814054d51f6ccc465f51208a874ca2216e79e9 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 26 Jul 2021 15:02:14 +0200 Subject: [PATCH 26/26] [skip ci] updated version information after successful test of v3.0.0-alpha4-171-g6a9c892d7 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 026bbcf9f..c4c59cb81 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha4-137-gb69b85754 +v3.0.0-alpha4-171-g6a9c892d7