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/make_cut.f90

@@ -32,7 +32,7 @@ recursive subroutine make_cut (level,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,&
 ! cell is further cut and make_cut is recursively called
 ! if we are in a cut cell and level is equal to levelmax, we call make_matrix
 ! with material properties that have been interpolated from the
-! various materail properties contribnuting to the cut cell
+! various material properties contribnuting to the cut cell
 
 ! level       : current level in cut cell algorithm. It varies between 0 and levelmax.
 ! levelapprox : used to improve the postitive volutme calculation by further division
@@ -58,7 +58,6 @@ recursive subroutine make_cut (level,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,&
 ! icut        : returned, 0 if homogeneous element - 1 if cut element
 ! ileaves     : current leaf number (useless except for debugging)
 
-
 !------------------------------------------------------------------------------|
 !((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
 !------------------------------------------------------------------------------|
@@ -103,7 +102,7 @@ double precision volmax
 double precision r(params%mpe),s(params%mpe),t(params%mpe)
 double precision aelp(params%mpe*ndof,params%mpe*ndof),belp(params%mpe*ndof)
 double precision h(params%mpe),vol_lsf0,prod
-double precision r0p,s0p,t0p,rstp
+double precision r0p,s0p,t0p,rstp,eps
 double precision viscosity,density,penal,expon,diffusivity,heat,activ,expan
 character (len=8) plasticity_type
 double precision plasticity_parameters(9)
@@ -114,9 +113,12 @@ integer matel
 double precision eviscosityp
 logical  is_plastic_temp,vbounded_temp
 double precision weight
+
 !-------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------
 
+eps=1.d-10
+
 if (level.eq.0) then
    ael=0.d0
    bel=0.d0
@@ -156,16 +158,13 @@ end if
 
 !=====[end new stuff]=====
 
-call make_matrix (params,ndof,aelp,belp,icon,x,y,z, &
-                  kfix,mat(matel)%viscosity,mat(matel)%density, &
-                  mat(matel)%penalty,mat(matel)%expon, &
-                  mat(matel)%activationenergy,mat(matel)%expansion, &
-                  mat(matel)%diffusivity,mat(matel)%heat, &
-                  mat(matel)%plasticity_type, &
-                  mat(matel)%plasticity_parameters, &
-                  u,v,w,temp,pressure,strain,is_plastic_temp,nnode, &
-                  f,r0,s0,t0,rst,ileaves,eviscosityp,vbounded_temp,threadinfo, &
-                  weight)
+call make_matrix (params,ndof,aelp,belp,icon,x,y,z,kfix,mat(matel)%viscosity,  &
+                 mat(matel)%density,mat(matel)%penalty,mat(matel)%expon,       &
+                  mat(matel)%activationenergy,mat(matel)%expansion,            &
+                  mat(matel)%diffusivity,mat(matel)%heat,                      &
+                  mat(matel)%plasticity_type,mat(matel)%plasticity_parameters, &
+                  u,v,w,temp,pressure,strain,is_plastic_temp,nnode,f,r0,s0,t0, &
+                  rst,ileaves,eviscosityp,vbounded_temp,threadinfo,weight)
 
 ael=ael+aelp/(8.d0**level)
 bel=bel+belp/(8.d0**level)
@@ -205,49 +204,119 @@ 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
-    density=0.d0
-    penal=0.d0
-    expon=0.d0
-    activ=0.d0
-    expan=0.d0
-    diffusivity=0.d0
-    heat=0.d0
-    do i=1,nlsf
-       matel=params%materialn(i)
-       viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
-       density=density+vol_lsf(i)*mat(matel)%density
-       penal=penal+vol_lsf(i)*mat(matel)%penalty
-       expon=expon+vol_lsf(i)*mat(matel)%expon
-       activ=activ+vol_lsf(i)*mat(matel)%activationenergy
-       expan=expan+vol_lsf(i)*mat(matel)%expansion
-       diffusivity=diffusivity+vol_lsf(i)/mat(matel)%diffusivity
-       heat=heat+vol_lsf(i)*mat(matel)%heat
-    enddo
     vol_lsf0=1.d0-sum(vol_lsf)
-    matel=params%materialn(0)
-    viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
-    density=density+vol_lsf0*mat(matel)%density
-    penal=penal+vol_lsf0*mat(matel)%penalty
-    expon=expon+vol_lsf0*mat(matel)%expon
-    activ=activ+vol_lsf0*mat(matel)%activationenergy
-    expan=expan+vol_lsf0*mat(matel)%expansion
-    heat=heat+vol_lsf0*mat(matel)%heat
-    diffusivity=diffusivity+vol_lsf0/mat(matel)%diffusivity
-    diffusivity=1.d0/diffusivity! note that some properties add geometrically not algebraically
-
-    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
-       endif
-    enddo
+    
+    if (vol_lsf0.gt.eps) then                                                   ! Element contains void, always use divFEM
+      viscosity=0.d0
+      density=0.d0
+      penal=0.d0
+      expon=0.d0
+      activ=0.d0
+      expan=0.d0
+      diffusivity=0.d0
+      heat=0.d0
+      do i=1,nlsf
+        matel=params%materialn(i)
+        viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
+        density=density+vol_lsf(i)*mat(matel)%density
+        penal=penal+vol_lsf(i)*mat(matel)%penalty
+        expon=expon+vol_lsf(i)*mat(matel)%expon
+        activ=activ+vol_lsf(i)*mat(matel)%activationenergy
+        expan=expan+vol_lsf(i)*mat(matel)%expansion
+        diffusivity=diffusivity+vol_lsf(i)/mat(matel)%diffusivity
+        heat=heat+vol_lsf(i)*mat(matel)%heat
+      enddo
+      matel=params%materialn(0)
+      viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
+      density=density+vol_lsf0*mat(matel)%density
+      penal=penal+vol_lsf0*mat(matel)%penalty
+      expon=expon+vol_lsf0*mat(matel)%expon
+      activ=activ+vol_lsf0*mat(matel)%activationenergy
+      expan=expan+vol_lsf0*mat(matel)%expansion
+      heat=heat+vol_lsf0*mat(matel)%heat
+      diffusivity=diffusivity+vol_lsf0/mat(matel)%diffusivity
+      diffusivity=1.d0/diffusivity                                              ! Note that some properties add geometrically not algebraically
+    elseif (params%matrule.eq.1) then                                           ! Assign material properties of volumetric majority material
+      matel=params%materialn(sum(maxloc(vol_lsf)))
+      viscosity=mat(matel)%viscosity
+      density=mat(matel)%density
+      penal=mat(matel)%penalty
+      expon=mat(matel)%expon
+      activ=mat(matel)%activationenergy
+      expan=mat(matel)%expansion
+      diffusivity=mat(matel)%diffusivity
+      heat=mat(matel)%heat
+    elseif (params%matrule.eq.2) then                                           ! Assign material properties of volumetric minority material
+      matel=params%materialn(sum(minloc(vol_lsf)))
+      viscosity=mat(matel)%viscosity
+      density=mat(matel)%density
+      penal=mat(matel)%penalty
+      expon=mat(matel)%expon
+      activ=mat(matel)%activationenergy
+      expan=mat(matel)%expansion
+      diffusivity=mat(matel)%diffusivity
+      heat=mat(matel)%heat
+    else
+      if (params%matrule.ne.0) print *,'Invalid matrule value, using divFEM'
+      viscosity=0.d0
+      density=0.d0
+      penal=0.d0
+      expon=0.d0
+      activ=0.d0
+      expan=0.d0
+      diffusivity=0.d0
+      heat=0.d0
+      do i=1,nlsf
+        matel=params%materialn(i)
+        viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
+        density=density+vol_lsf(i)*mat(matel)%density
+        penal=penal+vol_lsf(i)*mat(matel)%penalty
+        expon=expon+vol_lsf(i)*mat(matel)%expon
+        activ=activ+vol_lsf(i)*mat(matel)%activationenergy
+        expan=expan+vol_lsf(i)*mat(matel)%expansion
+        diffusivity=diffusivity+vol_lsf(i)/mat(matel)%diffusivity
+        heat=heat+vol_lsf(i)*mat(matel)%heat
+      enddo
+      matel=params%materialn(0)
+      viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
+      density=density+vol_lsf0*mat(matel)%density
+      penal=penal+vol_lsf0*mat(matel)%penalty
+      expon=expon+vol_lsf0*mat(matel)%expon
+      activ=activ+vol_lsf0*mat(matel)%activationenergy
+      expan=expan+vol_lsf0*mat(matel)%expansion
+      heat=heat+vol_lsf0*mat(matel)%heat
+      diffusivity=diffusivity+vol_lsf0/mat(matel)%diffusivity
+      diffusivity=1.d0/diffusivity                                              ! Note that some properties add geometrically not algebraically
+    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
+    else
+      matel=params%materialn(0)
+      plasticity_type=mat(matel)%plasticity_type
+      plasticity_parameters=mat(matel)%plasticity_parameters
+    endif
+
+    !volmax=0.d0                                                                 ! Plasticity always uses a majority rule
+    !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
+    !   endif
+    !enddo
 
 !    if (plasticity_type/='No') then
 !       print *,vol_lsf
@@ -313,12 +382,10 @@ do kk=1,2
          t0p=t0+rst*(t(1)+1.d0)/2.d0
          rstp=rst/2.d0
          levelp=level+1
-         call make_cut (levelp,levelmax, &
-                           ndof,ael,bel,icon,x,y,z, &
-                           kfix,mat, &
-                           u,v,w,temp,pressure,strain,is_plastic,nnode, &
-                           f,lsfp,nlsf,r0p,s0p,t0p,rstp,jcut,ileaves,eviscosity,vbounded,params,threadinfo, &
-                           weightel)
+         call make_cut (levelp,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,u,v,w,&
+                       temp,pressure,strain,is_plastic,nnode,f,lsfp,nlsf,r0p,  &
+                       s0p,t0p,rstp,jcut,ileaves,eviscosity,vbounded,params,   &
+                       threadinfo,weightel)
       enddo
    enddo
 enddo
@@ -330,4 +397,4 @@ return
 end
 
 !-------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
+!-------------------------------------------------------------------------------
\ No newline at end of file