Added option to assign material properties using a majority or minority rule,...

Added option to assign material properties using a majority or minority rule, to complement what is done with DivFEM

src/pressure_cut.f90

@@ -16,10 +16,10 @@
 !------------------------------------------------------------------------------|
 !------------------------------------------------------------------------------|
 
-recursive subroutine pressure_cut (params,level,levelmax,icon,   &
-                                   x,y,z,mat,materialn,                   &
-                                   u,v,w,temp,pressure,strain,nnode,           &
-                                   lsf,nlsf,r0,s0,t0,rst,icut,ileaves,eviscosity)
+recursive subroutine pressure_cut (params,level,levelmax,icon,x,y,z,mat,       &
+                                   materialn,u,v,w,temp,pressure,strain,nnode, &
+                                   lsf,nlsf,r0,s0,t0,rst,icut,ileaves,         &
+                                   eviscosity)
                                    
 
 !------------------------------------------------------------------------------|
@@ -79,7 +79,7 @@ double precision eviscosity
 integer err,i,j,k,ii,jj,kk,levelp,matel,jcut
 double precision r(params%mpe),s(params%mpe),t(params%mpe)
 double precision h(params%mpe),vol_lsf0
-double precision r0p,s0p,t0p,rstp,dt
+double precision r0p,s0p,t0p,rstp,dt,eps
 double precision viscosity,density,penal,expon,diffusivity,heat
 double precision pressurep,prod,volmax
 double precision,dimension(:,:),allocatable::lsfp
@@ -90,6 +90,8 @@ double precision plasticity_parameters(9)
 !------------------------------------------------------------------------------|
 !------------------------------------------------------------------------------|
 
+eps=1.d-10
+
 if (level.eq.0) then
    pressure=0.d0
 endif
@@ -106,16 +108,13 @@ do i=1,nlsf
    endif
 enddo
 
-call make_pressure (params,icon,x,y,z, &
-                    mat(matel)%viscosity, &
-                    mat(matel)%penalty,mat(matel)%expon, &
-                    u,v,w,temp,pressurep,strain,nnode, &
-                    r0,s0,t0,rst,mat(matel)%plasticity_type,&
-                    mat(matel)%plasticity_parameters,eviscosity)
+call make_pressure (params,icon,x,y,z,mat(matel)%viscosity,mat(matel)%penalty, &
+                   mat(matel)%expon,u,v,w,temp,pressurep,strain,nnode,r0,s0,t0,&
+                   rst,mat(matel)%plasticity_type,&
+                   mat(matel)%plasticity_parameters,eviscosity)
 
 pressure=pressure+pressurep/(8.d0**level)
 
-
 icut=0
 
 return
@@ -128,7 +127,7 @@ icut=1
 ! we use an approximate algorithm
 
 if (level.eq.levelmax) then
