Skip to content

Commit bd0f307

Browse files
MarDiehlMarDiehl
committed
Merge branch 'implement_node_BC_mesh' into 'development'
Node BC for the mesh solver See merge request damask/DAMASK!1067
2 parents 316609c + 7c61502 commit bd0f307

File tree

4 files changed

+139
-97
lines changed

4 files changed

+139
-97
lines changed

‎src/mesh/DAMASK_mesh.f90‎

Lines changed: 46 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -68,16 +68,22 @@ program DAMASK_mesh
6868
loadcase_string
6969
integer :: &
7070
stagItMax, & !< max number of field level staggered iterations
71-
maxCutBack !< max number of cutbacks
71+
maxCutBack, & !< max number of cutbacks
72+
skipBCTagDigits, & !< number of characters to skip to print BC tag (T descriptor)
73+
BCTag !< tag value read from YAML load file
74+
integer, allocatable, dimension(:) :: &
75+
knownBCTags !< array of read BC tags (check for duplicates)
76+
7277

7378
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
7479
type(tSolutionState), allocatable, dimension(:) :: solres
75-
PetscInt :: faceSet, currentFaceSet, dimPlex
80+
PetscInt :: boundary, dimPlex
7681
PetscErrorCode :: err_PETSc
7782
external :: &
7883
quit
7984
character(len=:), allocatable :: &
80-
fileContent, fname
85+
fileContent, fname, tagPrintFormat
86+
character(len=6) :: BC_elem
8187

8288

8389
!--------------------------------------------------------------------------------------------------
@@ -114,34 +120,40 @@ program DAMASK_mesh
114120
load_steps => load%get_list('loadstep')
115121

116122
allocate(loadCases(size(load_steps)))
117-
118-
119123
do l = 1, size(load_steps)
120124
load_step => load_steps%get_dict(l)
121125
step_bc => load_step%get_dict('boundary_conditions')
122126
step_mech => step_bc%get_list('mechanical')
123127
allocate(loadCases(l)%mechBC(mesh_Nboundaries))
124128
loadCases(l)%mechBC(:)%nComponents = dimPlex !< X, Y (, Z) displacements
125-
do faceSet = 1, mesh_Nboundaries
126-
allocate(loadCases(l)%mechBC(faceSet)%Value(dimPlex), source = 0.0_pREAL)
127-
allocate(loadCases(l)%mechBC(faceSet)%Mask(dimPlex), source = .false.)
129+
do boundary = 1_pPETSCINT, mesh_Nboundaries
130+
allocate(loadCases(l)%mechBC(boundary)%dot_u(dimPlex), source = 0.0_pREAL)
131+
allocate(loadCases(l)%mechBC(boundary)%active(dimPlex), source = .false.)
128132
end do
129133

