Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/python/tamaas/nonlinear_solvers/akantu_solver.py
797f71620eb8features/declarative

Last Change
Raw File

akantu_solver.py
View Options

"""
@author Lucas Frérot <lucas.frerot@epfl.ch>
@section LICENSE
Copyright (©) 2018 EPFL (Ecole Polytechnique Fédérale de
Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
Solides)
Tamaas is free software: you can redistribute it and/or modify it under the
terms of the GNU Lesser General Public License as published by the Free
Software Foundation, either version 3 of the License, or (at your option) any
later version.
Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
details.
You should have received a copy of the GNU Lesser General Public License
along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
"""
from .. import Residual, EPSolver
from ..dumpers._helper import makeTractionPeriodic
import akantu as ak
import numpy as np
import scipy.interpolate as ip
__all__ = ["AkantuSolver"]
class FromTraction:
"""Functor for Neumann boundary conditions"""
def __init__(self, interpolation):
self.interpolation = interpolation
def operator(self, quad_point, dual, coord, normals):
val = self.interpolation(*(coord[:-1]))
dual[2] = val[0]
class PoincareFEM_plast(Residual):
def __init__(self, tm_model, ak_model):
super().__init__(tm_model)
self.tm_model = tm_model
self.discretization = tm_model.getDiscretization().copy()
self.previous_pressure = np.zeros(self.discretization[1:])
self.residual_displacement = np.zeros(self.discretization[1:])
self.discretization[1] += 1
self.discretization[2] += 1
self.coordinates = [np.linspace(0, L, N)
for L, N in zip(tm_model.getSystemSize(),
self.discretization)]
self.model = ak_model
self.bot_face = np.where(self.model.getMesh().getNodes()[:, 2] == 0)
self.ids = \
{'inelastic_strain':
'previous_solid_mechanics_model:0:'
+ 'plastic_linear_isotropic_hardening:inelastic_strain',
'iso_hardening':
'previous_solid_mechanics_model:0:'
+ 'plastic_linear_isotropic_hardening:iso_hardening'}
self.savePreviousState()
def apply(self, field, out=None):
interpolation = ip.interp2d(self.coordinates[1],
self.coordinates[2],
makeTractionPeriodic(field))
functor = FromTraction(interpolation)
self.model.getExternalForce()[...] = 0
self.model.applyNeumannBC(functor, "bot")
self.model.solveStep()
displacements = self.model.getDisplacement()[self.bot_face, 2]
if out is not None:
# All these shenanigans because of node ordering
bot_points = self.model.getMesh().getNodes()[self.bot_face, :-1]
bot_points = bot_points.reshape(
(self.discretization[1] * self.discretization[2], 2))
x, y = np.meshgrid(*self.coordinates[1:], indexing='ij')
disp = ip.griddata(bot_points, displacements.T, (x, y))
out[...] = disp.reshape(self.discretization[1:])[:-1, :-1]
out -= out.mean()
def savePreviousState(self):
self.inelastic_strain = \
self.model.getMaterial(0).getArrayReal(
self.ids['inelastic_strain'], ak._tetrahedron_4).copy()
self.hardening = \
self.model.getMaterial(0).getArrayReal(
self.ids['iso_hardening'], ak._tetrahedron_4).copy()
self.previous_pressure[...] = self.tm_model.getTraction()[:, :, 2]
def restorePreviousState(self):
inelastic_strain = \
self.model.getMaterial(0).getArrayReal(
self.ids['inelastic_strain'], ak._tetrahedron_4)
hardening = \
self.model.getMaterial(0).getArrayReal(
self.ids['iso_hardening'], ak._tetrahedron_4)
inelastic_strain[...] = self.inelastic_strain
hardening[...] = self.hardening
def updateState(self, increment):
"""
Update the plastic deformation state
- Here we disregard the increment
"""
self.restorePreviousState()
self.apply(self.tm_model.getTraction()[:, :, 2])
self.apply(np.zeros_like(self.residual_displacement),
self.residual_displacement)
self.apply(self.tm_model.getTraction()[:, :, 2])
self.savePreviousState()
def getStrainIncrement(self):
return np.array()
class AkantuSolver(EPSolver):
"""Non-linear finite elements solver"""
def __init__(self, tm_model, material_file, mesh_name):
"""Initialize Akantu SolidMechanicsModel"""
ak.parseInput(material_file)
self.tm_model = tm_model
self.name = mesh_name
# Loading mesh
self.mesh = ak.Mesh(3)
self.mesh.read(mesh_name)
# Creating Akantu model
self.model = ak.SolidMechanicsModel(self.mesh)
self.model.initFull(_analysis_method=ak._static)
# Poincare operator/dummy_residual
self.poincare = PoincareFEM_plast(tm_model, self.model)
super().__init__(self.poincare)
# Setting boundary conditions
position = self.mesh.getNodes()
blocked_dofs = self.model.getBlockedDOFs()
epsilon = 1e-10
for i in range(self.mesh.getNbNodes()):
if np.abs(position[i, 2] - 1) < epsilon: # Blocking top surface
blocked_dofs[i, 2] = True
if np.abs(position[i, 0]) < epsilon:
blocked_dofs[i, 1] = True
if np.abs(position[i, 1]) < epsilon:
blocked_dofs[i, 0] = True
if np.abs(position[i, 0]) < epsilon \
or np.abs(position[i, 0] - 1.) < epsilon:
blocked_dofs[i, 0] = True
if np.abs(position[i, 1]) < epsilon \
or np.abs(position[i, 1] - 1.) < epsilon:
blocked_dofs[i, 1] = True
solver = self.model.getNonLinearSolver("static")
solver.set('max_iterations', 200)
solver.set('threshold', 1e-6)
solver.set('convergence_type', ak._scc_residual)
self.surface_correction = tm_model.getDisplacement()[0, :, :, 2]
self.previous_solution = np.zeros_like(self.surface_correction)
def solve(self):
"""Apply and release pressure to get the residual displacements"""
# this is the pressure increment
pressure = self.tm_model.getTraction()[:, :, 2]
self.poincare.restorePreviousState()
self.poincare.apply(pressure + self.poincare.previous_pressure)
residual_disp = np.zeros_like(self.surface_correction)
self.poincare.apply(np.zeros_like(pressure), residual_disp)
self.surface_correction[...] = \
residual_disp - self.poincare.residual_displacement
def reset(self):
self.previous_solution[...] = 0
c4science · Help