dumpers.py
f0619f7851b3master

dumpers.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/>.

			import numpy as np
			import scipy.interpolate as ip
			import tamaas as tm
			import os
			import sys

			try:
			import uvw
			except ImportError:
			print("uvw is not installed. You can install it with "
			+ "'pip install --user uvw'")
			sys.exit(1)


			def makePeriodic(data):
			data = np.append(data, np.expand_dims(data[:, 0, ...], axis=1), axis=1)
			data = np.append(data, np.expand_dims(data[:, :, 0, ...], axis=2), axis=2)
			return data


			def makeTractionPeriodic(data):
			data = np.append(data, np.expand_dims(data[0, ...], axis=0), axis=0)
			data = np.append(data, np.expand_dims(data[:, 0, ...], axis=1), axis=1)
			return data


			class UVWDumper(tm.ModelDumper):
			"""Dumper for elasto-plastic calculations"""
			def __init__(self, residual, filename, dump_numpys=False):
			tm.ModelDumper.__init__(self)
			self.residual = residual
			self.count = 0
			self.filename = filename
			self.dump_numpys = dump_numpys

			if not os.path.exists('paraview'):
			os.mkdir('paraview')
			if self.dump_numpys and not os.path.exists('numpys'):
			os.mkdir('numpys')

			def dump(self, model):
			self.uvw_dump(model)
			if self.dump_numpys:
			self.numpy_dump(model)
			self.count += 1

			def uvw_dump(self, model):
			"""Dump displacements, plastic deformations and stresses"""
			discretization = model.getDiscretization().copy()
			discretization[1] += 1
			discretization[2] += 1

			coordinates = [np.linspace(0, L, N)
			for L, N in zip(model.getSystemSize(),
			discretization)]
			# Correct order of coordinate dimensions
			dimension_indices = [1, 2, 0]

			# Creating rectilinear grid with correct order for components
			grid = uvw.RectilinearGrid(
			"paraview/{}_{:04d}.vtr".format(self.filename, self.count),
			(coordinates[i] for i in dimension_indices))

			# Point data arrays
			stress = uvw.DataArray(makePeriodic(self.residual.getStress()),
			dimension_indices, 'stress')
			plastic_strain = uvw.DataArray(
			makePeriodic(self.residual.getPlasticStrain()),
			dimension_indices, 'plastic_strain')
			displacement = uvw.DataArray(makePeriodic(model.getDisplacement()),
			dimension_indices, 'displacement')
			residual = uvw.DataArray(makePeriodic(self.residual.getVector()),
			dimension_indices, 'residual')

			grid.addPointData(stress)
			grid.addPointData(plastic_strain)
			grid.addPointData(displacement)
			grid.addPointData(residual)

			grid.write()

			def numpy_dump(self, model):
			normal_traction = makeTractionPeriodic(model.getTraction()[:, :, 2])
			normal_displacement = makeTractionPeriodic(
			model.getDisplacement()[0, :, :, 2])
			np.savez('numpys/{}_{:04d}.npz'.format(self.filename, self.count),
			traction=normal_traction,
			displacement=normal_displacement)


			class AkantuDumper(tm.ModelDumper):
			def __init__(self, ak_solver, filename, extra_dumps=False):
			tm.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)

tamaas/python/dumpers.pyf0619f7851b3master

dumpers.pyView Options

tamaas/python/dumpers.py
f0619f7851b3master

dumpers.py
View Options