DAMASK with grid solvers  Revision: v2.0.3-2204-gdb1f2151
The Düsseldorf Advanced Material Simulation Kit with Grid Solvers
homogenization_mech_RGC.f90
Go to the documentation of this file.
1 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
2 # 1 "<built-in>"
3 # 1 "<command-line>"
4 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
5 !--------------------------------------------------------------------------------------------------
12 !--------------------------------------------------------------------------------------------------
13 submodule(homogenization) homogenization_mech_rgc
14  use rotations
15 
16  type :: tparameters
17  integer, dimension(:), allocatable :: &
18  Nconstituents
19  real(pReal) :: &
20  xiAlpha, &
21  ciAlpha
22  real(pReal), dimension(:), allocatable :: &
23  dAlpha, &
24  angles
25  integer :: &
26  of_debug = 0
27  character(len=pStringLen), allocatable, dimension(:) :: &
28  output
29  end type tparameters
30 
31  type :: trgcstate
32  real(pReal), pointer, dimension(:) :: &
33  work, &
34  penaltyEnergy
35  real(pReal), pointer, dimension(:,:) :: &
36  relaxationVector
37  end type trgcstate
38 
39  type :: trgcdependentstate
40  real(pReal), allocatable, dimension(:) :: &
41  volumeDiscrepancy, &
42  relaxationRate_avg, &
43  relaxationRate_max
44  real(pReal), allocatable, dimension(:,:) :: &
45  mismatch
46  real(pReal), allocatable, dimension(:,:,:) :: &
47  orientation
48  end type trgcdependentstate
49 
50  type :: tnumerics_rgc
51  real(pReal) :: &
52  atol, & !< absolute tolerance of RGC residuum
53  rtol, & !< relative tolerance of RGC residuum
54  absMax, & !< absolute maximum of RGC residuum
55  relMax, & !< relative maximum of RGC residuum
56  pPert, & !< perturbation for computing RGC penalty tangent
57  xSmoo, & !< RGC penalty smoothing parameter (hyperbolic tangent)
58  viscPower, & !< power (sensitivity rate) of numerical viscosity in RGC scheme, Default 1.0e0: Newton viscosity (linear model)
59  viscModus, & !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied
60  refRelaxRate, & !< reference relaxation rate in RGC viscosity
61  maxdRelax, & !< threshold of maximum relaxation vector increment (if exceed this then cutback)
62  maxVolDiscr, & !< threshold of maximum volume discrepancy allowed
63  volDiscrMod, & !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
64  volDiscrPow
65  end type tnumerics_rgc
66 
67  type(tparameters), dimension(:), allocatable :: &
68  param
69  type(tRGCstate), dimension(:), allocatable :: &
70  state, &
71  state0
72  type(tRGCdependentState), dimension(:), allocatable :: &
73  dependentState
74  type(tNumerics_RGC) :: &
75  num ! numerics parameters. Better name?
76 
77 contains
78 
79 !--------------------------------------------------------------------------------------------------
81 !--------------------------------------------------------------------------------------------------
82 module subroutine mech_rgc_init
83 
84  integer :: &
85  Ninstance, &
86  h, &
87  NofMyHomog, &
88  sizeState, nIntFaceTot
89 
90  write(6,'(/,a)') ' <<<+- homogenization_'//homogenization_rgc_label//' init -+>>>'; flush(6)
91 
92  write(6,'(/,a)') ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
93  write(6,'(a)') ' https://doi.org/10.1007/s12289-009-0619-1'
94 
95  write(6,'(/,a)') ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
96  write(6,'(a)') ' https://doi.org/10.1088/0965-0393/18/1/015006'
97 
98  ninstance = count(homogenization_type == homogenization_rgc_id)
99  if (iand(debug_level(debug_homogenization),debug_levelbasic) /= 0) &
100  write(6,'(a16,1x,i5,/)') '# instances:',ninstance
101 
102  allocate(param(ninstance))
103  allocate(state(ninstance))
104  allocate(state0(ninstance))
105  allocate(dependentstate(ninstance))
106 
107  num%atol = config_numerics%getFloat('atol_rgc', defaultval=1.0e+4_preal)
108  num%rtol = config_numerics%getFloat('rtol_rgc', defaultval=1.0e-3_preal)
109  num%absMax = config_numerics%getFloat('amax_rgc', defaultval=1.0e+10_preal)
110  num%relMax = config_numerics%getFloat('rmax_rgc', defaultval=1.0e+2_preal)
111  num%pPert = config_numerics%getFloat('perturbpenalty_rgc', defaultval=1.0e-7_preal)
112  num%xSmoo = config_numerics%getFloat('relvantmismatch_rgc', defaultval=1.0e-5_preal)
113  num%viscPower = config_numerics%getFloat('viscositypower_rgc', defaultval=1.0e+0_preal)
114  num%viscModus = config_numerics%getFloat('viscositymodulus_rgc', defaultval=0.0e+0_preal)
115  num%refRelaxRate = config_numerics%getFloat('refrelaxationrate_rgc',defaultval=1.0e-3_preal)
116  num%maxdRelax = config_numerics%getFloat('maxrelaxationrate_rgc',defaultval=1.0e+0_preal)
117  num%maxVolDiscr = config_numerics%getFloat('maxvoldiscrepancy_rgc',defaultval=1.0e-5_preal)
118  num%volDiscrMod = config_numerics%getFloat('voldiscrepancymod_rgc',defaultval=1.0e+12_preal)
119  num%volDiscrPow = config_numerics%getFloat('dicrepancypower_rgc', defaultval=5.0_preal)
120 
121  if (num%atol <= 0.0_preal) call io_error(301,ext_msg='absTol_RGC')
122  if (num%rtol <= 0.0_preal) call io_error(301,ext_msg='relTol_RGC')
123  if (num%absMax <= 0.0_preal) call io_error(301,ext_msg='absMax_RGC')
124  if (num%relMax <= 0.0_preal) call io_error(301,ext_msg='relMax_RGC')
125  if (num%pPert <= 0.0_preal) call io_error(301,ext_msg='pPert_RGC')
126  if (num%xSmoo <= 0.0_preal) call io_error(301,ext_msg='xSmoo_RGC')
127  if (num%viscPower < 0.0_preal) call io_error(301,ext_msg='viscPower_RGC')
128  if (num%viscModus < 0.0_preal) call io_error(301,ext_msg='viscModus_RGC')
129  if (num%refRelaxRate <= 0.0_preal) call io_error(301,ext_msg='refRelaxRate_RGC')
130  if (num%maxdRelax <= 0.0_preal) call io_error(301,ext_msg='maxdRelax_RGC')
131  if (num%maxVolDiscr <= 0.0_preal) call io_error(301,ext_msg='maxVolDiscr_RGC')
132  if (num%volDiscrMod < 0.0_preal) call io_error(301,ext_msg='volDiscrMod_RGC')
133  if (num%volDiscrPow <= 0.0_preal) call io_error(301,ext_msg='volDiscrPw_RGC')
134 
135  do h = 1, size(homogenization_type)
136  if (homogenization_type(h) /= homogenization_rgc_id) cycle
137  associate(prm => param(homogenization_typeinstance(h)), &
138  stt => state(homogenization_typeinstance(h)), &
139  st0 => state0(homogenization_typeinstance(h)), &
140  dst => dependentstate(homogenization_typeinstance(h)), &
141  config => config_homogenization(h))
142 
143 
144 
145 
146 
147 
148 
149  prm%output = config%getStrings('(output)',defaultval=emptystringarray)
150 
151  prm%Nconstituents = config%getInts('clustersize',requiredsize=3)
152  if (homogenization_ngrains(h) /= product(prm%Nconstituents)) &
153  call io_error(211,ext_msg='clustersize ('//homogenization_rgc_label//')')
154 
155  prm%xiAlpha = config%getFloat('scalingparameter')
156  prm%ciAlpha = config%getFloat('overproportionality')
157 
158  prm%dAlpha = config%getFloats('grainsize', requiredsize=3)
159  prm%angles = config%getFloats('clusterorientation',requiredsize=3)
160 
161  nofmyhomog = count(material_homogenizationat == h)
162  nintfacetot = 3*( (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) &
163  + prm%Nconstituents(1)*(prm%Nconstituents(2)-1)*prm%Nconstituents(3) &
164  + prm%Nconstituents(1)*prm%Nconstituents(2)*(prm%Nconstituents(3)-1))
165  sizestate = nintfacetot &
166  + size(['avg constitutive work ','average penalty energy'])
167 
168  homogstate(h)%sizeState = sizestate
169  allocate(homogstate(h)%state0 (sizestate,nofmyhomog), source=0.0_preal)
170  allocate(homogstate(h)%subState0(sizestate,nofmyhomog), source=0.0_preal)
171  allocate(homogstate(h)%state (sizestate,nofmyhomog), source=0.0_preal)
172 
173  stt%relaxationVector => homogstate(h)%state(1:nintfacetot,:)
174  st0%relaxationVector => homogstate(h)%state0(1:nintfacetot,:)
175  stt%work => homogstate(h)%state(nintfacetot+1,:)
176  stt%penaltyEnergy => homogstate(h)%state(nintfacetot+2,:)
177 
178  allocate(dst%volumeDiscrepancy( nofmyhomog))
179  allocate(dst%relaxationRate_avg( nofmyhomog))
180  allocate(dst%relaxationRate_max( nofmyhomog))
181  allocate(dst%mismatch( 3,nofmyhomog))
182 
183 !--------------------------------------------------------------------------------------------------
184 ! assigning cluster orientations
185  dependentstate(homogenization_typeinstance(h))%orientation = spread(eu2om(prm%angles*inrad),3,nofmyhomog)
186  !dst%orientation = spread(eu2om(prm%angles*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason)
187 
188  end associate
189 
190  enddo
191 
192 end subroutine mech_rgc_init
193 
194 
195 !--------------------------------------------------------------------------------------------------
197 !--------------------------------------------------------------------------------------------------
198 module subroutine mech_rgc_partitiondeformation(f,avgf,instance,of)
199 
200  real(pReal), dimension (:,:,:), intent(out) :: F
201 
202  real(pReal), dimension (3,3), intent(in) :: avgF
203  integer, intent(in) :: &
204  instance, &
205  of
206 
207  real(pReal), dimension(3) :: aVect,nVect
208  integer, dimension(4) :: intFace
209  integer, dimension(3) :: iGrain3
210  integer :: iGrain,iFace,i,j
211 
212  associate(prm => param(instance))
213 
214 !--------------------------------------------------------------------------------------------------
215 ! compute the deformation gradient of individual grains due to relaxations
216  f = 0.0_preal
217  do igrain = 1,product(prm%Nconstituents)
218  igrain3 = grain1to3(igrain,prm%Nconstituents)
219  do iface = 1,6
220  intface = getinterface(iface,igrain3) ! identifying 6 interfaces of each grain
221  avect = relaxationvector(intface,instance,of) ! get the relaxation vectors for each interface from global relaxation vector array
222  nvect = interfacenormal(intface,instance,of)
223  forall (i=1:3,j=1:3) &
224  f(i,j,igrain) = f(i,j,igrain) + avect(i)*nvect(j) ! calculating deformation relaxations due to interface relaxation
225  enddo
226  f(1:3,1:3,igrain) = f(1:3,1:3,igrain) + avgf ! resulting relaxed deformation gradient
227 
228 # 234 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
229  enddo
230 
231  end associate
232 
233 end subroutine mech_rgc_partitiondeformation
234 
235 
236 !--------------------------------------------------------------------------------------------------
238 ! "happy" with result
239 !--------------------------------------------------------------------------------------------------
240 module procedure mech_rgc_updatestate
241 
242  integer, dimension(4) :: intFaceN,intFaceP,faceID
243  integer, dimension(3) :: nGDim,iGr3N,iGr3P
244  integer :: instance,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, of
245  real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD
246  real(pReal), dimension(3,size(P,3)) :: NN,devNull
247  real(pReal), dimension(3) :: normP,normN,mornP,mornN
248  real(pReal) :: residMax,stresMax
249  logical :: error
250  real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
251  real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
252 
253 
254 
255 
256 
257  zerotimestep: if(deq0(dt)) then
258  mech_rgc_updatestate = .true. ! pretend everything is fine and return
259  return
260  endif zerotimestep
261 
262  instance = homogenization_typeinstance(material_homogenizationat(el))
263  of = material_homogenizationmemberat(ip,el)
264 
265  associate(stt => state(instance), st0 => state0(instance), dst => dependentstate(instance), prm => param(instance))
266 
267 !--------------------------------------------------------------------------------------------------
268 ! get the dimension of the cluster (grains and interfaces)
269  ngdim = prm%Nconstituents
270  ngrain = product(ngdim)
271  nintfacetot = (ngdim(1)-1)*ngdim(2)*ngdim(3) &
272  + ngdim(1)*(ngdim(2)-1)*ngdim(3) &
273  + ngdim(1)*ngdim(2)*(ngdim(3)-1)
274 
275 !--------------------------------------------------------------------------------------------------
276 ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
277  allocate(resid(3*nintfacetot), source=0.0_preal)
278  allocate(tract(nintfacetot,3), source=0.0_preal)
279  relax = stt%relaxationVector(:,of)
280  drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of)
281 
282 # 296 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
283 
284 !--------------------------------------------------------------------------------------------------
285 ! computing interface mismatch and stress penalty tensor for all interfaces of all grains
286  call stresspenalty(r,nn,avgf,f,ip,el,instance,of)
287 
288 !--------------------------------------------------------------------------------------------------
289 ! calculating volume discrepancy and stress penalty related to overall volume discrepancy
290  call volumepenalty(d,dst%volumeDiscrepancy(of),avgf,f,ngrain,instance,of)
291 
292 # 320 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
293 
294 !------------------------------------------------------------------------------------------------
295 ! computing the residual stress from the balance of traction at all (interior) interfaces
296  do inum = 1,nintfacetot
297  faceid = interface1to4(inum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index)
298 
299 !--------------------------------------------------------------------------------------------------
300 ! identify the left/bottom/back grain (-|N)
301  igr3n = faceid(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index)
302  igrn = grain3to1(igr3n,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
303  intfacen = getinterface(2*faceid(1),igr3n)
304  normn = interfacenormal(intfacen,instance,of)
305 
306 !--------------------------------------------------------------------------------------------------
307 ! identify the right/up/front grain (+|P)
308  igr3p = igr3n
309  igr3p(faceid(1)) = igr3n(faceid(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index)
310  igrp = grain3to1(igr3p,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
311  intfacep = getinterface(2*faceid(1)-1,igr3p)
312  normp = interfacenormal(intfacep,instance,of)
313 
314 !--------------------------------------------------------------------------------------------------
315 ! compute the residual of traction at the interface (in local system, 4-dimensional index)
316  do i = 1,3
317  tract(inum,i) = sign(num%viscModus*(abs(drelax(i+3*(inum-1)))/(num%refRelaxRate*dt))**num%viscPower, &
318  drelax(i+3*(inum-1))) ! contribution from the relaxation viscosity
319  do j = 1,3
320  tract(inum,i) = tract(inum,i) + (p(i,j,igrp) + r(i,j,igrp) + d(i,j,igrp))*normp(j) & ! contribution from material stress P, mismatch penalty R, and volume penalty D projected into the interface
321  + (p(i,j,igrn) + r(i,j,igrn) + d(i,j,igrn))*normn(j)
322  resid(i+3*(inum-1)) = tract(inum,i) ! translate the local residual into global 1-dimensional residual array
323  enddo
324  enddo
325 
326 
327 
328 
329 
330 
331 
332 
333  enddo
334 
335 !--------------------------------------------------------------------------------------------------
336 ! convergence check for stress residual
337  stresmax = maxval(abs(p)) ! get the maximum of first Piola-Kirchhoff (material) stress
338  residmax = maxval(abs(tract)) ! get the maximum of the residual
339 
340 # 380 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
341 
342  mech_rgc_updatestate = .false.
343 
344 !--------------------------------------------------------------------------------------------------
345 ! If convergence reached => done and happy
346  if (residmax < num%rtol*stresmax .or. residmax < num%atol) then
347  mech_rgc_updatestate = .true.
348 
349 
350 
351 
352 
353 !--------------------------------------------------------------------------------------------------
354 ! compute/update the state for postResult, i.e., all energy densities computed by time-integration
355  do igrain = 1,product(prm%Nconstituents)
356  do i = 1,3;do j = 1,3
357  stt%work(of) = stt%work(of) &
358  + p(i,j,igrain)*(f(i,j,igrain) - f0(i,j,igrain))/real(ngrain,preal)
359  stt%penaltyEnergy(of) = stt%penaltyEnergy(of) &
360  + r(i,j,igrain)*(f(i,j,igrain) - f0(i,j,igrain))/real(ngrain,preal)
361  enddo; enddo
362  enddo
363 
364  dst%mismatch(1:3,of) = sum(nn,2)/real(ngrain,preal)
365  dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nintfacetot,preal)
366  dst%relaxationRate_max(of) = maxval(abs(drelax))/dt
367 
368 # 420 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
369 
370  return
371 
372 !--------------------------------------------------------------------------------------------------
373 ! if residual blows-up => done but unhappy
374  elseif (residmax > num%relMax*stresmax .or. residmax > num%absMax) then ! try to restart when residual blows up exceeding maximum bound
375  mech_rgc_updatestate = [.true.,.false.] ! with direct cut-back
376 
377 
378 
379 
380 
381 
382  return
383 
384  else ! proceed with computing the Jacobian and state update
385 
386 
387 
388 
389 
390  endif
391 
392 !---------------------------------------------------------------------------------------------------
393 ! construct the global Jacobian matrix for updating the global relaxation vector array when
394 ! convergence is not yet reached ...
395 
396 !--------------------------------------------------------------------------------------------------
397 ! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix"
398  allocate(smatrix(3*nintfacetot,3*nintfacetot), source=0.0_preal)
399  do inum = 1,nintfacetot
400  faceid = interface1to4(inum,param(instance)%Nconstituents) ! assembling of local dPdF into global Jacobian matrix
401 
402 !--------------------------------------------------------------------------------------------------
403 ! identify the left/bottom/back grain (-|N)
404  igr3n = faceid(2:4) ! identifying the grain ID in local coordinate sytem
405  igrn = grain3to1(igr3n,param(instance)%Nconstituents) ! translate into global grain ID
406  intfacen = getinterface(2*faceid(1),igr3n) ! identifying the connecting interface in local coordinate system
407  normn = interfacenormal(intfacen,instance,of)
408  do iface = 1,6
409  intfacen = getinterface(iface,igr3n) ! identifying all interfaces that influence relaxation of the above interface
410  mornn = interfacenormal(intfacen,instance,of)
411  imun = interface4to1(intfacen,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index
412  if (imun > 0) then ! get the corresponding tangent
413  do i=1,3; do j=1,3; do k=1,3; do l=1,3
414  smatrix(3*(inum-1)+i,3*(imun-1)+j) = smatrix(3*(inum-1)+i,3*(imun-1)+j) &
415  + dpdf(i,k,j,l,igrn)*normn(k)*mornn(l)
416  enddo;enddo;enddo;enddo
417 ! projecting the material tangent dPdF into the interface
418 ! to obtain the Jacobian matrix contribution of dPdF
419  endif
420  enddo
421 
422 !--------------------------------------------------------------------------------------------------
423 ! identify the right/up/front grain (+|P)
424  igr3p = igr3n
425  igr3p(faceid(1)) = igr3n(faceid(1))+1 ! identifying the grain ID in local coordinate sytem
426  igrp = grain3to1(igr3p,param(instance)%Nconstituents) ! translate into global grain ID
427  intfacep = getinterface(2*faceid(1)-1,igr3p) ! identifying the connecting interface in local coordinate system
428  normp = interfacenormal(intfacep,instance,of)
429  do iface = 1,6
430  intfacep = getinterface(iface,igr3p) ! identifying all interfaces that influence relaxation of the above interface
431  mornp = interfacenormal(intfacep,instance,of)
432  imun = interface4to1(intfacep,param(instance)%Nconstituents) ! translate the interfaces ID into local 4-dimensional index
433  if (imun > 0) then ! get the corresponding tangent
434  do i=1,3; do j=1,3; do k=1,3; do l=1,3
435  smatrix(3*(inum-1)+i,3*(imun-1)+j) = smatrix(3*(inum-1)+i,3*(imun-1)+j) &
436  + dpdf(i,k,j,l,igrp)*normp(k)*mornp(l)
437  enddo;enddo;enddo;enddo
438  endif
439  enddo
440  enddo
441 
442 # 503 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
443 
444 !--------------------------------------------------------------------------------------------------
445 ! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical
446 ! perturbation method) "pmatrix"
447  allocate(pmatrix(3*nintfacetot,3*nintfacetot), source=0.0_preal)
448  allocate(p_relax(3*nintfacetot), source=0.0_preal)
449  allocate(p_resid(3*nintfacetot), source=0.0_preal)
450 
451  do ipert = 1,3*nintfacetot
452  p_relax = relax
453  p_relax(ipert) = relax(ipert) + num%pPert ! perturb the relaxation vector
454  stt%relaxationVector(:,of) = p_relax
455  call graindeformation(pf,avgf,instance,of) ! rain deformation from perturbed state
456  call stresspenalty(pr,devnull, avgf,pf,ip,el,instance,of) ! stress penalty due to interface mismatch from perturbed state
457  call volumepenalty(pd,devnull(1,1), avgf,pf,ngrain,instance,of) ! stress penalty due to volume discrepancy from perturbed state
458 
459 !--------------------------------------------------------------------------------------------------
460 ! computing the global stress residual array from the perturbed state
461  p_resid = 0.0_preal
462  do inum = 1,nintfacetot
463  faceid = interface1to4(inum,param(instance)%Nconstituents) ! identifying the interface ID in local coordinate system (4-dimensional index)
464 
465 !--------------------------------------------------------------------------------------------------
466 ! identify the left/bottom/back grain (-|N)
467  igr3n = faceid(2:4) ! identify the grain ID in local coordinate system (3-dimensional index)
468  igrn = grain3to1(igr3n,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
469  intfacen = getinterface(2*faceid(1),igr3n) ! identify the interface ID of the grain
470  normn = interfacenormal(intfacen,instance,of)
471 
472 !--------------------------------------------------------------------------------------------------
473 ! identify the right/up/front grain (+|P)
474  igr3p = igr3n
475  igr3p(faceid(1)) = igr3n(faceid(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index)
476  igrp = grain3to1(igr3p,param(instance)%Nconstituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
477  intfacep = getinterface(2*faceid(1)-1,igr3p) ! identify the interface ID of the grain
478  normp = interfacenormal(intfacep,instance,of)
479 
480 !--------------------------------------------------------------------------------------------------
481 ! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state
482 ! at all interfaces
483  do i = 1,3; do j = 1,3
484  p_resid(i+3*(inum-1)) = p_resid(i+3*(inum-1)) + (pr(i,j,igrp) - r(i,j,igrp))*normp(j) &
485  + (pr(i,j,igrn) - r(i,j,igrn))*normn(j) &
486  + (pd(i,j,igrp) - d(i,j,igrp))*normp(j) &
487  + (pd(i,j,igrn) - d(i,j,igrn))*normn(j)
488  enddo; enddo
489  enddo
490  pmatrix(:,ipert) = p_resid/num%pPert
491  enddo
492 
493 # 563 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
494 
495 !--------------------------------------------------------------------------------------------------
496 ! ... of the numerical viscosity traction "rmatrix"
497  allocate(rmatrix(3*nintfacetot,3*nintfacetot),source=0.0_preal)
498  do i=1,3*nintfacetot
499  rmatrix(i,i) = num%viscModus*num%viscPower/(num%refRelaxRate*dt)* & ! tangent due to numerical viscosity traction appears
500  (abs(drelax(i))/(num%refRelaxRate*dt))**(num%viscPower - 1.0_preal) ! only in the main diagonal term
501  enddo
502 
503 # 582 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
504 
505 !--------------------------------------------------------------------------------------------------
506 ! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix
507  allocate(jmatrix(3*nintfacetot,3*nintfacetot)); jmatrix = smatrix + pmatrix + rmatrix
508 
509 # 597 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
510 
511 !--------------------------------------------------------------------------------------------------
512 ! computing the update of the state variable (relaxation vectors) using the Jacobian matrix
513  allocate(jnverse(3*nintfacetot,3*nintfacetot),source=0.0_preal)
514  call math_invert(jnverse,error,jmatrix)
515 
516 # 613 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
517 
518 !--------------------------------------------------------------------------------------------------
519 ! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
520  drelax = 0.0_preal
521  do i = 1,3*nintfacetot;do j = 1,3*nintfacetot
522  drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable
523  enddo; enddo
524  stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration
525  if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
526  mech_rgc_updatestate = [.true.,.false.]
527  !$OMP CRITICAL (write2out)
528  write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
529  write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax))
530  flush(6)
531  !$OMP END CRITICAL (write2out)
532  endif
533 
534 # 640 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
535 
536  end associate
537 
538  contains
539  !------------------------------------------------------------------------------------------------
541  !------------------------------------------------------------------------------------------------
542  subroutine stresspenalty(rPen,nMis,avgF,fDef,ip,el,instance,of)
543 
544  real(pReal), dimension (:,:,:), intent(out) :: rPen
545  real(pReal), dimension (:,:), intent(out) :: nMis
546 
547  real(pReal), dimension (:,:,:), intent(in) :: fDef
548  real(pReal), dimension (3,3), intent(in) :: avgF
549  integer, intent(in) :: ip,el,instance,of
550 
551  integer, dimension (4) :: intFace
552  integer, dimension (3) :: iGrain3,iGNghb3,nGDim
553  real(pReal), dimension (3,3) :: gDef,nDef
554  real(pReal), dimension (3) :: nVect,surfCorr
555  real(pReal), dimension (2) :: Gmoduli
556  integer :: iGrain,iGNghb,iFace,i,j,k,l
557  real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
558  real(pReal), parameter :: nDefToler = 1.0e-10_preal
559 
560 
561 
562 
563  ngdim = param(instance)%Nconstituents
564  rpen = 0.0_preal
565  nmis = 0.0_preal
566 
567  !----------------------------------------------------------------------------------------------
568  ! get the correction factor the modulus of penalty stress representing the evolution of area of
569  ! the interfaces due to deformations
570 
571  surfcorr = surfacecorrection(avgf,instance,of)
572 
573  associate(prm => param(instance))
574 
575 # 688 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/homogenization_mech_RGC.f90"
576 
577  !-----------------------------------------------------------------------------------------------
578  ! computing the mismatch and penalty stress tensor of all grains
579  grainloop: do igrain = 1,product(prm%Nconstituents)
580  gmoduli = equivalentmoduli(igrain,ip,el)
581  mugrain = gmoduli(1) ! collecting the equivalent shear modulus of grain
582  bggrain = gmoduli(2) ! and the lengthh of Burgers vector
583  igrain3 = grain1to3(igrain,prm%Nconstituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position
584 
585  interfaceloop: do iface = 1,6
586  intface = getinterface(iface,igrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain
587  nvect = interfacenormal(intface,instance,of)
588  ignghb3 = igrain3 ! identify the neighboring grain across the interface
589  ignghb3(abs(intface(1))) = ignghb3(abs(intface(1))) &
590  + int(real(intface(1),preal)/real(abs(intface(1)),preal))
591  where(ignghb3 < 1) ignghb3 = ngdim
592  where(ignghb3 >ngdim) ignghb3 = 1
593  ignghb = grain3to1(ignghb3,prm%Nconstituents) ! get the ID of the neighboring grain
594  gmoduli = equivalentmoduli(ignghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor
595  mugnghb = gmoduli(1)
596  bggnghb = gmoduli(2)
597  gdef = 0.5_preal*(fdef(1:3,1:3,ignghb) - fdef(1:3,1:3,igrain)) ! difference/jump in deformation gradeint across the neighbor
598 
599  !-------------------------------------------------------------------------------------------
600  ! compute the mismatch tensor of all interfaces
601  ndefnorm = 0.0_preal
602  ndef = 0.0_preal
603  do i = 1,3; do j = 1,3
604  do k = 1,3; do l = 1,3
605  ndef(i,j) = ndef(i,j) - nvect(k)*gdef(i,l)*math_levicivita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient
606  enddo; enddo
607  ndefnorm = ndefnorm + ndef(i,j)**2.0_preal ! compute the norm of the mismatch tensor
608  enddo; enddo
609  ndefnorm = max(ndeftoler,sqrt(ndefnorm)) ! approximation to zero mismatch if mismatch is zero (singularity)
610  nmis(abs(intface(1)),igrain) = nmis(abs(intface(1)),igrain) + ndefnorm ! total amount of mismatch experienced by the grain (at all six interfaces)
611 
612 
613 
614 
615 
616 
617 
618 
619  !-------------------------------------------------------------------------------------------
620  ! compute the stress penalty of all interfaces
621  do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3
622  rpen(i,j,igrain) = rpen(i,j,igrain) + 0.5_preal*(mugrain*bggrain + mugnghb*bggnghb)*prm%xiAlpha &
623  *surfcorr(abs(intface(1)))/prm%dAlpha(abs(intface(1))) &
624  *cosh(prm%ciAlpha*ndefnorm) &
625  *0.5_preal*nvect(l)*ndef(i,k)/ndefnorm*math_levicivita(k,l,j) &
626  *tanh(ndefnorm/num%xSmoo)
627  enddo; enddo;enddo; enddo
628  enddo interfaceloop
629 
630 
631 
632 
633 
634 
635 
636  enddo grainloop
637 
638  end associate
639 
640  end subroutine stresspenalty
641 
642 
643  !------------------------------------------------------------------------------------------------
645  !------------------------------------------------------------------------------------------------
646  subroutine volumepenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of)
647 
648  real(pReal), dimension (:,:,:), intent(out) :: vPen ! stress-like penalty due to volume
649  real(pReal), intent(out) :: vDiscrep ! total volume discrepancy
650 
651  real(pReal), dimension (:,:,:), intent(in) :: fDef ! deformation gradients
652  real(pReal), dimension (3,3), intent(in) :: fAvg ! overall deformation gradient
653  integer, intent(in) :: &
654  Ngrain, &
655  instance, &
656  of
657 
658  real(pReal), dimension(size(vPen,3)) :: gVol
659  integer :: i
660 
661  !----------------------------------------------------------------------------------------------
662  ! compute the volumes of grains and of cluster
663  vdiscrep = math_det33(favg) ! compute the volume of the cluster
664  do i = 1,ngrain
665  gvol(i) = math_det33(fdef(1:3,1:3,i)) ! compute the volume of individual grains
666  vdiscrep = vdiscrep - gvol(i)/real(ngrain,preal) ! calculate the difference/dicrepancy between
667  ! the volume of the cluster and the the total volume of grains
668  enddo
669 
670  !----------------------------------------------------------------------------------------------
671  ! calculate the stress and penalty due to volume discrepancy
672  vpen = 0.0_preal
673  do i = 1,ngrain
674  vpen(:,:,i) = -1.0_preal/real(ngrain,preal)*num%volDiscrMod*num%volDiscrPow/num%maxVolDiscr* &
675  sign((abs(vdiscrep)/num%maxVolDiscr)**(num%volDiscrPow - 1.0),vdiscrep)* &
676  gvol(i)*transpose(math_inv33(fdef(:,:,i)))
677 
678 
679 
680 
681 
682 
683 
684 
685  enddo
686 
687  end subroutine volumepenalty
688 
689 
690  !--------------------------------------------------------------------------------------------------
692  ! deformation
693  !--------------------------------------------------------------------------------------------------
694  function surfacecorrection(avgF,instance,of)
695 
696  real(pReal), dimension(3) :: surfaceCorrection
697 
698  real(pReal), dimension(3,3), intent(in) :: avgF
699  integer, intent(in) :: &
700  instance, &
701  of
702  real(pReal), dimension(3,3) :: invC
703  real(pReal), dimension(3) :: nVect
704  real(pReal) :: detF
705  integer :: i,j,iBase
706  logical :: error
707 
708  call math_invert33(invc,detf,error,matmul(transpose(avgf),avgf))
709 
710  surfacecorrection = 0.0_preal
711  do ibase = 1,3
712  nvect = interfacenormal([ibase,1,1,1],instance,of)
713  do i = 1,3; do j = 1,3
714  surfacecorrection(ibase) = surfacecorrection(ibase) + invc(i,j)*nvect(i)*nvect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal
715  enddo; enddo
716  surfacecorrection(ibase) = sqrt(surfacecorrection(ibase))*detf ! get the surface correction factor (area contraction/enlargement)
717  enddo
718 
719  end function surfacecorrection
720 
721 
722  !--------------------------------------------------------------------------------------------------
724  !--------------------------------------------------------------------------------------------------
725  function equivalentmoduli(grainID,ip,el)
726 
727  real(pReal), dimension(2) :: equivalentModuli
728 
729  integer, intent(in) :: &
730  grainID,&
731  ip, & !< integration point number
732  el
733  real(pReal), dimension(6,6) :: elasTens
734  real(pReal) :: &
735  cEquiv_11, &
736  cEquiv_12, &
737  cEquiv_44
738 
739  elastens = constitutive_homogenizedc(grainid,ip,el)
740 
741  !----------------------------------------------------------------------------------------------
742  ! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005)
743  cequiv_11 = (elastens(1,1) + elastens(2,2) + elastens(3,3))/3.0_preal
744  cequiv_12 = (elastens(1,2) + elastens(2,3) + elastens(3,1) + &
745  elastens(1,3) + elastens(2,1) + elastens(3,2))/6.0_preal
746  cequiv_44 = (elastens(4,4) + elastens(5,5) + elastens(6,6))/3.0_preal
747  equivalentmoduli(1) = 0.2_preal*(cequiv_11 - cequiv_12) + 0.6_preal*cequiv_44
748 
749  !----------------------------------------------------------------------------------------------
750  ! obtain the length of Burgers vector (could be model dependend)
751  equivalentmoduli(2) = 2.5e-10_preal
752 
753  end function equivalentmoduli
754 
755 
756  !--------------------------------------------------------------------------------------------------
758  ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
759  !--------------------------------------------------------------------------------------------------
760  subroutine graindeformation(F, avgF, instance, of)
761 
762  real(pReal), dimension(:,:,:), intent(out) :: F
763 
764  real(pReal), dimension(:,:), intent(in) :: avgF
765  integer, intent(in) :: &
766  instance, &
767  of
768 
769  real(pReal), dimension(3) :: aVect,nVect
770  integer, dimension(4) :: intFace
771  integer, dimension(3) :: iGrain3
772  integer :: iGrain,iFace,i,j
773 
774  !-------------------------------------------------------------------------------------------------
775  ! compute the deformation gradient of individual grains due to relaxations
776 
777  associate(prm => param(instance))
778 
779  f = 0.0_preal
780  do igrain = 1,product(prm%Nconstituents)
781  igrain3 = grain1to3(igrain,prm%Nconstituents)
782  do iface = 1,6
783  intface = getinterface(iface,igrain3)
784  avect = relaxationvector(intface,instance,of)
785  nvect = interfacenormal(intface,instance,of)
786  forall (i=1:3,j=1:3) &
787  f(i,j,igrain) = f(i,j,igrain) + avect(i)*nvect(j) ! effective relaxations
788  enddo
789  f(1:3,1:3,igrain) = f(1:3,1:3,igrain) + avgf ! relaxed deformation gradient
790  enddo
791 
792  end associate
793 
794  end subroutine graindeformation
795 
796 end procedure mech_rgc_updatestate
797 
798 
799 !--------------------------------------------------------------------------------------------------
801 !--------------------------------------------------------------------------------------------------
802 module subroutine mech_rgc_averagestressanditstangent(avgp,davgpdavgf,p,dpdf,instance)
803 
804  real(pReal), dimension (3,3), intent(out) :: avgP
805  real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF
806 
807  real(pReal), dimension (:,:,:), intent(in) :: P
808  real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF
809  integer, intent(in) :: instance
810 
811  avgp = sum(p,3) /real(product(param(instance)%Nconstituents),preal)
812  davgpdavgf = sum(dpdf,5)/real(product(param(instance)%Nconstituents),preal)
813 
814 end subroutine mech_rgc_averagestressanditstangent
815 
816 
817 !--------------------------------------------------------------------------------------------------
819 !--------------------------------------------------------------------------------------------------
820 module subroutine mech_rgc_results(instance,group)
821 
822  integer, intent(in) :: instance
823  character(len=*), intent(in) :: group
824 
825  integer :: o
826 
827  associate(stt => state(instance), dst => dependentstate(instance), prm => param(instance))
828  outputsloop: do o = 1,size(prm%output)
829  select case(trim(prm%output(o)))
830  case('constitutivework')
831  call results_writedataset(group,stt%work,'W',&
832  'work density','J/m³')
833  case('magnitudemismatch')
834  call results_writedataset(group,dst%mismatch,'N',&
835  'average mismatch tensor','1')
836  case('penaltyenergy')
837  call results_writedataset(group,stt%penaltyEnergy,'R',&
838  'mismatch penalty density','J/m³')
839  case('volumediscrepancy')
840  call results_writedataset(group,dst%volumeDiscrepancy,'Delta_V',&
841  'volume discrepancy','m³')
842  case('maximumrelaxrate')
843  call results_writedataset(group,dst%relaxationrate_max,'max_alpha_dot',&
844  'maximum relaxation rate','m/s')
845  case('averagerelaxrate')
846  call results_writedataset(group,dst%relaxationrate_avg,'avg_alpha_dot',&
847  'average relaxation rate','m/s')
848  end select
849  enddo outputsloop
850  end associate
851 
852 end subroutine mech_rgc_results
853 
854 
855 !--------------------------------------------------------------------------------------------------
857 !--------------------------------------------------------------------------------------------------
858 pure function relaxationvector(intFace,instance,of)
859 
860  real(pReal), dimension (3) :: relaxationVector
861 
862  integer, intent(in) :: instance,of
863  integer, dimension(4), intent(in) :: intFace
864 
865  integer :: iNum
866 
867 !--------------------------------------------------------------------------------------------------
868 ! collect the interface relaxation vector from the global state array
869 
870  inum = interface4to1(intface,param(instance)%Nconstituents) ! identify the position of the interface in global state array
871  if (inum > 0) then
872  relaxationvector = state(instance)%relaxationVector((3*inum-2):(3*inum),of)
873  else
874  relaxationvector = 0.0_preal
875  endif
876 
877 end function relaxationvector
878 
879 
880 !--------------------------------------------------------------------------------------------------
882 !--------------------------------------------------------------------------------------------------
883 pure function interfacenormal(intFace,instance,of)
884 
885  real(pReal), dimension(3) :: interfaceNormal
886 
887  integer, dimension(4), intent(in) :: intFace
888  integer, intent(in) :: &
889  instance, &
890  of
891 
892  integer :: nPos
893 
894 !--------------------------------------------------------------------------------------------------
895 ! get the normal of the interface, identified from the value of intFace(1)
896  interfacenormal = 0.0_preal
897  npos = abs(intface(1)) ! identify the position of the interface in global state array
898  interfacenormal(npos) = real(intface(1)/abs(intface(1)),preal) ! get the normal vector w.r.t. cluster axis
899 
900  interfacenormal = matmul(dependentstate(instance)%orientation(1:3,1:3,of),interfacenormal) ! map the normal vector into sample coordinate system (basis)
901 
902 end function interfacenormal
903 
904 
905 !--------------------------------------------------------------------------------------------------
907 !--------------------------------------------------------------------------------------------------
908 pure function getinterface(iFace,iGrain3)
909 
910  integer, dimension(4) :: getInterface
911 
912  integer, dimension(3), intent(in) :: iGrain3
913  integer, intent(in) :: iFace
914 
915  integer :: iDir
916 
917  idir = (int(real(iface-1,preal)/2.0_preal)+1)*(-1)**iface
918  getinterface(1) = idir
919 
920 !--------------------------------------------------------------------------------------------------
921 ! identify the interface position by the direction of its normal
922  getinterface(2:4) = igrain3
923  if (idir < 0) getinterface(1-idir) = getinterface(1-idir)-1 ! to have a correlation with coordinate/position in real space
924 
925 end function getinterface
926 
927 
928 !--------------------------------------------------------------------------------------------------
930 !--------------------------------------------------------------------------------------------------
931 pure function grain1to3(grain1,nGDim)
932 
933  integer, dimension(3) :: grain1to3
934 
935  integer, intent(in) :: grain1
936  integer, dimension(3), intent(in) :: nGDim
937 
938  grain1to3 = 1 + [mod((grain1-1), ngdim(1)), &
939  mod((grain1-1)/ ngdim(1),ngdim(2)), &
940  (grain1-1)/(ngdim(1)*ngdim(2))]
941 
942 end function grain1to3
943 
944 
945 !--------------------------------------------------------------------------------------------------
947 !--------------------------------------------------------------------------------------------------
948 integer pure function grain3to1(grain3,nGDim)
949 
950  integer, dimension(3), intent(in) :: grain3
951  integer, dimension(3), intent(in) :: nGDim
952 
953  grain3to1 = grain3(1) &
954  + ngdim(1)*(grain3(2)-1) &
955  + ngdim(1)*ngdim(2)*(grain3(3)-1)
956 
957 end function grain3to1
958 
959 
960 !--------------------------------------------------------------------------------------------------
962 !--------------------------------------------------------------------------------------------------
963 integer pure function interface4to1(iFace4D, nGDim)
964 
965  integer, dimension(4), intent(in) :: iFace4D
966  integer, dimension(3), intent(in) :: nGDim
967 
968 
969  select case(abs(iface4d(1)))
970 
971  case(1)
972  if ((iface4d(2) == 0) .or. (iface4d(2) == ngdim(1))) then
973  interface4to1 = 0
974  else
975  interface4to1 = iface4d(3) + ngdim(2)*(iface4d(4)-1) &
976  + ngdim(2)*ngdim(3)*(iface4d(2)-1)
977  endif
978 
979  case(2)
980  if ((iface4d(3) == 0) .or. (iface4d(3) == ngdim(2))) then
981  interface4to1 = 0
982  else
983  interface4to1 = iface4d(4) + ngdim(3)*(iface4d(2)-1) &
984  + ngdim(3)*ngdim(1)*(iface4d(3)-1) &
985  + (ngdim(1)-1)*ngdim(2)*ngdim(3) ! total # of interfaces normal || e1
986  endif
987 
988  case(3)
989  if ((iface4d(4) == 0) .or. (iface4d(4) == ngdim(3))) then
990  interface4to1 = 0
991  else
992  interface4to1 = iface4d(2) + ngdim(1)*(iface4d(3)-1) &
993  + ngdim(1)*ngdim(2)*(iface4d(4)-1) &
994  + (ngdim(1)-1)*ngdim(2)*ngdim(3) & ! total # of interfaces normal || e1
995  + ngdim(1)*(ngdim(2)-1)*ngdim(3) ! total # of interfaces normal || e2
996  endif
997 
998  case default
999  interface4to1 = -1
1000 
1001  end select
1002 
1003 end function interface4to1
1004 
1005 
1006 !--------------------------------------------------------------------------------------------------
1008 !--------------------------------------------------------------------------------------------------
1009 pure function interface1to4(iFace1D, nGDim)
1010 
1011  integer, dimension(4) :: interface1to4
1012 
1013  integer, intent(in) :: iFace1D
1014  integer, dimension(3), intent(in) :: nGDim
1015  integer, dimension(3) :: nIntFace
1016 
1017 !--------------------------------------------------------------------------------------------------
1018 ! compute the total number of interfaces, which ...
1019  nintface = [(ngdim(1)-1)*ngdim(2)*ngdim(3), & ! ... normal || e1
1020  ngdim(1)*(ngdim(2)-1)*ngdim(3), & ! ... normal || e2
1021  ngdim(1)*ngdim(2)*(ngdim(3)-1)] ! ... normal || e3
1022 
1023 !--------------------------------------------------------------------------------------------------
1024 ! get the corresponding interface ID in 4D (normal and local position)
1025  if (iface1d > 0 .and. iface1d <= nintface(1)) then ! interface with normal || e1
1026  interface1to4(1) = 1
1027  interface1to4(3) = mod((iface1d-1),ngdim(2))+1
1028  interface1to4(4) = mod(int(real(iface1d-1,preal)/real(ngdim(2),preal)),ngdim(3))+1
1029  interface1to4(2) = int(real(iface1d-1,preal)/real(ngdim(2),preal)/real(ngdim(3),preal))+1
1030  elseif (iface1d > nintface(1) .and. iface1d <= (nintface(2) + nintface(1))) then ! interface with normal || e2
1031  interface1to4(1) = 2
1032  interface1to4(4) = mod((iface1d-nintface(1)-1),ngdim(3))+1
1033  interface1to4(2) = mod(int(real(iface1d-nintface(1)-1,preal)/real(ngdim(3),preal)),ngdim(1))+1
1034  interface1to4(3) = int(real(iface1d-nintface(1)-1,preal)/real(ngdim(3),preal)/real(ngdim(1),preal))+1
1035  elseif (iface1d > nintface(2) + nintface(1) .and. iface1d <= (nintface(3) + nintface(2) + nintface(1))) then ! interface with normal || e3
1036  interface1to4(1) = 3
1037  interface1to4(2) = mod((iface1d-nintface(2)-nintface(1)-1),ngdim(1))+1
1038  interface1to4(3) = mod(int(real(iface1d-nintface(2)-nintface(1)-1,preal)/real(ngdim(1),preal)),ngdim(2))+1
1039  interface1to4(4) = int(real(iface1d-nintface(2)-nintface(1)-1,preal)/real(ngdim(1),preal)/real(ngdim(2),preal))+1
1040  endif
1041 
1042 end function interface1to4
1043 
1044 
1045 end submodule homogenization_mech_rgc
rotations
rotation storage and conversion
Definition: rotations.f90:53
config
Reads in the material configuration from file.
Definition: config.f90:13
homogenization
homogenization manager, organizing deformation partitioning and stress homogenization
Definition: homogenization.f90:11
rotations::eu2om
pure real(preal) function, dimension(3, 3), public eu2om(eu)
Euler angles to orientation matrix.
Definition: rotations.f90:702