vesicles_in_rotating_boxΒΆ

../_images/vesicles_in_rotating_box.png
#/usr/bin/env python-mpacts
from mpacts.commands.monitors.progress import ProgressIndicator
from mpacts.contact.detectors.fixedlist import FixedListContactDetector
from mpacts.contact.detectors.multigrid import MultiGridContactDetector
from mpacts.contact.models.bending.commonedge import LinearBendingBetweenTriangles
from mpacts.contact.models.collision.linearpressure.linearpressure_rt import LinearPressureCoulombInt
from mpacts.contact.models.springs.elastic.elastic_basic import LinearSpringDamper
from mpacts.core.simulation import Simulation
from mpacts.geometrygenerators.meshobjects import IcoSphere
from mpacts.geometrygenerators.polyhedron import Connectivity
from mpacts.geometrygenerators.trianglegeometries import UnitIcoSphere
from mpacts.particles.particlecontainer import ParticleContainer
from mpacts.commands.time_evolution.rotatebody import RotateBodyCommand
import mpacts.geometrygenerators.quadgeometries as quadgen
import mpacts.particles as prt
import mpacts.commands.force.constraints.areaconservation as area
import mpacts.commands.force.constraints.volumeconservation as vol
import mpacts.geometrygenerators.transformations as transfo
from mpacts.core.valueproperties import Variable,VariableFunction

#-----------------------------------------------------------------------------------------------------------------------
dt = 6e-5
vesicle_positions = [(-0.25, 0, 0), (0.25, 0, 0)]
#-----------------------------------------------------------------------------------------------------------------------

mysim = Simulation("simulation", timestep=dt)
params = mysim('params')
endTime = Variable("endtime", params, value=2.5)
mass_vesicle = Variable("mass_vesicle", params, value=0.1)
r  = Variable("r_vesicle", params, value=0.2 )
subdivision = Variable("subdivision", params, value=3)
k  = Variable("k_spring", params, value=10)
c  = Variable('c_spring', params, value=2 )
kb = Variable("k_bend", params, value=0.025)
kv = Variable("bulk_modulus", params, value=100)
ka = Variable("ka", params, value=0.001)
kd = Variable("kd", params, value=0.001)
rmax = Variable("r_max", params, value=10 )

kc = Variable("k_contact", params, value=100)
c_vv = Variable("c_ves_ves", params, value=4)
c_vg = Variable("c_ves_groun", params, value=10)
ct = Variable("c_t", params, value=4)
ct_ground = Variable("c_t_ground", params, value=10)
mu = Variable('mu', params, value=0.25)

Lbox = Variable("L_box", params, value=1.2)
w = Variable("w", params, value=2 )

out_int = Variable("output_interval", params, value=0.015)
cd_ue  = Variable("cd_ue", params, value=20)
cd_kd = Variable("cd_kd", params, value=0.015)
g = Variable("g", params, value=(0,-9.81,0))

#-----------------------------------------------------------------------------------------------------------------------

box = ParticleContainer("box", prt.RigidBody.compose(
                       (prt.RigidQuad, "plane", "controlPoints")), mysim)

bx = box.add_particle()
(points, quads) = quadgen.create_box(Lbox.get_value(), Lbox.get_value(), Lbox.get_value())
bx.controlPoints.add_and_set_particles(xBF=points)
bx.plane.add_and_set_particles(vertexIndices=quads)

RotateBodyCommand("RotateBox", mysim, pc=box, particle_idx=0, w=w, axis = (0, 0, 1) )
box("plane").VTKWriterCmd(executeInterval=out_int)
vesicle = ParticleContainer("vesicle"
                           , prt.DeformableBody.compose( (prt.Node, "nodes")
                                                       , (prt.DeformableRoundedTriangle, "triangles", "nodes"))
                           , mysim)

vesicle("triangles").TriangleNormalsAndAreaCmd(x=vesicle('nodes')['x'])
vesicle("triangles").ComputeCurvatureCmd(maximal_radius=rmax, x=vesicle('nodes')['x'] )
vesicle("triangles").EncompassingSpheresCmd(x=vesicle('nodes')['x'])
area.AreaConservation(vesicle, vesicle("triangles"), ka=ka, kd=kd)
vol.VolumeConservation(vesicle, vesicle("triangles"), kv=kv)
vesicle("nodes").GravityCmd(g=g)
vesicle("nodes").TimeIntegration_ForwardEuler_Translation_Cmd()

writer = vesicle("triangles").VTKWriterCmd(executeInterval=out_int)
writer.select_all(True)
writer2 = vesicle.CSVWriterCmd(executeInterval=out_int)

spring = LinearSpringDamper( k=k
                           , c=c
                           , pc1=vesicle('nodes')
                           , pc2=vesicle('nodes'))

bending = LinearBendingBetweenTriangles( k_bend=kb
                                       , pc1=vesicle('triangles')
                                       , pc2=vesicle('triangles'))

cd_springs = FixedListContactDetector("CD_Springs", mysim, cmodel=spring)
cd_bending = FixedListContactDetector("CD_bending", mysim, cmodel=bending)

CM_ves_ves = LinearPressureCoulombInt( pc1=vesicle("triangles")
                                     , pc2=vesicle('triangles')
                                     , c_t=ct
                                     , mu=mu
                                     , k=kc
                                     , c=c_vv )

CM_ves_wall = LinearPressureCoulombInt( pc1=vesicle("triangles")
                                      , pc2=box('plane')
                                      , c_t=ct
                                      , mu=mu
                                      , k=kc
                                      , c=c_vg)

vesicleplaneCD = MultiGridContactDetector( "Vesicle-box-contact"
                                         , mysim
                                         , cmodel=CM_ves_wall
                                         , update_every=cd_ue
                                         , keep_distance=cd_kd)

vesiclevesicleCD = MultiGridContactDetector( "Vesicle-vesicle-contact"
                                           , mysim
                                           , cmodel=CM_ves_ves
                                           , update_every=cd_ue
                                           , keep_distance=cd_kd)

# Let's print something
ppcmd = ProgressIndicator("PrintProgress", mysim, print_interval=3)

#-----------------------------------------------------------------------------------------------------------------------

p = UnitIcoSphere(subdivision.get_value())
p.scaleVertices(r.get_value())

#p.optimize_bandwidth()#NOTE that if we enable this, the hash is not guaranteed anymore
#                       bewteen different versions of networkx, since the heuristic algorithm
#                       'reverse_cuthill_mckee_ordering' might result in a slightly different permutation.

con = Connectivity(p)
for point in vesicle_positions:
    ves = vesicle.add_particle()
    ves.nodes.add_and_set_particles( x=transfo.translate(p.vertices,point)
                                   , m=mass_vesicle.get_value()/p.nVertices() )
    ves.triangles.add_and_set_particles(vertexIndices=p.triangles)
    ves.add_connectivity(cd_springs, con.edgeCorners)
    ves.add_connectivity(cd_bending, con.edgeTriangles)

#print mysim.print_tree()
mysim.run_until( endTime.get_value())

#-----------------------------------------------------------------------------------------------------
#Verification that the output is identical to a reference value
from mpacts.tools.arrayhash import check_and_document_hash
check_and_document_hash( vesicle('nodes')['x'], 'fc157b3594d4840f5d786edbefc82e84', tolerance=1e-6)
#-----------------------------------------------------------------------------------------------------