Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/python/tamaas/dumpers/akantu_dumper.py
74b9233dc365features/netcdf-mpi

Last Change
Raw File

akantu_dumper.py
View Options

# @file
# @section LICENSE
#
# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program 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 Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
"""
Akantu related dumper
"""
from .. import ModelDumper
from ._helper import makeTractionPeriodic
import os
import numpy as np
import scipy.interpolate as ip
import uvw
__all__ = ["AkantuDumper"]
class AkantuDumper(ModelDumper):
"""Dumper used when plastic deformations are compute with Akantu"""
def __init__(self, ak_solver, filename, extra_dumps=False):
ModelDumper.__init__(self)
self.solver = ak_solver
self.model = self.solver.model
self.filename = filename
self.model.setBaseName(filename)
self.extra_dumps = extra_dumps
if not os.path.exists('paraview'):
os.mkdir('paraview')
if self.extra_dumps and not os.path.exists('numpys'):
os.mkdir('numpys')
fields = ["displacement",
"external_force",
"internal_force",
"inelastic_strain",
"iso_hardening",
"stress",
"blocked_dofs"]
for field in fields:
self.model.addDumpField(field)
self.bot_face = np.where(self.model.getMesh().getNodes()[:, 2] == 0)
self.count = 0
def dump(self, model):
self.akantu_dump(model)
if self.extra_dumps:
self.uvw_dump(model)
self.numpy_dump(model)
self.count += 1
def uvw_dump(self, model):
self.discretization = np.array(model.getDiscretization()[1:], copy=True)
self.discretization += 1
self.coordinates = [np.linspace(0, L, N)
for L, N in zip(model.getSystemSize()[1:],
self.discretization)]
# Correct order of coordinate dimensions
dimension_indices = range(2)
# Creating rectilinear grid with correct order for components
grid = uvw.RectilinearGrid(
"paraview/{}_traction_{:04d}.vtr".format(self.filename,
self.count),
(self.coordinates[i] for i in dimension_indices))
traction = uvw.DataArray(makeTractionPeriodic(model.getTraction()),
dimension_indices, 'traction')
grid.addPointData(traction)
grid.write()
def akantu_dump(self, model):
self.model.dump()
def numpy_dump(self, model):
arrays = {}
arrays['traction'] = makeTractionPeriodic(model.getTraction()[:, :, 2])
displacements = self.model.getDisplacement()[self.bot_face, 2]
bot_points = self.model.getMesh().getNodes()[self.bot_face, :-1]
bot_points = bot_points.reshape(
(self.discretization[0] * self.discretization[1], 2))
x, y = np.meshgrid(*self.coordinates, indexing='ij')
disp = ip.griddata(bot_points,
displacements.T, (x, y))
disp = disp.reshape(self.discretization)
arrays['displacement'] = disp - disp.mean()
np.savez('numpys/{}_{:04d}.npz'.format(self.filename,
self.count), **arrays)
c4science · Help