134+
allocate(knownBCTags(size(step_mech)), source = -1)
130135
do m = 1, size(step_mech)
131136
mech_BC => step_mech%get_dict(m)
132-
currentFaceSet = -1
133-
do faceSet = 1, mesh_Nboundaries
134-
if (mesh_boundaries(faceSet) == mech_BC%get_asInt('tag')) currentFaceSet = faceSet
135-
end do
136-
if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
137+
BCTag = mech_BC%get_asInt('tag')
138+
boundary = findloc(mesh_boundariesIS, BCtag, dim = 1)
139+
if (boundary == 0) then !< tag not defined in mesh file
140+
call IO_error(error_ID = 837, ext_msg = 'BC tag '//mech_BC%get_asStr('tag')// &
141+
' refers to inexistent node/edge/face mesh group')
142+
else if (findloc(knownBCTags, BCTag, dim = 1) /= 0) then !< duplicated tag
143+
call IO_error(error_ID = 837, ext_msg = 'BC redefinition not allowed (tag '// &
144+
mech_BC%get_asStr('tag')//')')
145+
end if
146+
knownBCTags(m) = BCTag
137147
mech_u => mech_BC%get_list('dot_u')
138148
do component = 1, dimPlex
139149
if (mech_u%get_asStr(component) /= 'x') then
140-
loadCases(l)%mechBC(currentFaceSet)%Mask(component) = .true.
141-
loadCases(l)%mechBC(currentFaceSet)%Value(component) = mech_u%get_asReal(component)
150+
loadCases(l)%mechBC(boundary)%active(component) = .true.
151+
loadCases(l)%mechBC(boundary)%dot_u(component) = mech_u%get_asReal(component)
142152
end if
143153
end do
144154
end do
155+
deallocate(knownBCTags)
156+
145157
step_discretization => load_step%get_dict('discretization')
146158
loadCases(l)%t = step_discretization%get_asReal('t')
147159
loadCases(l)%N = step_discretization%get_asInt ('N')
@@ -157,28 +169,33 @@ program DAMASK_mesh
157169
!--------------------------------------------------------------------------------------------------
158170
! consistency checks and output of load case
159171
errorID = 0
172+
skipBCTagDigits = 4+6+1+floor(log10(real(maxval(mesh_boundariesIS))))+1+2 !< Indentation(4)+BC_elem(6)+blank(1)+descriptor T(1)+NumDigits(floor(..)+1)+2(start printing after 1 blank)
173+
allocate(character(len=skipBCTagDigits) :: tagPrintFormat)
174+
write(tagPrintFormat,'(i0)') skipBCTagDigits
160175
checkLoadcases: do l = 1, size(load_steps)
161176
write (loadcase_string, '(i0)' ) l
162177
print'(/,1x,a,1x,i0)', 'load case:', l
163178
if (.not. loadCases(l)%estimate_rate) &
164179
print'(2x,a)', 'drop guessing along trajectory'
165180
print'(2x,a)', 'Field '//trim(FIELD_MECH_label)
166181

167-
do faceSet = 1, mesh_Nboundaries
168-
do component = 1, dimPlex
169-
if (loadCases(l)%mechBC(faceSet)%Mask(component)) &
170-
print'(a,i2,a,i2,a,f12.7)', &
171-
' Face ', mesh_boundaries(faceSet), &
172-
' Component ', component, &
173-
' Value ', loadCases(l)%mechBC(faceSet)%Value(component)
182+
do boundary = 1_pPETSCINT, mesh_Nboundaries
183+
BC_elem = merge('Vertex', merge('Edge ', 'Face ', dimplex == 2_pPETSCINT), &
184+
mesh_boundariesIdx(boundary) == 1_pPETSCINT)
185+
do component = 1_pPETSCINT, dimPlex
186+
if (loadCases(l)%mechBC(boundary)%active(component)) &
187+
print'(4x,a,1x,i0,T'//tagPrintFormat//',a,1x,i1,1x,a,1x,f12.7)', &
188+
BC_elem, mesh_boundariesIS(boundary), &
189+
'Component', component, &
190+
'Value', loadCases(l)%mechBC(boundary)%dot_u(component)
174191
end do
175192
end do
176-
print'(2x,a,f12.6)', 'time: ', loadCases(l)%t
177-
if (loadCases(l)%N < 1) errorID = 835 ! non-positive incs count
178-
print'(2x,a,i5)', 'increments: ', loadCases(l)%N
179-
if (loadCases(l)%f_out < 1) errorID = 836 ! non-positive result frequency
180-
print'(2x,a,i5)', 'output frequency: ', &
181-
loadCases(l)%f_out
193+
194+
print'(2x,a,T21,g0.6)', 'time:', loadCases(l)%t
195+
if (loadCases(l)%N < 1) errorID = 835 ! non-positive incs count
196+
print'(2x,a,T21,i0)', 'increments:', loadCases(l)%N
197+
if (loadCases(l)%f_out < 1) errorID = 836 ! non-positive result frequency
198+
print'(2x,a,T21,i0)', 'output frequency:', loadCases(l)%f_out
182199
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
183200
end do checkLoadcases
184201

@@ -222,7 +239,7 @@ program DAMASK_mesh
222239
'Time', t, &
223240
's: Increment ', inc, '/', loadCases(l)%N,&
224241
'-', stepFraction, '/', subStepFactor**cutBackLevel,&
225-
' of load case ', l,'/',size(load_steps)
242+
' of load case ', l,'/', size(load_steps)
226243
write(incInfo,'(4(a,i0))') &
227244
'Increment ',totalIncsCounter,'/',sum(loadCases%N),&
228245
'-',stepFraction, '/', subStepFactor**cutBackLevel

‎src/mesh/FEM_utilities.f90‎

Lines changed: 14 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -48,8 +48,8 @@ module FEM_utilities
4848

4949
type, public :: tMechBC
5050
integer :: nComponents = 0
51-
real(pREAL), allocatable, dimension(:) :: Value
52-
logical, allocatable, dimension(:) :: Mask
51+
real(pREAL), allocatable, dimension(:) :: dot_u
52+
logical, allocatable, dimension(:) :: active
5353
end type tMechBC
5454

5555
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR<23)
@@ -161,35 +161,37 @@ subroutine utilities_projectBCValues(dm_local,solution_local,section,mechBC,Delt
161161
real(pREAL), intent(in) :: Delta_t
162162
PetscInt, intent(in) :: dimPlex
163163

164-
PetscInt :: component, face, bcSize, &
164+
PetscInt :: component, boundary, bcSize, &
165165
nBcPoints, point, dof, numDof, numComp, offset
166166
IS :: bcPointsIS
167167
PetscErrorCode :: err_PETSc
168168
PetscInt, pointer :: bcPoints(:)
169169
real(pREAL), pointer :: localArray(:)
170170

171+
character(len=11) :: setLabel
171172

172173

173-
do face = 1, mesh_Nboundaries; do component = 1, dimPlex
174-
if (mechBC(face)%Mask(component)) then
175-
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
176-
if (bcSize > 0) then
177-
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPointsIS,err_PETSc)
174+
do boundary = 1_pPETSCINT, mesh_Nboundaries; do component = 1_pPETSCINT, dimPlex
175+
if (mechBC(boundary)%active(component)) then
176+
setLabel = mesh_BCTypeLabel(mesh_boundariesIdx(boundary))
177+
call DMGetStratumSize(dm_local,setLabel,mesh_boundariesIS(boundary),bcSize,err_PETSc)
178+
if (bcSize > 0_pPETSCINT) then
179+
call DMGetStratumIS(dm_local,setLabel,mesh_boundariesIS(boundary),bcPointsIS,err_PETSc)
178180
CHKERRQ(err_PETSc)
179181
call PetscSectionGetFieldComponents(section,0_pPETSCINT,numComp,err_PETSc)
180182
CHKERRQ(err_PETSc)
181183
call ISGetSize(bcPointsIS,nBcPoints,err_PETSc)
182184
CHKERRQ(err_PETSc)
183-
if (nBcPoints > 0) call ISGetIndices(bcPointsIS,bcPoints,err_PETSc)
185+
if (nBcPoints > 0_pPETSCINT) call ISGetIndices(bcPointsIS,bcPoints,err_PETSc)
184186
call VecGetArray(solution_local,localArray,err_PETSc)
185187
CHKERRQ(err_PETSc)
186-
do point = 1, nBcPoints
188+
do point = 1_pPETSCINT, nBcPoints
187189
call PetscSectionGetFieldDof(section,bcPoints(point),0_pPETSCINT,numDof,err_PETSc)
188190
CHKERRQ(err_PETSc)
189191
call PetscSectionGetFieldOffset(section,bcPoints(point),0_pPETSCINT,offset,err_PETSc)
190192
CHKERRQ(err_PETSc)
191193
do dof = offset+component, offset+numDof, numComp
192-
localArray(dof) = localArray(dof) + mechBC(face)%Value(component)*Delta_t
194+
localArray(dof) = localArray(dof) + mechBC(boundary)%dot_u(component)*Delta_t
193195
end do
194196
end do
195197
call VecRestoreArray(solution_local,localArray,err_PETSc)
@@ -198,7 +200,7 @@ subroutine utilities_projectBCValues(dm_local,solution_local,section,mechBC,Delt
198200
CHKERRQ(err_PETSc)
199201
call VecAssemblyEnd(solution_local, err_PETSc)
200202
CHKERRQ(err_PETSc)
201-
if (nBcPoints > 0) call ISRestoreIndices(bcPointsIS,bcPoints,err_PETSc)
203+
if (nBcPoints > 0_pPETSCINT) call ISRestoreIndices(bcPointsIS,bcPoints,err_PETSc)
202204
CHKERRQ(err_PETSc)
203205
call ISDestroy(bcPointsIS,err_PETSc)
204206
CHKERRQ(err_PETSc)

‎src/mesh/discretization_mesh.f90‎

Lines changed: 50 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,14 @@ module discretization_mesh
3838
private
3939

4040
PetscInt, public, protected :: &
41-
mesh_Nboundaries
41+
mesh_Nboundaries !< Number of defined BC (total)
42+
43+
PetscInt, dimension(2), public, protected :: &
44+
mesh_BCTypeSetSize !< Number of each BC type (1: Vertex, 2: Faces)
45+
46+
PetscInt, dimension(:), allocatable, public, protected :: &
47+
mesh_boundariesIS, & !< Index Set (tag values) of BC in mesh file
48+
mesh_boundariesIdx !< BC Type index (BCType_Vertex, BCType_Face)
4249

4350
PetscInt, public, protected :: &
4451
mesh_NcpElemsGlobal
@@ -51,10 +58,14 @@ module discretization_mesh
5158
mesh_maxNips !< max number of IPs in any CP element
5259
!!!! END DEPRECATED !!!!!
5360

54-
DM, public :: geomMesh
61+
PetscInt, parameter :: &
62+
BCTypeFace = 1_pPETSCINT, &
63+
BCTypeVertex = 2_pPETSCINT
5564

56-
PetscInt, dimension(:), allocatable, public, protected :: &
57-
mesh_boundaries
65+
character(len=11), dimension(2), public, protected :: &
66+
mesh_BCTypeLabel = ['Face Sets ', 'Vertex Sets'] !< PETSc default BC labels
67+
68+
DM, public :: geomMesh
5869

5970
real(pREAL), dimension(:,:), allocatable :: &
6071
mesh_ipVolume, & !< volume associated with IP (initially!)
@@ -63,6 +74,7 @@ module discretization_mesh
6374
real(pREAL), dimension(:,:,:), allocatable :: &
6475
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
6576

77+
6678
#if PETSC_VERSION_MINOR < 16
6779
external :: &
6880
DMDestroy
@@ -85,9 +97,10 @@ subroutine discretization_mesh_init()
8597
j
8698
PetscSF :: sf
8799
DM :: globalMesh
88-
PetscInt :: nFaceSets, Nboundaries, NelemsGlobal, Nelems
89-
PetscInt, pointer, dimension(:) :: pFaceSets
90-
IS :: faceSetIS
100+
PetscInt :: NelemsGlobal, Nelems
101+
PetscInt :: label
102+
IS :: setIS
103+
PetscInt, pointer, dimension(:) :: pSets
91104
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>=18)
92105
PetscQuadrature :: quadrature
93106
#endif
@@ -99,7 +112,7 @@ subroutine discretization_mesh_init()
99112
num_solver, &
100113
num_mesh
101114
PetscInt :: p_i
102-
integer:: dim !< integration order (quadrature rule)
115+
integer:: dim
103116
type(tvec) :: coords_node0
104117
real(pREAL), pointer, dimension(:) :: &
105118
mesh_node0_temp
@@ -118,6 +131,7 @@ subroutine discretization_mesh_init()
118131
num_mesh => num_solver%get_dict('mesh',defaultVal=emptyDict)
119132
p_i = num_mesh%get_asInt('p_i',defaultVal=2)
120133

