			#!/usr/bin/python3

			from scipy.optimize import minimize
			import pylibmultiscale as lm
			import subprocess
			import numpy as np

			lm.loadModules()
			Dim = 3

			# create the managers
			lm.FilterManager.getManager(Dim)
			lm.ActionManager.getManager(Dim)
			lm.DomainMultiScale.getManager()

			# create a communicator and the groups
			comm = lm.Communicator.getCommunicator()
			comm.reset()
			nb_proc = comm.getNumberFreeProcs()
			comm.addGroup('md', nb_proc)
			fe_group = lm.Communicator.getGroup('fe')

			# 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.reset()
			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. # poisson's ratio
			]
			""")


			# create akantu domain
			dom = lm.DomainAkantu3('fe', fe_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
			epot = dom.potentialEnergy()
			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)
			print(p[:, 0]-1 < -1e-8)
			print(p[:, 0]-1)
			print(blocked_nodes)
			print(p[blocked_nodes, 0])

			# 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, :]
			return dom.potentialEnergy()


			ret = minimize(objective_function, u.flatten(),
			jac=objective_gradient, method='CG')
			print(ret)
			print(ret.x.reshape(u.shape))
			x = ret.x.reshape(u.shape)
			u[free_nodes, :] = x[free_nodes, :]
			dumper.dump()

libmultiscale/examples/single_scale_continuum_akantu/lm.py
3b768f5b4c8bmaster

lm.py
View Options

libmultiscale/examples/single_scale_continuum_akantu/lm.py3b768f5b4c8bmaster

lm.pyView Options

libmultiscale/examples/single_scale_continuum_akantu/lm.py
3b768f5b4c8bmaster

lm.py
View Options