-   allocate (vol_lsf(nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vol_lsf in make_cut$')
+   allocate (vol_lsf(nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vol_lsf in pressure_cut$')
    do i=1,nlsf
       call compute_positive_volume (lsf(1,i),vol_lsf(i),params%levelapprox)
    enddo
@@ -144,40 +143,68 @@ if (level.eq.levelmax) then
        vol_lsf(i)=max(vol_lsf(i),0.d0)
        vol_lsf(i)=min(vol_lsf(i),1.d0)
     enddo
-
-   viscosity=0.d0
-   penal=0.d0
-   expon=0.d0
-   do i=1,nlsf
-      matel=materialn(i)
-      viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
-      penal=penal+vol_lsf(i)*mat(matel)%penalty
-      expon=expon+vol_lsf(i)*mat(matel)%expon
-   enddo
-   vol_lsf0=1.d0-sum(vol_lsf)
-   matel=materialn(0)
-   viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
-   penal=penal+vol_lsf0*mat(matel)%penalty
-   expon=expon+vol_lsf0*mat(matel)%expon
-
-   volmax=0.d0
-   matel=params%materialn(0)
-   plasticity_type=mat(matel)%plasticity_type
-   plasticity_parameters=mat(matel)%plasticity_parameters
-   do i=1,nlsf
-      if (vol_lsf(i).gt.volmax) then
-         volmax=vol_lsf(i)
-         matel=params%materialn(i)
-         plasticity_type=mat(matel)%plasticity_type
-         plasticity_parameters=mat(matel)%plasticity_parameters
+    vol_lsf0=1.d0-sum(vol_lsf)
+
+    if (vol_lsf0.gt.eps) then                                                   ! Element contains void, always use divFEM
+      viscosity=0.d0
+      penal=0.d0
+      expon=0.d0
+      do i=1,nlsf
+        matel=materialn(i)
+        viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
+        penal=penal+vol_lsf(i)*mat(matel)%penalty
+        expon=expon+vol_lsf(i)*mat(matel)%expon
+      enddo
+      matel=materialn(0)
+      viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
+      penal=penal+vol_lsf0*mat(matel)%penalty
+      expon=expon+vol_lsf0*mat(matel)%expon
+    elseif (params%matrule.eq.1) then                                           ! Assign material properties of volumetric majority material
+      matel=params%materialn(sum(maxloc(vol_lsf)))
+      viscosity=mat(matel)%viscosity
+      penal=mat(matel)%penalty
+      expon=mat(matel)%expon
+    elseif (params%matrule.eq.2) then                                           ! Assign material properties of volumetric minority material
+      matel=params%materialn(sum(minloc(vol_lsf)))
+      viscosity=mat(matel)%viscosity
+      penal=mat(matel)%penalty
+      expon=mat(matel)%expon
+    else
+      if (params%matrule.ne.0) print *,'Invalid matrule value, using divFEM'
+      viscosity=0.d0
+      penal=0.d0
+      expon=0.d0
+      do i=1,nlsf
+        matel=materialn(i)
+        viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
+        penal=penal+vol_lsf(i)*mat(matel)%penalty
+        expon=expon+vol_lsf(i)*mat(matel)%expon
+      enddo
+      matel=materialn(0)
+      viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
+      penal=penal+vol_lsf0*mat(matel)%penalty
+      expon=expon+vol_lsf0*mat(matel)%expon
+    endif
+
+    if (maxval(vol_lsf).gt.0.d0) then
+      if (params%matrule.eq.2) then
+        matel=params%materialn(sum(minloc(vol_lsf)))
+        plasticity_type=mat(matel)%plasticity_type
+        plasticity_parameters=mat(matel)%plasticity_parameters
+      else
+        matel=params%materialn(sum(maxloc(vol_lsf)))
+        plasticity_type=mat(matel)%plasticity_type
+        plasticity_parameters=mat(matel)%plasticity_parameters
       endif
-   enddo
-
-
-   call make_pressure (params,icon,x,y,z, &
-                       viscosity,penal,expon, &
-                       u,v,w,temp,pressurep,strain,nnode, &
-                       r0,s0,t0,rst,plasticity_type,plasticity_parameters,eviscosity)
+    else
+      matel=params%materialn(0)
+      plasticity_type=mat(matel)%plasticity_type
+      plasticity_parameters=mat(matel)%plasticity_parameters
+    endif
+
+   call make_pressure (params,icon,x,y,z,viscosity,penal,expon,u,v,w,temp,     &
+                      pressurep,strain,nnode,r0,s0,t0,rst,plasticity_type,     &
+                      plasticity_parameters,eviscosity)
    pressure=pressure+pressurep/(8.d0**level)
    deallocate (vol_lsf)
    return
@@ -185,7 +212,7 @@ endif
 
 ! If we are not at the bottom level, we keep dividing
 
-allocate (lsfp(params%mpe,nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsfp in make_cut$')
+allocate (lsfp(params%mpe,nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsfp in pressure_cut$')
 
 do kk=1,2
    t(1:4)=-1.d0+float(kk-1)
@@ -225,11 +252,9 @@ do kk=1,2
          t0p=t0+rst*(t(1)+1.d0)/2.d0
          rstp=rst/2.d0
          levelp=level+1
-         call pressure_cut (params,levelp,levelmax, &
-                            icon,x,y,z, &
-                            mat,materialn, &
-                            u,v,w,temp,pressure,strain,nnode, &
-                            lsfp,nlsf,r0p,s0p,t0p,rstp,jcut,ileaves,eviscosity)
+         call pressure_cut (params,levelp,levelmax,icon,x,y,z,mat,materialn,u, &
+                           v,w,temp,pressure,strain,nnode,lsfp,nlsf,r0p,s0p,   &
+                           t0p,rstp,jcut,ileaves,eviscosity)
       enddo
    enddo
 enddo
@@ -239,4 +264,4 @@ deallocate (lsfp)
 return
 end
 !------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
\ No newline at end of file