Page MenuHomec4science
Files F99439754

lm_dynamic.py
No OneTemporary
Actions

File Metadata

Created
Fri, Jan 24, 14:12

lm_dynamic.py
View Options

#!/usr/bin/python3
from scipy.optimize import minimize
import pylibmultiscale as lm
import subprocess
import numpy as np
lm.loadModules()
Dim = 3
lm.spatial_dimension.fset(Dim)
# create a communicator and the groups
comm = lm.Communicator.getCommunicator()
all_group = lm.Communicator.getGroup('all')
# create the geometries
cube = lm.Cube(Dim, 'fe_geom')
cube.params.bbox = [0, 1, 0, 1, 0, 1]
cube.init()
gManager = lm.GeometryManager.getManager()
gManager.addGeometry(cube)
# create Akantu mesh with gmsh
with open("mesh.geo", 'w') as f:
f.write("""
lc = 1;
Point(1) = {0., 0., 0., lc};
Point(2) = {0., 1., 0., lc};
Point(3) = {1., 1., 0., lc};
Point(4) = {1., 0., 0., lc};
Point(5) = {0., 0., 1., lc};
Point(6) = {0., 1., 1., lc};
Point(7) = {1., 1., 1., lc};
Point(8) = {1., 0., 1., lc};
Line(1) = {5, 8};
Line(2) = {8, 7};
Line(3) = {7, 6};
Line(4) = {6, 5};
Line(5) = {8, 4};
Line(6) = {4, 3};
Line(7) = {3, 7};
Line(8) = {4, 1};
Line(9) = {1, 2};
Line(10) = {2, 3};
Line(11) = {1, 5};
Line(12) = {6, 2};
Curve Loop(1) = {10, 7, 3, 12};
Plane Surface(1) = {1};
Curve Loop(2) = {6, -10, -9, -8};
Plane Surface(2) = {2};
Curve Loop(3) = {11, -4, 12, -9};
Plane Surface(3) = {3};
Curve Loop(4) = {2, 3, 4, 1};
Plane Surface(4) = {4};
Curve Loop(5) = {6, 7, -2, 5};
Plane Surface(5) = {5};
Curve Loop(6) = {8, 11, 1, 5};
Plane Surface(6) = {6};
Surface Loop(1) = {4, 5, 2, 1, 3, 6};
Volume(1) = {1};
""")
subprocess.call("gmsh -3 -o mesh.msh mesh.geo", shell=True)
# create the material file for Akantu
with open("material.dat", 'w') as f:
f.write("""
material elastic [
name = steel
rho = 1 # density
E = 1 # young's modulus
nu = 0.3 # poisson's ratio
]
""")
# create akantu domain
dom = lm.DomainAkantu3('fe', all_group)
dom.params.material_filename = "material.dat"
dom.params.mesh_filename = "mesh.msh"
# dom.params.pbc = [1, 0, 1]
dom.params.domain_geometry = "fe_geom"
dom.params.timestep = 1.
dom.init()
# create a visualization dumper with Paraview
dumper = lm.DumperParaview("aka")
dumper.params.input = dom
dumper.params.disp = True
dumper.params.vel = True
dumper.params.force = True
dumper.init()
# get the associated fields
u = dom.primal() # displacement
p = dom.getField(lm.position0)
acc = dom.acceleration()
# displace the object
u[:, 0] = p[:, 0] * .1
dom.updateGradient()
dumper.dump()
def performStep():
time_step = 1
u = dom.primal() # displacement
v = dom.primalTimeDerivative()
acc = dom.acceleration()
# predictor
v += 0.5 * time_step * acc
u += time_step * v
# treats syncing for parallelism
dom.enforceCompatibility()
# update force&acceleration
dom.updateGradient()
dom.updateAcceleration()
# corrector
v += 0.5 * time_step * acc
for i in range(0, 100):
performStep()
dumper.dump()

Event Timeline

c4science · Help