diff --git a/python/tamaas/dumpers/__init__.py b/python/tamaas/dumpers/__init__.py
index 42d35b6..2d29c6e 100644
--- a/python/tamaas/dumpers/__init__.py
+++ b/python/tamaas/dumpers/__init__.py
@@ -1,341 +1,342 @@
# -*- mode:python; coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-2020 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 .
"""
Dumpers for the class tamaas.Model
"""
from __future__ import print_function
from sys import stderr
from os import makedirs
import os.path
import numpy as np
from .. import ModelDumper, model_type, mpi, type_traits
from ._helper import step_dump, directory_dump, local_slice, _is_surface_field
def _get_attributes(model):
"Get model attributes"
return {
'model_type': str(model.type),
- 'system_size': model.getSystemSize(),
- 'discretization': model.getDiscretization()
+ 'system_size': model.system_size,
+ 'discretization': model.global_shape,
}
+
class FieldDumper(ModelDumper):
"""Abstract dumper for python classes using fields"""
postfix = ''
extension = ''
name_format = "{basename}{postfix}.{extension}"
def __init__(self, basename, *fields, **kwargs):
"""Construct with desired fields"""
super(FieldDumper, self).__init__()
self.basename = basename
self.fields = list(fields)
self.all_fields = kwargs.get('all_fields', False)
def add_field(self, field):
"""Add another field to the dump"""
if field not in self.fields:
self.fields.append(field)
def dump_to_file(self, file_descriptor, model):
"""Dump to a file (name or handle)"""
def get_fields(self, model):
"""Get the desired fields"""
if not self.all_fields:
requested_fields = self.fields
else:
requested_fields = model.getFields()
return {field: model.getField(field) for field in requested_fields}
def dump(self, model):
"Dump model"
self.dump_to_file(self.file_path, model)
@property
def file_path(self):
"""Get the default filename"""
return self.name_format.format(basename=self.basename,
postfix=self.postfix,
extension=self.extension)
@directory_dump('numpys')
@step_dump
class NumpyDumper(FieldDumper):
"""Dumper to compressed numpy files"""
extension = 'npz'
def dump_to_file(self, file_descriptor, model):
"""Saving to compressed multi-field Numpy format"""
if mpi.size() > 1:
raise RuntimeError("NumpyDumper does not function "
"at all in parallel")
np.savez_compressed(file_descriptor, attrs=_get_attributes(model),
**self.get_fields(model))
try:
import h5py
@directory_dump('hdf5')
@step_dump
class H5Dumper(FieldDumper):
"""Dumper to HDF5 file format"""
extension = 'h5'
def dump_to_file(self, file_descriptor, model):
"""Saving to HDF5 with metadata about the model"""
# Setup for MPI
if not h5py.get_config().mpi and mpi.size() > 1:
raise RuntimeError("HDF5 does not have MPI support")
if mpi.size() > 1:
from mpi4py import MPI # noqa
mpi_args = dict(driver='mpio', comm=MPI.COMM_WORLD)
comp_args = {} # compression does not work in parallel
else:
mpi_args = {}
comp_args = dict(compression='gzip', compression_opts=7)
with h5py.File(file_descriptor, 'w', **mpi_args) as handle:
# Writing data
for name, field in self.get_fields(model).items():
shape = list(field.shape)
if mpi.size() > 1:
xdim = 0 if _is_surface_field(field, model) else 1
shape[xdim] = MPI.COMM_WORLD.allreduce(shape[xdim])
dset = handle.create_dataset(name, shape, field.dtype,
**comp_args)
dset[local_slice(field, model)] = field
# Writing metadata
for name, attr in _get_attributes(model).items():
handle.attrs[name] = attr
except ImportError:
pass
try:
import uvw # noqa
import uvw.parallel
@directory_dump('paraview')
@step_dump
class UVWDumper(FieldDumper):
"""Dumper to VTK files for elasto-plastic calculations"""
extension = 'vtr'
forbidden_fields = ['traction', 'gap']
def dump_to_file(self, file_descriptor, model):
"""Dump displacements, plastic deformations and stresses"""
if mpi.size() > 1:
raise RuntimeError("UVWDumper does not function "
"properly in parallel")
bdim = len(model.getBoundaryDiscretization())
# Local MPI size
lsize = model.getDiscretization()
gsize = mpi.global_shape(model.getBoundaryDiscretization())
gshape = gsize
if len(lsize) > bdim:
gshape = [model.getDiscretization()[0]] + gshape
# Space coordinates
coordinates = [np.linspace(0, L, N, endpoint=False)
for L, N in zip(model.getSystemSize(), gshape)]
# If model has subsurfce domain, z-coordinate is always first
dimension_indices = np.arange(bdim)
if len(lsize) > bdim:
dimension_indices += 1
dimension_indices = np.concatenate((dimension_indices, [0]))
coordinates[0] = \
np.linspace(0, model.getSystemSize()[0], gshape[0])
offset = np.zeros_like(dimension_indices)
offset[0] = mpi.local_offset(gsize)
rectgrid = uvw.RectilinearGrid if mpi.size() == 1 \
else uvw.parallel.PRectilinearGrid
# Creating rectilinear grid with correct order for components
coordlist = [coordinates[i][o:o+lsize[i]]
for i, o in zip(dimension_indices, offset)]
grid = rectgrid(
file_descriptor, coordlist,
compression=True,
offsets=offset,
)
# Iterator over fields we want to dump
fields_it = filter(lambda t: t[0] not in self.forbidden_fields,
self.get_fields(model).items())
# We make fields periodic for visualization
for name, field in fields_it:
array = uvw.DataArray(field, dimension_indices, name)
grid.addPointData(array)
grid.write()
@directory_dump('paraview')
class UVWGroupDumper(FieldDumper):
"Dumper to ParaViewData files"
extension = 'pvd'
def __init__(self, basename, *fields, **kwargs):
"""Construct with desired fields"""
super(UVWGroupDumper, self).__init__(basename, *fields, **kwargs)
subdir = os.path.join('paraview', basename + '-VTR')
if not os.path.exists(subdir):
makedirs(subdir)
self.uvw_dumper = UVWDumper(
os.path.join(basename + '-VTR', basename), *fields, **kwargs
)
self.group = uvw.ParaViewData(self.file_path, compression=True)
def dump_to_file(self, file_descriptor, model):
self.group.addFile(
self.uvw_dumper.file_path.replace('paraview/', ''),
timestep=self.uvw_dumper.count
)
self.group.write()
self.uvw_dumper.dump(model)
except ImportError as error:
print(error, file=stderr)
try:
from netCDF4 import Dataset
@directory_dump('netcdf')
class NetCDFDumper(FieldDumper):
"""Dumper to netCDF4 files"""
extension = "nc"
boundary_fields = ['traction', 'gap']
def _file_setup(self, grp, model):
grp.createDimension('frame', None)
# Local dimensions
model_dim = len(model.getDiscretization())
voigt_dim = type_traits[model.type].voigt
self._vec = grp.createDimension('spatial', model_dim)
self._tens = grp.createDimension('Voigt', voigt_dim)
self.model_info = model.getDiscretization(), model.type
# Create boundary dimensions
for label, size, length in zip(
"xy",
model.getBoundaryDiscretization(),
model.getBoundarySystemSize()
):
grp.createDimension(label, size)
coord = grp.createVariable(label, 'f8', (label,))
coord[:] = np.linspace(0, length, size, endpoint=False)
self._create_variables(
grp, model,
lambda f: _is_surface_field(f[1], model),
model.getBoundaryDiscretization(), "xy"
)
# Create volume dimension
if model.type in {model_type.volume_1d, model_type.volume_2d}:
size = model.getDiscretization()[0]
grp.createDimension("z", size)
coord = grp.createVariable("z", 'f8', ("z",))
coord[:] = np.linspace(0, model.getSystemSize()[0], size)
self._create_variables(
grp, model,
lambda f: not _is_surface_field(f[1], model),
model.getDiscretization(), "zxy"
)
self.has_setup = True
def dump_to_file(self, file_descriptor, model):
if mpi.size() > 1:
raise RuntimeError("NetCDFDumper does not function "
"properly in parallel")
mode = 'a' if os.path.isfile(file_descriptor) \
and getattr(self, 'has_setup', False) else 'w'
with Dataset(file_descriptor, mode,
format='NETCDF4_CLASSIC') as rootgrp:
if rootgrp.dimensions == {}:
self._file_setup(rootgrp, model)
if self.model_info != (model.getDiscretization(), model.type):
raise Exception("Unexpected model {}".format(model))
self._dump_generic(rootgrp, model)
def _create_variables(self, grp, model, predicate,
shape, dimensions):
field_dim = len(shape)
fields = list(filter(predicate, self.get_fields(model).items()))
dim_labels = list(dimensions[:field_dim])
for label, data in fields:
local_dim = []
# If we have an extra component
if data.ndim > field_dim:
if data.shape[-1] == self._tens.size:
local_dim = [self._tens.name]
elif data.shape[-1] == self._vec.size:
local_dim = [self._vec.name]
grp.createVariable(label, 'f8',
['frame'] + dim_labels + local_dim,
zlib=True)
def _dump_generic(self, grp, model):
fields = self.get_fields(model).items()
new_frame = grp.dimensions['frame'].size
for label, data in fields:
var = grp[label]
slice_in_global = (new_frame,) + local_slice(data, model)
var[slice_in_global] = np.array(data, dtype=np.double)
except ImportError:
pass
diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
index d5b569f..62b6942 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,311 +1,316 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "model.hh"
#include "adhesion_functional.hh"
#include "functional.hh"
#include "integral_operator.hh"
#include "model_dumper.hh"
#include "model_extensions.hh"
#include "model_factory.hh"
#include "numpy.hh"
#include "residual.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
using namespace py::literals;
/// Wrap functional classes
void wrapFunctionals(py::module& mod) {
py::class_,
functional::wrap::PyFunctional>
func(mod, "Functional");
func.def(py::init<>())
.def("computeF", &functional::Functional::computeF)
.def("computeGradF", &functional::Functional::computeGradF);
py::class_ adh(mod, "AdhesionFunctional",
func);
adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
&functional::AdhesionFunctional::setParameters)
// legacy wrapper code
.def("setParameters", &functional::AdhesionFunctional::setParameters);
py::class_(
mod, "ExponentialAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
py::class_(
mod, "MaugisAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
py::class_(
mod, "SquaredExponentialAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
}
template
std::unique_ptr> instanciateFromNumpy(numpy& num) {
std::unique_ptr> result = nullptr;
switch (num.ndim()) {
case 2:
result = std::make_unique>>(num);
return result;
case 3:
result = std::make_unique>>(num);
return result;
case 4:
result = std::make_unique>>(num);
return result;
default:
TAMAAS_EXCEPTION("instanciateFromNumpy expects the last dimension of numpy "
"array to be the number of components");
}
}
/// Wrap IntegralOperator
void wrapIntegralOperator(py::module& mod) {
py::class_(mod, "IntegralOperator")
.def("apply",
[](IntegralOperator& op, numpy input, numpy output) {
auto in = instanciateFromNumpy(input);
auto out = instanciateFromNumpy(output);
op.apply(*in, *out);
})
.def("updateFromModel", &IntegralOperator::updateFromModel)
.def("getModel", &IntegralOperator::getModel,
py::return_value_policy::reference)
.def("getKind", &IntegralOperator::getKind)
.def("getType", &IntegralOperator::getType);
py::enum_(mod, "integration_method")
.value("linear", integration_method::linear)
.value("cutoff", integration_method::cutoff);
}
/// Wrap BEEngine classes
void wrapBEEngine(py::module& mod) {
py::class_(mod, "BEEngine")
.def("solveNeumann", &BEEngine::solveNeumann)
.def("solveDirichlet", &BEEngine::solveDirichlet)
.def("getModel", &BEEngine::getModel, py::return_value_policy::reference)
.def("registerNeumann", &BEEngine::registerNeumann)
.def("registerDirichlet", &BEEngine::registerDirichlet);
}
template
void wrapModelTypeTrait(py::module& mod) {
using trait = model_type_traits;
py::class_(mod, trait::repr)
.def_property_readonly_static("dimension",
[](py::object) { return trait::dimension; })
.def_property_readonly_static(
"components", [](py::object) { return trait::components; })
.def_property_readonly_static(
"boundary_dimension",
[](py::object) { return trait::boundary_dimension; })
.def_property_readonly_static("voigt",
[](py::object) { return trait::voigt; })
.def_property_readonly_static("indices",
[](py::object) { return trait::indices; });
}
/// Wrap Models
void wrapModelClass(py::module& mod) {
py::enum_(mod, "model_type")
.value("basic_1d", model_type::basic_1d)
.value("basic_2d", model_type::basic_2d)
.value("surface_1d", model_type::surface_1d)
.value("surface_2d", model_type::surface_2d)
.value("volume_1d", model_type::volume_1d)
.value("volume_2d", model_type::volume_2d);
auto trait_mod = mod.def_submodule("_type_traits");
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
py::class_(mod, "Model")
.def_property_readonly("type", &Model::getType)
.def("setElasticity", &Model::setElasticity, "E"_a, "nu"_a)
.def_property("E", &Model::getYoungModulus, &Model::setYoungModulus)
.def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio)
.def("getHertzModulus", &Model::getHertzModulus)
.def("getYoungModulus", &Model::getYoungModulus)
.def("getShearModulus", &Model::getShearModulus)
.def("getPoissonRatio", &Model::getPoissonRatio)
.def("getTraction", (GridBase & (Model::*)()) & Model::getTraction,
py::return_value_policy::reference_internal)
.def("getDisplacement",
(GridBase & (Model::*)()) & Model::getDisplacement,
py::return_value_policy::reference_internal)
.def("getSystemSize", &Model::getSystemSize)
.def("getDiscretization", &Model::getDiscretization)
.def("getBoundarySystemSize", &Model::getBoundarySystemSize)
.def("getBoundaryDiscretization", &Model::getBoundaryDiscretization)
.def("solveNeumann", &Model::solveNeumann)
.def("solveDirichlet", &Model::solveDirichlet)
.def("dump", &Model::dump)
.def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>())
.def("setDumper",
[](Model&, const ModelDumper&) {
TAMAAS_EXCEPTION("setDumper() is not a member of Model; use "
"addDumper() instead");
})
.def("getBEEngine", &Model::getBEEngine,
py::return_value_policy::reference_internal)
.def("getIntegralOperator", &Model::getIntegralOperator,
"operator_name"_a, py::return_value_policy::reference_internal)
.def("registerField",
[](Model& m, std::string name, numpy field) {
auto f = instanciateFromNumpy(field);
m.registerField(name, std::move(f));
},
"field_name"_a, "field"_a, py::keep_alive<1, 3>())
.def("getField",
(GridBase & (Model::*)(const std::string&)) & Model::getField,
"field_name"_a, py::return_value_policy::reference_internal)
.def("getFields", &Model::getFields)
.def("applyElasticity",
[](Model& model, numpy stress, numpy strain) {
auto out = instanciateFromNumpy(stress);
auto in = instanciateFromNumpy(strain);
model.applyElasticity(*out, *in);
})
.def("__repr__",
[](const Model& m) {
std::stringstream ss;
ss << m;
return ss.str();
})
.def("__getitem__",
[](const Model& m, std::string key) -> const GridBase& {
try {
return m.getField(key);
} catch (std::out_of_range&) {
throw py::key_error();
}
},
py::return_value_policy::reference_internal)
.def("__contains__",
[](const Model& m, std::string key) {
const auto fields = m.getFields();
return std::find(fields.begin(), fields.end(), key) !=
fields.end();
},
py::keep_alive<0, 1>())
.def("__iter__",
[](const Model& m) {
const auto& fields = m.getFieldsMap();
return py::make_key_iterator(fields.cbegin(), fields.cend());
},
- py::keep_alive<0, 1>());
+ py::keep_alive<0, 1>())
+ .def_property_readonly("shape", &Model::getDiscretization)
+ .def_property_readonly("global_shape", &Model::getGlobalDiscretization)
+ .def_property_readonly("boundary_shape",
+ &Model::getBoundaryDiscretization)
+ .def_property_readonly("system_size", &Model::getSystemSize);
py::class_>(
mod, "ModelDumper")
.def(py::init<>())
.def("dump", &ModelDumper::dump, "model"_a)
.def("__lshift__",
[](ModelDumper& dumper, Model& model) { dumper << model; });
}
/// Wrap factory for models
void wrapModelFactory(py::module& mod) {
py::class_(mod, "ModelFactory")
.def_static("createModel", &ModelFactory::createModel, "model_type"_a,
"system_size"_a, "discretization"_a)
.def_static("createResidual", &ModelFactory::createResidual,
"model_type"_a, "sigma_y"_a, "hardening"_a = 0.)
.def_static("registerVolumeOperators",
&ModelFactory::registerVolumeOperators, "model"_a);
}
/// Wrap residual class
void wrapResidual(py::module& mod) {
// TODO adapt to n-dim
py::class_(mod, "Residual")
.def(py::init())
.def("computeResidual",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidual(*in);
})
.def("computeStress",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeStress(*in);
})
.def("updateState",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.updateState(*in);
})
.def("computeResidualDisplacement",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidualDisplacement(*in);
})
.def("applyTangent",
[](Residual& res, numpy& output, numpy& input,
numpy& current_strain_inc) {
auto out = instanciateFromNumpy(output);
auto in = instanciateFromNumpy(input);
auto inc = instanciateFromNumpy(current_strain_inc);
res.applyTangent(*out, *in, *inc);
},
"output"_a, "input"_a, "current_strain_increment"_a)
.def("getVector", &Residual::getVector,
py::return_value_policy::reference_internal)
.def("getPlasticStrain", &Residual::getPlasticStrain,
py::return_value_policy::reference_internal)
.def("getStress", &Residual::getStress,
py::return_value_policy::reference_internal)
.def("setIntegrationMethod", &Residual::setIntegrationMethod, "method"_a,
"cutoff"_a = 1e-12)
.def_property("yield_stress", &Residual::getYieldStress,
&Residual::setYieldStress)
.def_property("hardening_modulus", &Residual::getHardeningModulus,
&Residual::setHardeningModulus)
.def_property_readonly("model", &Residual::getModel);
}
void wrapModel(py::module& mod) {
wrapBEEngine(mod);
wrapModelClass(mod);
wrapModelFactory(mod);
wrapFunctionals(mod);
wrapResidual(mod);
wrapIntegralOperator(mod);
}
} // namespace wrap
} // namespace tamaas