134+
call PetscOptionsInsertString(PETSC_NULL_OPTIONS, ' -dm_plex_gmsh_mark_vertices ', err_PETSc)
121135
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>16)
122136
call DMPlexCreateFromFile(PETSC_COMM_WORLD,CLI_geomFile,'n/a',PETSC_TRUE,globalMesh,err_PETSc)
123137
#else
@@ -128,16 +142,8 @@ subroutine discretization_mesh_init()
128142
CHKERRQ(err_PETSc)
129143
call DMGetStratumSize(globalMesh,'depth',dimPlex,NelemsGlobal,err_PETSc)
130144
CHKERRQ(err_PETSc)
131-
mesh_NcpElemsGlobal = int(NelemsGlobal)
145+
mesh_NcpElemsGlobal = NelemsGlobal
132146

133-
! get number of IDs in face sets (for boundary conditions?)
134-
call DMGetLabelSize(globalMesh,'Face Sets',Nboundaries,err_PETSc)
135-
CHKERRQ(err_PETSc)
136-
mesh_Nboundaries = int(Nboundaries)
137-
call MPI_Bcast(mesh_Nboundaries,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
138-
call parallelization_chkerr(err_MPI)
139-
call MPI_Bcast(mesh_NcpElemsGlobal,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
140-
call parallelization_chkerr(err_MPI)
141147
dim = int(dimPlex)
142148
call MPI_Bcast(dim,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
143149
dimPlex = int(dim,pPETSCINT)
@@ -150,19 +156,35 @@ subroutine discretization_mesh_init()
150156
end if
151157
CHKERRQ(err_PETSc)
152158

153-
allocate(mesh_boundaries(mesh_Nboundaries), source = 0_pPETSCINT)
154-
call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,err_PETSc)
155-
CHKERRQ(err_PETSc)
156-
call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,err_PETSc)
157-
CHKERRQ(err_PETSc)
158-
if (nFaceSets > 0) then
159-
call ISGetIndices(faceSetIS,pFaceSets,err_PETSc)
159+
mesh_Nboundaries = 0_pPETSCINT
160+
allocate(mesh_boundariesIS(mesh_Nboundaries))
161+
mesh_BCTypeSetSize = 0_pPETSCINT
162+
do label = 1_pPETSCINT, size(mesh_BCTypeLabel)
163+
call DMGetLabelSize(globalMesh,mesh_BCTypeLabel(label),mesh_BCTypeSetSize(label),err_PETSc)
160164
CHKERRQ(err_PETSc)
161-
mesh_boundaries(1:nFaceSets) = pFaceSets
165+
call DMGetLabelIdIS(globalMesh,mesh_BCTypeLabel(label),setIS,err_PETSc)
162166
CHKERRQ(err_PETSc)
163-
call ISRestoreIndices(faceSetIS,pFaceSets,err_PETSc)
164-
end if
165-
call MPI_Bcast(mesh_boundaries,int(mesh_Nboundaries),MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
167+
if (mesh_BCTypeSetSize(label) > 0_pPETSCINT) then
168+
call ISGetIndices(setIS,pSets,err_PETSc)
169+
CHKERRQ(err_PETSc)
170+
mesh_boundariesIS = [mesh_boundariesIS, pSets]
171+
call ISRestoreIndices(setIS,pSets,err_PETSc)
172+
CHKERRQ(err_PETSc)
173+
mesh_Nboundaries = mesh_Nboundaries+mesh_BCTypeSetSize(label)
174+
end if
175+
end do
176+
allocate(mesh_boundariesIdx(mesh_Nboundaries), source = BCTypeFace)
177+
mesh_boundariesIdx(mesh_BCTypeSetSize(BCTypeFace)+1_pPETSCINT:mesh_Nboundaries) = BCTypeVertex
178+
179+
call MPI_Bcast(mesh_BCTypeSetSize,int(size(mesh_BCTypeSetSize),kind=MPI_INTEGER_KIND),MPI_INTEGER,0_MPI_INTEGER_KIND, &
180+
MPI_COMM_WORLD,err_MPI)
181+
call parallelization_chkerr(err_MPI)
182+
call MPI_Bcast(mesh_Nboundaries,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
183+
call parallelization_chkerr(err_MPI)
184+
call MPI_Bcast(mesh_boundariesIdx,int(mesh_Nboundaries,kind=MPI_INTEGER_KIND),MPI_INTEGER,0_MPI_INTEGER_KIND, &
185+
MPI_COMM_WORLD,err_MPI)
186+
call parallelization_chkerr(err_MPI)
187+
call MPI_Bcast(mesh_boundariesIS,int(mesh_Nboundaries),MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
166188
call parallelization_chkerr(err_MPI)
167189

168190
call DMDestroy(globalMesh,err_PETSc)

0 commit comments

Comments
 (0)