Page MenuHomec4science
Files F99397101

lm_static.py
No OneTemporary
Actions

File Metadata

Created
Fri, Jan 24, 04:50

lm_static.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()
dumper.dump()
# get the associated fields
u = dom.primal() # displacement
f = dom.gradient() # total forces
p = dom.getField(lm.position0)
# select boundary conditions
free_nodes = p[:, 0] > 1e-8
free_nodes *= p[:, 0]-1 < -1e-8
blocked_nodes = np.logical_not(free_nodes)
# displace the object
u[blocked_nodes, 0] = p[blocked_nodes, 0] * 1.
dom.updateGradient()
dumper.dump()
def objective_gradient(x):
x = x.reshape(u.shape)
u[free_nodes, :] = x[free_nodes, :]
dom.updateGradient()
f[blocked_nodes, :] = 0
return -f.flatten()
def objective_function(x):
x = x.reshape(u.shape)
u[free_nodes, :] = x[free_nodes, :]
dom.updateGradient()
return dom.potentialEnergy()
ret = minimize(objective_function, u.flatten(),
jac=objective_gradient)
print(ret)
print(ret.x.reshape(u.shape))
x = ret.x.reshape(u.shape)
u[free_nodes, :] = x[free_nodes, :]
dumper.dump()
# testing the result
assert (np.fabs(u[:, 0] - p[:, 0] * 1.) < 1e-10).any()

Event Timeline

c4science · Help