diff --git a/cmake/check_builtin_expect.cc b/cmake/check_builtin_expect.cc
index c1e0fa4e3..3698dbd44 100644
--- a/cmake/check_builtin_expect.cc
+++ b/cmake/check_builtin_expect.cc
@@ -1,26 +1,26 @@
/**
* Copyright (©) 2022-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
#include
int main() {
if (__builtin_expect(true, 1)) {
std::cout << "has __builtin_expect" << std::endl;
}
}
diff --git a/cmake/check_constexpr_map.cc b/cmake/check_constexpr_map.cc
index 0c8c2101f..6e7e2bf9c 100644
--- a/cmake/check_constexpr_map.cc
+++ b/cmake/check_constexpr_map.cc
@@ -1,35 +1,35 @@
/**
* Copyright (©) 2022-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
#include
#include "aka_constexpr_map.hh"
int main() {
akantu::details::static_switch_dispatch(
std::tuple,
std::integral_constant>{},
[](auto && type) {
std::cout << std::decay_t::value << std::endl;
},
2, [](auto && /*type*/) { std::cout << "Default" << std::endl; },
std::make_index_sequence<2>{});
}
diff --git a/cmake/check_gnu_unlikely.cc b/cmake/check_gnu_unlikely.cc
index a567a1ad9..dbf35c049 100644
--- a/cmake/check_gnu_unlikely.cc
+++ b/cmake/check_gnu_unlikely.cc
@@ -1,26 +1,26 @@
/**
* Copyright (©) 2022-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
#include
int main() {
if (true) [[gnu::likely]] {
std::cout << "has __builtin_expect" << std::endl;
}
}
diff --git a/cmake/libakantu/v3/printers.py b/cmake/libakantu/v3/printers.py
index 652152366..2210e16fe 100755
--- a/cmake/libakantu/v3/printers.py
+++ b/cmake/libakantu/v3/printers.py
@@ -1,370 +1,356 @@
#!/usr/bin/env python3
__copyright__ = "Copyright (©) 2016-2023 EPFL (Ecole Polytechnique Fédérale" \
" de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \
" en Mécanique des Solides)"
__license__ = "LGPLv3"
-"""printers.py: gdb pretty printers."""
-
-__author__ = "Nicolas Richart"
-__credits__ = [
- "Nicolas Richart ",
-]
-__copyright__ = (
- "Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale"
- " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation"
- " en Mécanique des Solides)"
-)
-__license__ = "LGPLv3"
-
-
#
# Inspired from boost's pretty printers from
# Rüdiger Sonderfeld
# and from Pretty-printers for libstc++ from Free Software Foundation, Inc.
#
import gdb
import re
import itertools
try:
import gdb.printing
except ImportError:
raise RuntimeError("Your GDB is too old")
class AkantuPrinter(object):
regex = None
@classmethod
def supports(cls, typename):
# print('{0} ~= {1}'.format(typename, cls.regex), file=sys.stderr)
return cls.regex.search(typename)
@classmethod
def get_basic_type(cls, value):
"""Determines the type associated to a value"""
_type = value.type
# If it points to a reference, get the reference.
if _type.code == gdb.TYPE_CODE_REF:
_type = _type.target()
# Get the unqualified type, stripped of typedefs.
_type = _type.unqualified().strip_typedefs()
return _type.tag
__akantu_pretty_printers__ = \
gdb.printing.RegexpCollectionPrettyPrinter("libakantu-v3")
def register_pretty_printer(pretty_printer):
"""Registers a Pretty Printer"""
__akantu_pretty_printers__.add_printer(
pretty_printer.name, pretty_printer.regex, pretty_printer
)
return pretty_printer
@register_pretty_printer
class AkaArrayPrinter(AkantuPrinter):
"""Pretty printer for akantu::Array"""
regex = re.compile(
"^akantu::Array<(.*?), (true|false)>$|^akantu::ArrayDataLayer<(.*?), \(akantu::ArrayAllocationType\)1>$" # noqa
)
name = "akantu::Array"
def __init__(self, value):
self.typename = self.get_basic_type(value)
self.value = value
self.ptr = self.value["values"]
self.size = int(self.value["size_"])
self.nb_component = int(self.value["nb_component"])
def to_string(self):
m = self.regex.search(self.typename)
return "Array<{0}>({1}, {2}) stored at {3}".format(
m.group(1), self.size, self.nb_component, self.ptr
)
def display_hint(self):
return "array"
def children(self):
_ptr = self.ptr
for i in range(self.size):
for j in range(self.nb_component):
yield ("[{0}, {1}]".format(i, j), (_ptr + j).dereference())
_ptr = _ptr + self.nb_component
class RbtreeIterator(object):
"""
Turn an RB-tree-based container (std::map, std::set etc.) into
a Python iterable object.
"""
def __init__(self, rbtree):
self.size = rbtree["_M_t"]["_M_impl"]["_M_node_count"]
self.node = rbtree["_M_t"]["_M_impl"]["_M_header"]["_M_left"]
self.count = 0
def __iter__(self):
return self
def __len__(self):
return int(self.size)
def __next__(self):
if self.count == self.size:
raise StopIteration
result = self.node
self.count = self.count + 1
if self.count < self.size:
# Compute the next node.
node = self.node
if node.dereference()["_M_right"]:
node = node.dereference()["_M_right"]
while node.dereference()["_M_left"]:
node = node.dereference()["_M_left"]
else:
parent = node.dereference()["_M_parent"]
while node == parent.dereference()["_M_right"]:
node = parent
parent = parent.dereference()["_M_parent"]
if node.dereference()["_M_right"] != parent:
node = parent
self.node = node
return result
def get_value_from_aligned_membuf(buf, valtype):
"""Returns the value held in a __gnu_cxx::__aligned_membuf."""
return buf["_M_storage"].address.cast(valtype.pointer()).dereference()
def get_value_from_Rb_tree_node(node):
"""Returns the value held in an _Rb_tree_node<_Val>"""
try:
member = node.type.fields()[1].name
if member == "_M_value_field":
# C++03 implementation, node contains the value as a member
return node["_M_value_field"]
elif member == "_M_storage":
# C++11 implementation, node stores value in __aligned_membuf
valtype = node.type.template_argument(0)
return get_value_from_aligned_membuf(node["_M_storage"], valtype)
except: # noqa: E722
pass
raise ValueError("Unsupported implementation for %s" % str(node.type))
def find_type(orig, name):
typ = orig.strip_typedefs()
while True:
# Strip cv-qualifiers. PR 67440.
search = "%s::%s" % (typ.unqualified(), name)
try:
return gdb.lookup_type(search)
except RuntimeError:
pass
# The type was not found, so try the superclass. We only need
# to check the first superclass, so we don't bother with
# anything fancier here.
field = typ.fields()[0]
if not field.is_base_class:
raise ValueError("Cannot find type %s::%s" % (str(orig), name))
typ = field.type
@register_pretty_printer
class AkaElementTypeMapArrayPrinter(AkantuPrinter):
"""Pretty printer for akantu::ElementTypeMap>"""
regex = re.compile(
"^akantu::ElementTypeMap\*, akantu::(.*?)>$" # noqa
) # noqa: E501,W605
name = "akantu::ElementTypeMapArray"
# Turn an RbtreeIterator into a pretty-print iterator.
class _rb_iter(object):
def __init__(self, rbiter, type, ghost_type):
self.rbiter = rbiter
self.count = 0
self.type = type
self.ghost_type = ghost_type
def __iter__(self):
return self
def __next__(self):
if self.count % 2 == 0:
n = next(self.rbiter)
n = n.cast(self.type).dereference()
n = get_value_from_Rb_tree_node(n)
self.pair = n
item = "{0}:{1}".format(n["first"], self.ghost_type)
else:
item = self.pair["second"].dereference()
result = ("[{0}]".format(self.count), item)
self.count = self.count + 1
return result
def _iter(self, not_ghost, ghost, type):
iter_size = (len(not_ghost), len(ghost))
it = self._rb_iter(not_ghost, type, "_not_ghost")
for _ in range(iter_size[0] * 2):
yield next(it)
it = self._rb_iter(ghost, type, "_ghost")
for _ in range(iter_size[1] * 2):
yield next(it)
raise StopIteration
def __init__(self, value):
self.typename = self.get_basic_type(value)
self.value = value
self.data = self.value["data"]
self.ghost_data = self.value["ghost_data"]
def to_string(self):
m = self.regex.search(self.typename)
return "ElementTypMapArray<{0}> with " "{1} _not_ghost and {2} _ghost".format( # noqa
m.group(1),
len(RbtreeIterator(self.data)),
len(RbtreeIterator(self.ghost_data)))
def children(self):
# m = self.regex.search(self.typename)
rep_type = find_type(self.data.type, "_Rep_type")
node = find_type(rep_type, "_Link_type")
node = node.strip_typedefs()
return itertools.chain(
self._rb_iter(RbtreeIterator(self.data), node, "_not_ghost"),
self._rb_iter(RbtreeIterator(self.ghost_data), node, "_ghost"),
)
def display_hint(self):
return "map"
# @register_pretty_printer
class AkaTensorPrinter(AkantuPrinter):
"""Pretty printer for akantu::Tensor"""
regex = re.compile("^akantu::Tensor<(.*), +(.*), +(.*)>$")
name = "akantu::Tensor"
value = None
typename = ""
ptr = None
dims = []
ndims = 0
def pretty_print(self):
def ij2str(i, j, m):
return "{0}".format((self.ptr + m * j + i).dereference())
def line(i, m, n):
return "[{0}]".format(
", ".join((ij2str(i, j, m) for j in range(n))))
m = int(self.dims[0])
if self.ndims == 1:
n = 1
else:
n = int(self.dims[1])
return "[{0}]".format(", ".join(line(i, m, n) for i in range(m)))
def __init__(self, value):
self.typename = self.get_basic_type(value)
self.value = value
self.ptr = self.value["values"]
self.dims = self.value["n"]
def children(self):
yield ("values", self.pretty_print())
yield ("wrapped", self.value["wrapped"])
@register_pretty_printer
class AkaVectorPrinter(AkaTensorPrinter):
"""Pretty printer for akantu::Vector"""
regex = re.compile("^akantu::Vector<(.*)>$")
name = "akantu::Vector"
n = 0
ptr = 0x0
def __init__(self, value):
super(AkaVectorPrinter, self).__init__(value)
self.ndims = 1
def to_string(self):
m = self.regex.search(self.typename)
return "Vector<{0}>({1}) [{2}]".format(
m.group(1), int(self.dims[0]), str(self.ptr)
)
@register_pretty_printer
class AkaMatrixPrinter(AkaTensorPrinter):
"""Pretty printer for akantu::Matrix."""
regex = re.compile("^akantu::Matrix<(.*)>$")
name = "akantu::Matrix"
def __init__(self, value):
super(AkaMatrixPrinter, self).__init__(value)
self.ndims = 2
def to_string(self):
m = self.regex.search(self.typename)
return "Matrix<%s>(%d, %d) [%s]" % (
m.group(1),
int(self.dims[0]),
int(self.dims[1]),
str(self.ptr),
)
@register_pretty_printer
class AkaElementPrinter(AkantuPrinter):
"""Pretty printer for akantu::Element."""
regex = re.compile("^akantu::Element$")
name = "akantu::Element"
def __init__(self, value):
self.typename = self.get_basic_type(value)
self.value = value
self.element = self.value["element"]
self.eltype = self.value["type"]
self.ghost_type = self.value["ghost_type"]
self._ek_not_defined = gdb.parse_and_eval("akantu::_not_defined")
self._casper = gdb.parse_and_eval("akantu::_casper")
self._max_uint = gdb.parse_and_eval("(akantu::Int) -1")
def to_string(self):
if (
(self.element == self._max_uint)
and (self.eltype == self._ek_not_defined)
and (self.ghost_type == self._casper)
):
return "ElementNull"
return "Element({0}, {1}, {2})".format(
self.element, self.eltype, self.ghost_type
)
def register_akantu_printers(obj):
"""Register Akantu Pretty Printers."""
gdb.printing.register_pretty_printer(obj, __akantu_pretty_printers__)
diff --git a/cmake/material_lister.cc b/cmake/material_lister.cc
index 8dbc28da3..1b8cc7c91 100644
--- a/cmake/material_lister.cc
+++ b/cmake/material_lister.cc
@@ -1,43 +1,43 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "material_list.hh"
#include
#include
int main(__attribute__((unused)) int argc,
__attribute__((unused)) char * argv[]) {
#define PRINT_OUT_OPTIONS(r, data, i, elem) \
<< ":" << BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 0, elem))
#define PRINT_OUT(r, data, elem) \
std::cout << BOOST_PP_STRINGIZE(BOOST_PP_ARRAY_ELEM(0, elem)) << ":" \
<< BOOST_PP_STRINGIZE(BOOST_PP_ARRAY_ELEM(1, elem)) BOOST_PP_IF( \
BOOST_PP_EQUAL(3, BOOST_PP_ARRAY_SIZE(elem)), \
BOOST_PP_SEQ_FOR_EACH_I(PRINT_OUT_OPTIONS, _, \
BOOST_PP_ARRAY_ELEM(2, elem)), ) \
<< std::endl;
BOOST_PP_SEQ_FOR_EACH(PRINT_OUT, "toto", AKANTU_MATERIAL_LIST);
return 0;
}
diff --git a/cmake/semver.py b/cmake/semver.py
index d39f08480..19c3c06c5 100644
--- a/cmake/semver.py
+++ b/cmake/semver.py
@@ -1,240 +1,239 @@
#!/usr/bin/env python3
__copyright__ = "Copyright (©) 2022-2023 EPFL (Ecole Polytechnique Fédérale" \
" de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \
" en Mécanique des Solides)"
__license__ = "LGPLv3"
-"""module to get the VERSION from a git repository."""
+
import os
import re
import sys
import subprocess
def run_git_command(args):
"""Run git commands and capture outputs."""
git_dir = os.path.realpath(
os.path.join(os.path.dirname(__file__),
os.pardir))
cmd = ["git"] + args
p = subprocess.Popen(
cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=git_dir
)
stdout = p.communicate()[0].strip().decode()
if p.returncode != 0:
return None, p.returncode
return stdout, p.returncode
def _split_git_describe(describe):
describe_mo = re.search(
r"^(?P.+)"
r"-(?P\d+)"
r"-g(?P[0-9a-f]+)"
r"(-(?Pdirty))?$",
describe,
)
if describe_mo:
pieces = {}
pieces["tag"] = describe_mo.group("tag")
# distance: number of commits since tag
pieces["distance"] = int(describe_mo.group("distance"))
# commit: short hex revision ID
pieces["short"] = describe_mo.group("short")
if describe_mo.group("dirty"):
pieces["dirty"] = describe_mo.group("dirty")
return pieces
return None
def _eprint(*args, **kwargs):
print(*args, file=sys.stderr, **kwargs)
semver_re = re.compile(
r"^(?P0|[1-9]\d*)"
r"(\.(?P0|[1-9]\d*))?"
r"(\.(?P0|[1-9]\d*))?"
r"(?:-(?P(?:0|[1-9]\d*|\d*[a-zA-Z-][0-9a-zA-Z-]*)"
r"(?:\.(?:0|[1-9]\d*|\d*[a-zA-Z-][0-9a-zA-Z-]*))*))?"
r"(?:\+(?P[0-9a-zA-Z-]+(?:\.[0-9a-zA-Z-]+)*))?$"
)
def _parse_semver(version):
semver_mo = semver_re.search(version)
if not semver_mo:
return {}
pieces = {}
for p in ["major", "minor", "patch", "prerelease", "build"]:
if semver_mo.group(p):
pieces[p] = semver_mo.group(p)
return pieces
def get_git_version():
"""Get the version from the git repository."""
out, rc = run_git_command(["rev-parse", "--git-dir"])
if rc != 0:
return None
git_describe, rc = run_git_command(
["describe", "--tags", "--dirty", "--always", "--match", "v*"]
)
_eprint(f"git describe {git_describe}")
if "-g" in git_describe:
# TAG-DISTANCE-gHEX
pieces = _split_git_describe(git_describe)
else:
# tag only or no tag and hash
pieces = {"tag": git_describe}
# major.minor.patch-prerelease+build
if not pieces or ("tag" not in pieces):
return None
semver_mo = semver_re.search(pieces["tag"][1:])
if not semver_mo:
return pieces
for p in ["major", "minor", "patch", "prerelease", "build"]:
if semver_mo.group(p):
pieces[p] = semver_mo.group(p)
return pieces
def get_git_attributes_version():
"""Get the version from the attributes set a `git archive`."""
file_dir = os.path.dirname(os.path.realpath(os.path.abspath(__file__)))
attributes = None
pieces = None
with open(os.path.join(file_dir, "git_info")) as fh:
describe_re = re.compile(r"describe: ([^$].*[^$])")
for line in fh:
mo = describe_re.search(line)
if mo:
attributes = mo.group(1)
break
if attributes:
pieces = _split_git_describe(attributes)
return pieces
def get_ci_version():
"""Get extra information from CI context."""
pieces = None
if "CI_AKANTU_INSTALL_PREFIX" not in os.environ:
return None
ci_akantu_install_prefix = os.environ["CI_AKANTU_INSTALL_PREFIX"]
akantu_dir = os.path.join(ci_akantu_install_prefix,
"lib", "cmake", "Akantu")
cmake_config = os.path.join(akantu_dir, "AkantuConfig.cmake")
if not os.path.exists(cmake_config):
return None
version = None
with open(cmake_config, "r") as fh:
version_re = re.compile(r"^set\(AKANTU_VERSION (.*)\)$")
for line in fh:
version_mo = version_re.search(line)
if version_mo:
version = version_mo.group(1)
break
if not version:
return None
pieces = _parse_semver(version)
return pieces
def get_version_file():
"""Get the version directly from the VERSION file."""
version_path = os.path.join(
os.path.realpath(
os.path.abspath(os.path.join(os.path.dirname(__file__),
os.path.pardir))
),
"VERSION",
)
if not os.path.exists(version_path):
return None
version = None
with open(version_path, "r") as fh:
version = fh.readline()
if not version:
return None
pieces = _parse_semver(version)
return pieces
def get_version():
"""Combine all the version determination functions."""
pieces = None
if not pieces:
pieces = get_ci_version()
_eprint(f"pieces from ci_version: {pieces}")
if not pieces:
pieces = get_git_version()
_eprint(f"pieces from git_version: {pieces}")
if not pieces:
pieces = get_git_attributes_version()
_eprint(f"pieces from git attributes: {pieces}")
if not pieces:
pieces = get_version_file()
_eprint(f"pieces from version file: {pieces}")
if not pieces:
raise Exception("No version could be determined")
semver_build = []
if "build" in pieces:
semver_build = [pieces["build"]]
if "distance" in pieces:
semver_build.extend([str(pieces["distance"]), "g" + pieces["short"]])
if "dirty" in pieces and pieces["dirty"]:
semver_build.append(pieces["dirty"])
if semver_build:
pieces["build_part"] = "+" + ".".join(semver_build)
else:
pieces["build_part"] = ""
if "prerelease" in pieces:
pieces["prerelease"] = "-" + pieces["prerelease"]
else:
pieces["prerelease"] = ""
semver = "{major}.{minor}.{patch}{prerelease}{build_part}".format(**pieces)
if "CI_MERGE_REQUEST_ID" in os.environ:
semver = "{}.mr{}".format(semver, os.environ["CI_MERGE_REQUEST_ID"])
return semver
if __name__ == "__main__":
print(get_version())
-
diff --git a/examples/boundary_conditions/predefined_bc/predefined_bc.cc b/examples/boundary_conditions/predefined_bc/predefined_bc.cc
index 3f0dea6ff..b1ac57964 100644
--- a/examples/boundary_conditions/predefined_bc/predefined_bc.cc
+++ b/examples/boundary_conditions/predefined_bc/predefined_bc.cc
@@ -1,52 +1,52 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
Mesh mesh(2);
mesh.read("square.msh");
// model initialization
SolidMechanicsModel model(mesh);
model.initFull();
// Dirichlet boundary conditions
model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "Fixed_x");
model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "Fixed_y");
// output in a paraview file
model.setBaseName("plate");
model.addDumpFieldVector("displacement");
model.addDumpField("blocked_dofs");
model.addDumpField("external_force");
model.dump();
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/boundary_conditions/python_user_defined_bc/python_user_defined_bc.cc b/examples/boundary_conditions/python_user_defined_bc/python_user_defined_bc.cc
index 9b58320bc..2b20b484b 100644
--- a/examples/boundary_conditions/python_user_defined_bc/python_user_defined_bc.cc
+++ b/examples/boundary_conditions/python_user_defined_bc/python_user_defined_bc.cc
@@ -1,86 +1,86 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
using namespace akantu;
class PYBIND11_EXPORT SineBoundary : public BC::Dirichlet::DirichletFunctor {
public:
SineBoundary(Real amplitude, Real phase) {
py_module = py::module::import("boundary_condition");
py_sin_boundary = py_module.attr("SinBoundary")(amplitude, phase);
}
public:
inline void operator()(__attribute__((unused)) UInt node,
Vector & flags, Vector & primal,
const Vector & coord) const {
py_sin_boundary.attr("compute")(primal, coord, flags);
}
protected:
py::object py_sin_boundary;
py::module py_module;
};
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
py::scoped_interpreter guard{};
Int spatial_dimension = 2;
Mesh mesh(spatial_dimension);
mesh.read("fine_mesh.msh");
SolidMechanicsModel model(mesh);
/// model initialization
model.initFull();
/// boundary conditions
Vector traction(2, 0.2);
SineBoundary sin_boundary(.2, 10.);
model.applyBC(sin_boundary, "Fixed_x");
model.applyBC(BC::Dirichlet::FixedValue(0., _y), "Fixed_y");
model.applyBC(BC::Neumann::FromTraction(traction), "Traction");
// output a paraview file with the boundary conditions
model.setBaseName("plate");
model.addDumpFieldVector("displacement");
model.addDumpFieldVector("external_force");
model.addDumpField("blocked_dofs");
model.dump();
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/boundary_conditions/user_defined_bc/user_defined_bc.cc b/examples/boundary_conditions/user_defined_bc/user_defined_bc.cc
index 51b576307..e805b5f61 100644
--- a/examples/boundary_conditions/user_defined_bc/user_defined_bc.cc
+++ b/examples/boundary_conditions/user_defined_bc/user_defined_bc.cc
@@ -1,78 +1,78 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
class SineBoundary : public BC::Dirichlet::DirichletFunctor {
public:
SineBoundary(Real amp, Real phase, BC::Axis ax = _x)
: DirichletFunctor(ax), amplitude(amp), phase(phase) {}
public:
inline void operator()(__attribute__((unused)) UInt node,
Vector & flags, Vector & primal,
const Vector & coord) const {
DIRICHLET_SANITY_CHECK;
flags(axis) = true;
primal(axis) = -amplitude * std::sin(phase * coord(1));
}
protected:
Real amplitude;
Real phase;
};
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
Int spatial_dimension = 2;
Mesh mesh(spatial_dimension);
mesh.read("fine_mesh.msh");
SolidMechanicsModel model(mesh);
/// model initialization
model.initFull();
/// boundary conditions
Vector traction{.2, .2};
model.applyBC(SineBoundary(.2, 10., _x), "Fixed_x");
model.applyBC(BC::Dirichlet::FixedValue(0., _y), "Fixed_y");
model.applyBC(BC::Neumann::FromTraction(traction), "Traction");
// output a paraview file with the boundary conditions
model.setBaseName("plate");
model.addDumpFieldVector("displacement");
model.addDumpFieldVector("external_force");
model.addDumpField("blocked_dofs");
model.dump();
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/cohesive_element/cohesive_extrinsic/cohesive_extrinsic.cc b/examples/cohesive_element/cohesive_extrinsic/cohesive_extrinsic.cc
index 2704a87d9..5c455e103 100644
--- a/examples/cohesive_element/cohesive_extrinsic/cohesive_extrinsic.cc
+++ b/examples/cohesive_element/cohesive_extrinsic/cohesive_extrinsic.cc
@@ -1,118 +1,118 @@
/**
* Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model_cohesive.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
const Int spatial_dimension = 2;
const Int max_steps = 1000;
Mesh mesh(spatial_dimension);
mesh.read("triangle.msh");
SolidMechanicsModelCohesive model(mesh);
/// model initialization
model.initFull(_analysis_method = _explicit_lumped_mass,
_is_extrinsic = true);
Real time_step = model.getStableTimeStep() * 0.05;
model.setTimeStep(time_step);
std::cout << "Time step: " << time_step << std::endl;
CohesiveElementInserter & inserter = model.getElementInserter();
inserter.setLimit(_y, 0.30, 0.20);
model.updateAutomaticInsertion();
Array & position = mesh.getNodes();
Array & velocity = model.getVelocity();
Array & boundary = model.getBlockedDOFs();
Array & displacement = model.getDisplacement();
Int nb_nodes = mesh.getNbNodes();
/// boundary conditions
for (Int n = 0; n < nb_nodes; ++n) {
if (position(n, 1) > 0.99 || position(n, 1) < -0.99) {
boundary(n, 1) = true;
}
if (position(n, 0) > 0.99 || position(n, 0) < -0.99) {
boundary(n, 0) = true;
}
}
model.setBaseName("extrinsic");
model.addDumpFieldVector("displacement");
model.addDumpField("velocity");
model.addDumpField("acceleration");
model.addDumpField("internal_force");
model.addDumpField("stress");
model.addDumpField("grad_u");
model.dump();
/// initial conditions
Real loading_rate = 0.5;
Real disp_update = loading_rate * time_step;
for (Int n = 0; n < nb_nodes; ++n) {
velocity(n, 1) = loading_rate * position(n, 1);
}
/// Main loop
for (Int s = 1; s <= max_steps; ++s) {
/// update displacement on extreme nodes
for (Int n = 0; n < nb_nodes; ++n) {
if (position(n, 1) > 0.99 || position(n, 1) < -0.99) {
displacement(n, 1) += disp_update * position(n, 1);
}
}
model.checkCohesiveStress();
model.solveStep();
if (s % 10 == 0) {
model.dump();
std::cout << "passing step " << s << "/" << max_steps << std::endl;
}
}
Real Ed = model.getEnergy("dissipated");
Real Edt = 200 * std::sqrt(2);
std::cout << Ed << " " << Edt << std::endl;
if (Ed < Edt * 0.999 || Ed > Edt * 1.001 || std::isnan(Ed)) {
std::cout << "The dissipated energy is incorrect" << std::endl;
return EXIT_FAILURE;
}
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/cohesive_element/cohesive_extrinsic_ig_tg/cohesive_extrinsic_ig_tg.cc b/examples/cohesive_element/cohesive_extrinsic_ig_tg/cohesive_extrinsic_ig_tg.cc
index 421ff17ce..eb5110fb7 100644
--- a/examples/cohesive_element/cohesive_extrinsic_ig_tg/cohesive_extrinsic_ig_tg.cc
+++ b/examples/cohesive_element/cohesive_extrinsic_ig_tg/cohesive_extrinsic_ig_tg.cc
@@ -1,141 +1,141 @@
/**
* Copyright (©) 2014-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model_cohesive.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
/* -------------------------------------------------------------------------- */
class Velocity : public BC::Dirichlet::DirichletFunctor {
public:
explicit Velocity(SolidMechanicsModel & model, Real vel, BC::Axis ax = _x)
: DirichletFunctor(ax), model(model), vel(vel) {
disp = vel * model.getTimeStep();
}
public:
inline void operator()(UInt node, Vector & /*flags*/,
Vector & disp,
const Vector & coord) const {
Real sign = std::signbit(coord(axis)) ? -1. : 1.;
disp(axis) += sign * this->disp;
model.getVelocity()(node, axis) = sign * vel;
}
private:
SolidMechanicsModel & model;
Real vel, disp;
};
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
const Int spatial_dimension = 2;
const Int max_steps = 1000;
Mesh mesh(spatial_dimension);
mesh.read("square.msh");
SolidMechanicsModelCohesive model(mesh);
MaterialCohesiveRules rules{{{"btop", "bbottom"}, "tg_cohesive"},
{{"btop", "btop"}, "ig_cohesive"},
{{"bbottom", "bbottom"}, "ig_cohesive"}};
/// model initialization
auto cohesive_material_selector =
std::make_shared(model, rules);
auto bulk_material_selector =
std::make_shared>("physical_names",
model);
auto && current_selector = model.getMaterialSelector();
cohesive_material_selector->setFallback(bulk_material_selector);
bulk_material_selector->setFallback(current_selector);
model.setMaterialSelector(cohesive_material_selector);
model.initFull(_analysis_method = _explicit_lumped_mass,
_is_extrinsic = true);
Real time_step = model.getStableTimeStep() * 0.05;
model.setTimeStep(time_step);
std::cout << "Time step: " << time_step << std::endl;
model.assembleMassLumped();
auto & position = mesh.getNodes();
auto & velocity = model.getVelocity();
model.applyBC(BC::Dirichlet::FlagOnly(_y), "top");
model.applyBC(BC::Dirichlet::FlagOnly(_y), "bottom");
model.applyBC(BC::Dirichlet::FlagOnly(_x), "left");
model.applyBC(BC::Dirichlet::FlagOnly(_x), "right");
model.setBaseName("extrinsic");
model.addDumpFieldVector("displacement");
model.addDumpField("velocity");
model.addDumpField("acceleration");
model.addDumpField("internal_force");
model.addDumpField("stress");
model.addDumpField("grad_u");
model.addDumpField("material_index");
model.dump();
/// initial conditions
Real loading_rate = 0.1;
// bar_height = 2
Real VI = loading_rate * 2 * 0.5;
for (auto && data : zip(make_view(position, spatial_dimension),
make_view(velocity, spatial_dimension))) {
std::get<1>(data) = loading_rate * std::get<0>(data);
}
model.dump();
Velocity vely(model, VI, _y);
Velocity velx(model, VI, _x);
/// Main loop
for (Int s = 1; s <= max_steps; ++s) {
model.applyBC(vely, "top");
model.applyBC(vely, "bottom");
model.applyBC(velx, "left");
model.applyBC(velx, "right");
model.checkCohesiveStress();
model.solveStep();
if (s % 10 == 0) {
model.dump();
std::cout << "passing step " << s << "/" << max_steps << std::endl;
}
}
return 0;
}
diff --git a/examples/cohesive_element/cohesive_intrinsic/cohesive_intrinsic.cc b/examples/cohesive_element/cohesive_intrinsic/cohesive_intrinsic.cc
index fc673bd3b..0c760b0ea 100644
--- a/examples/cohesive_element/cohesive_intrinsic/cohesive_intrinsic.cc
+++ b/examples/cohesive_element/cohesive_intrinsic/cohesive_intrinsic.cc
@@ -1,127 +1,127 @@
/**
* Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "element_group.hh"
#include "mesh_iterators.hh"
#include "solid_mechanics_model_cohesive.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
static void updateDisplacement(SolidMechanicsModelCohesive &,
const ElementGroup &, Real);
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
const Int spatial_dimension = 2;
const Int max_steps = 350;
Mesh mesh(spatial_dimension);
mesh.read("triangle.msh");
SolidMechanicsModelCohesive model(mesh);
model.getElementInserter().setLimit(_x, -0.26, -0.24);
/// model initialization
model.initFull(_analysis_method = _explicit_lumped_mass,
_is_extrinsic = false);
Real time_step = model.getStableTimeStep() * 0.8;
model.setTimeStep(time_step);
std::cout << "Time step: " << time_step << std::endl;
Array & boundary = model.getBlockedDOFs();
Int nb_nodes = mesh.getNbNodes();
/// boundary conditions
for (Int dim = 0; dim < spatial_dimension; ++dim) {
for (Int n = 0; n < nb_nodes; ++n) {
boundary(n, dim) = true;
}
}
model.setBaseName("intrinsic");
model.addDumpFieldVector("displacement");
model.addDumpField("velocity");
model.addDumpField("acceleration");
model.addDumpField("stress");
model.addDumpField("grad_u");
model.addDumpField("external_force");
model.addDumpField("internal_force");
model.dump();
/// update displacement
auto && elements = mesh.createElementGroup("diplacement");
Vector barycenter(spatial_dimension);
for_each_element(
mesh,
[&](auto && el) {
mesh.getBarycenter(el, barycenter);
if (barycenter(_x) > -0.25)
elements.add(el, true);
},
_element_kind = _ek_regular);
Real increment = 0.01;
updateDisplacement(model, elements, increment);
/// Main loop
for (Int s = 1; s <= max_steps; ++s) {
model.solveStep();
updateDisplacement(model, elements, increment);
if (s % 1 == 0) {
model.dump();
std::cout << "passing step " << s << "/" << max_steps << std::endl;
}
}
Real Ed = model.getEnergy("dissipated");
Real Edt = 2 * sqrt(2);
std::cout << Ed << " " << Edt << std::endl;
if (Ed < Edt * 0.999 || Ed > Edt * 1.001 || std::isnan(Ed)) {
std::cout << "The dissipated energy is incorrect" << std::endl;
return EXIT_FAILURE;
}
finalize();
return EXIT_SUCCESS;
}
/* -------------------------------------------------------------------------- */
static void updateDisplacement(SolidMechanicsModelCohesive & model,
const ElementGroup & group, Real increment) {
Array & displacement = model.getDisplacement();
for (auto && node : group.getNodeGroup().getNodes()) {
displacement(node, 0) += increment;
}
}
diff --git a/examples/contact_mechanics/cohesive_contact_explicit_dynamic.cc b/examples/contact_mechanics/cohesive_contact_explicit_dynamic.cc
index 379def871..6f0ba6b25 100644
--- a/examples/contact_mechanics/cohesive_contact_explicit_dynamic.cc
+++ b/examples/contact_mechanics/cohesive_contact_explicit_dynamic.cc
@@ -1,148 +1,148 @@
/**
* Copyright (©) 2021-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "contact_mechanics_model.hh"
#include "coupler_solid_cohesive_contact.hh"
#include "solid_mechanics_model_cohesive.hh"
#include "surface_selector.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
const Int spatial_dimension = 2;
initialize("material-cohesive.dat", argc, argv);
Real time_step{0.};
Real time_factor = 0.1;
UInt max_steps = 25000;
Real max_displacement = 1e-3;
Mesh mesh(spatial_dimension);
mesh.read("cohesive-contact.msh");
CouplerSolidCohesiveContact coupler(mesh);
auto & solid = coupler.getSolidMechanicsModelCohesive();
auto & contact = coupler.getContactMechanicsModel();
auto && material_selector =
std::make_shared(solid);
material_selector->setFallback(solid.getMaterialSelector());
solid.setMaterialSelector(material_selector);
auto && surface_selector = std::make_shared(mesh);
contact.getContactDetector().setSurfaceSelector(surface_selector);
coupler.initFull(_analysis_method = _explicit_lumped_mass,
_is_extrinsic = true);
coupler.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "sides");
time_step = solid.getStableTimeStep();
time_step *= time_factor;
std::cout << "Time Step = " << time_step << "s (" << time_step << "s)"
<< std::endl;
coupler.setTimeStep(time_step);
coupler.setBaseName("cohesive-contact-explicit-dynamic");
coupler.addDumpFieldVector("displacement");
coupler.addDumpFieldVector("velocity");
coupler.addDumpFieldVector("normals");
coupler.addDumpField("blocked_dofs");
coupler.addDumpField("grad_u");
coupler.addDumpField("stress");
coupler.addDumpField("gaps");
coupler.addDumpField("areas");
auto & velocity = solid.getVelocity();
auto & gaps = contact.getGaps();
Real damping_ratio = 0.99;
auto increment = max_displacement / max_steps;
for (auto i : arange(max_steps)) {
coupler.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "loading");
coupler.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "fixed");
coupler.solveStep();
solid.checkCohesiveStress();
// damping velocities only along the contacting zone
for (auto && tuple : zip(gaps, make_view(velocity, spatial_dimension))) {
auto & gap = std::get<0>(tuple);
auto & vel = std::get<1>(tuple);
if (gap > 0) {
vel *= damping_ratio;
}
}
// dumping energies
if (i % 1000 == 0) {
Real epot = solid.getEnergy("potential");
Real ekin = solid.getEnergy("kinetic");
std::cerr << i << "," << i * increment << "," << epot << "," << ekin
<< "," << epot + ekin << "," << std::endl;
}
if (i % 1000 == 0) {
coupler.dump();
}
}
for (auto i : arange(max_steps)) {
solid.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "loading");
solid.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "fixed");
coupler.solveStep();
coupler.checkCohesiveStress();
// damping velocities only along the contacting zone
for (auto && tuple : zip(gaps, make_view(velocity, spatial_dimension))) {
auto & gap = std::get<0>(tuple);
auto & vel = std::get<1>(tuple);
if (gap > 0) {
vel *= damping_ratio;
}
}
// dumping energies
if (i % 1000 == 0) {
Real epot = solid.getEnergy("potential");
Real ekin = solid.getEnergy("kinetic");
std::cerr << i << "," << i * increment << "," << epot << "," << ekin
<< "," << epot + ekin << "," << std::endl;
}
if (i % 1000 == 0) {
coupler.dump();
}
}
}
diff --git a/examples/contact_mechanics/contact_explicit_dynamic.cc b/examples/contact_mechanics/contact_explicit_dynamic.cc
index 422527458..03e7fb776 100644
--- a/examples/contact_mechanics/contact_explicit_dynamic.cc
+++ b/examples/contact_mechanics/contact_explicit_dynamic.cc
@@ -1,121 +1,121 @@
/**
* Copyright (©) 2013-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "contact_mechanics_model.hh"
#include "coupler_solid_contact.hh"
#include "non_linear_solver.hh"
#include "solid_mechanics_model.hh"
#include "surface_selector.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
const Int spatial_dimension = 2;
initialize("material.dat", argc, argv);
Real time_step;
Real time_factor = 0.1;
UInt max_steps = 20000;
Real max_displacement = 5e-3;
Mesh mesh(spatial_dimension);
mesh.read("hertz.msh");
CouplerSolidContact coupler(mesh);
auto & solid = coupler.getSolidMechanicsModel();
auto & contact = coupler.getContactMechanicsModel();
auto && selector = std::make_shared>(
"physical_names", solid);
solid.setMaterialSelector(selector);
coupler.initFull(_analysis_method = _explicit_lumped_mass);
auto && surface_selector = std::make_shared(mesh);
contact.getContactDetector().setSurfaceSelector(surface_selector);
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "fixed");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "fixed");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "loading");
solid.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "symmetry");
time_step = solid.getStableTimeStep();
time_step *= time_factor;
std::cout << "Time Step = " << time_step << "s (" << time_step << "s)"
<< std::endl;
coupler.setTimeStep(time_step);
coupler.setBaseName("contact-explicit-dynamic");
coupler.addDumpFieldVector("displacement");
coupler.addDumpFieldVector("velocity");
coupler.addDumpFieldVector("normals");
coupler.addDumpFieldVector("contact_force");
coupler.addDumpFieldVector("external_force");
coupler.addDumpFieldVector("internal_force");
coupler.addDumpField("gaps");
coupler.addDumpField("areas");
coupler.addDumpField("blocked_dofs");
coupler.addDumpField("grad_u");
coupler.addDumpField("stress");
auto & velocity = solid.getVelocity();
auto & gaps = contact.getGaps();
Real damping_ratio = 0.99;
auto increment = max_displacement / max_steps;
for (auto i : arange(max_steps)) {
solid.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "loading");
coupler.solveStep();
// damping velocities only along the contacting zone
for (auto && tuple : zip(gaps, make_view(velocity, spatial_dimension))) {
auto & gap = std::get<0>(tuple);
auto & vel = std::get<1>(tuple);
if (gap > 0) {
vel *= damping_ratio;
}
}
// dumping energies
if (i % 1000 == 0) {
Real epot = solid.getEnergy("potential");
Real ekin = solid.getEnergy("kinetic");
std::cerr << i << "," << i * increment << "," << epot << "," << ekin
<< "," << epot + ekin << "," << std::endl;
}
if (i % 1000 == 0) {
coupler.dump();
}
}
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/contact_mechanics/contact_explicit_static.cc b/examples/contact_mechanics/contact_explicit_static.cc
index 7b6b33f8f..6ebd407cc 100644
--- a/examples/contact_mechanics/contact_explicit_static.cc
+++ b/examples/contact_mechanics/contact_explicit_static.cc
@@ -1,84 +1,84 @@
/**
* Copyright (©) 2013-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "contact_mechanics_model.hh"
#include "coupler_solid_contact.hh"
#include "non_linear_solver.hh"
#include "solid_mechanics_model.hh"
#include "surface_selector.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
const Int spatial_dimension = 2;
initialize("material.dat", argc, argv);
Mesh mesh(spatial_dimension);
mesh.read("hertz.msh");
CouplerSolidContact coupler(mesh);
auto & solid = coupler.getSolidMechanicsModel();
auto & contact = coupler.getContactMechanicsModel();
auto && selector = std::make_shared>(
"physical_names", solid);
solid.setMaterialSelector(selector);
coupler.initFull(_analysis_method = _static);
auto && surface_selector = std::make_shared(mesh);
contact.getContactDetector().setSurfaceSelector(surface_selector);
coupler.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "fixed");
coupler.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "fixed");
coupler.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "loading");
coupler.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "symmetry");
coupler.setBaseName("contact-explicit-static");
coupler.addDumpFieldVector("displacement");
coupler.addDumpFieldVector("normals");
coupler.addDumpFieldVector("contact_force");
coupler.addDumpFieldVector("external_force");
coupler.addDumpFieldVector("internal_force");
coupler.addDumpField("gaps");
coupler.addDumpField("areas");
coupler.addDumpField("blocked_dofs");
coupler.addDumpField("grad_u");
coupler.addDumpField("stress");
auto max_steps = 100u;
for (auto _ [[gnu::unused]] : arange(max_steps)) {
auto increment = 1e-4;
coupler.applyBC(BC::Dirichlet::IncrementValue(-increment, _y), "loading");
coupler.solveStep();
coupler.dump();
}
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/embedded/embedded.cc b/examples/embedded/embedded.cc
index 826f7e978..e699a5eb7 100644
--- a/examples/embedded/embedded.cc
+++ b/examples/embedded/embedded.cc
@@ -1,88 +1,88 @@
/**
* Copyright (©) 2015-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "embedded_interface_model.hh"
#include "non_linear_solver.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
const Int dim = 2;
// Loading the concrete mesh
Mesh mesh(dim);
mesh.read("concrete.msh");
// Loading the reinforcement mesh
Mesh reinforcement_mesh(dim, "reinforcement_mesh");
// Exception is raised because reinforcement
// mesh contains only segments, i.e. 1D elements
try {
reinforcement_mesh.read("reinforcement.msh");
} catch (debug::Exception & e) {
}
// Model creation
EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim);
model.initFull(EmbeddedInterfaceModelOptions(_static));
// Boundary conditions
model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "XBlocked");
model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "YBlocked");
Vector force(dim);
force(0) = 0.0;
force(1) = -1.0;
model.applyBC(BC::Neumann::FromTraction(force), "Force");
// Dumping the concrete
model.setBaseName("concrete");
model.addDumpFieldVector("displacement");
model.addDumpFieldVector("external_force");
model.addDumpFieldVector("internal_force");
model.addDumpFieldTensor("stress");
// Dumping the reinforcement
model.setBaseNameToDumper("reinforcement", "reinforcement");
model.addDumpFieldTensorToDumper(
"reinforcement", "stress_embedded"); // dumping stress in reinforcement
auto & solver = model.getNonLinearSolver();
solver.set("max_iterations", 1);
solver.set("threshold", 1e-6);
solver.set("convergence_type", SolveConvergenceCriteria::_residual);
model.solveStep();
// Dumping model
model.dump();
model.dump("reinforcement");
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/explicit/explicit_dynamic.cc b/examples/explicit/explicit_dynamic.cc
index 96f0d7a40..cbc98f7df 100644
--- a/examples/explicit/explicit_dynamic.cc
+++ b/examples/explicit/explicit_dynamic.cc
@@ -1,98 +1,98 @@
/**
* Copyright (©) 2014-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
const Int spatial_dimension = 3;
const Real pulse_width = 2.;
const Real A = 0.01;
Real time_step;
Real time_factor = 0.8;
Int max_steps = 1000;
Mesh mesh(spatial_dimension);
if (Communicator::getStaticCommunicator().whoAmI() == 0)
mesh.read("bar.msh");
SolidMechanicsModel model(mesh);
/// model initialization
model.initFull(_analysis_method = _explicit_lumped_mass);
time_step = model.getStableTimeStep();
std::cout << "Time Step = " << time_step * time_factor << "s (" << time_step
<< "s)" << std::endl;
time_step *= time_factor;
model.setTimeStep(time_step);
/// boundary and initial conditions
Array & displacement = model.getDisplacement();
const Array & nodes = mesh.getNodes();
for (Int n = 0; n < mesh.getNbNodes(); ++n) {
Real x = nodes(n) - 2;
// Sinus * Gaussian
Real L = pulse_width;
Real k = 0.1 * 2 * M_PI * 3 / L;
displacement(n) = A * sin(k * x) * exp(-(k * x) * (k * x) / (L * L));
}
std::ofstream energy;
energy.open("energy.csv");
energy << "id,rtime,epot,ekin,tot" << std::endl;
model.setBaseName("explicit_dynamic");
model.addDumpField("displacement");
model.addDumpField("velocity");
model.addDumpField("acceleration");
model.addDumpField("stress");
model.dump();
for (Int s = 1; s <= max_steps; ++s) {
model.solveStep();
Real epot = model.getEnergy("potential");
Real ekin = model.getEnergy("kinetic");
energy << s << "," << s * time_step << "," << epot << "," << ekin << ","
<< epot + ekin << "," << std::endl;
if (s % 10 == 0)
std::cout << "passing step " << s << "/" << max_steps << std::endl;
model.dump();
}
energy.close();
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/heat_transfer/heat_transfer_dynamics_2d.cc b/examples/heat_transfer/heat_transfer_dynamics_2d.cc
index 03fff9563..3f6cd9ad4 100644
--- a/examples/heat_transfer/heat_transfer_dynamics_2d.cc
+++ b/examples/heat_transfer/heat_transfer_dynamics_2d.cc
@@ -1,92 +1,92 @@
/**
* Copyright (©) 2011-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "heat_transfer_model.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
const Int spatial_dimension = 2;
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
// create mesh
Mesh mesh(spatial_dimension);
mesh.read("square.msh");
HeatTransferModel model(mesh);
// initialize everything
model.initFull();
// get stable time step
Real time_step = model.getStableTimeStep() * 0.8;
std::cout << "time step is:" << time_step << std::endl;
model.setTimeStep(time_step);
// boundary conditions
const Array & nodes = model.getFEEngine().getMesh().getNodes();
Array & boundary = model.getBlockedDOFs();
Array & temperature = model.getTemperature();
double length = 1.;
auto nb_nodes = model.getFEEngine().getMesh().getNbNodes();
for (Int i = 0; i < nb_nodes; ++i) {
temperature(i) = 100.;
Real dx = nodes(i, 0) - length / 4.;
Real dy = 0.0;
Real dz = 0.0;
if (spatial_dimension > 1)
dy = nodes(i, 1) - length / 4.;
if (spatial_dimension == 3)
dz = nodes(i, 2) - length / 4.;
Real d = sqrt(dx * dx + dy * dy + dz * dz);
if (d < 0.1) {
boundary(i) = true;
temperature(i) = 300.;
}
}
model.setBaseName("heat_transfer_square2d");
model.addDumpField("temperature");
model.addDumpField("temperature_rate");
model.addDumpField("internal_heat_rate");
// main loop
int max_steps = 15000;
for (int i = 0; i < max_steps; i++) {
model.solveStep();
if (i % 100 == 0)
model.dump();
if (i % 10 == 0)
std::cout << "Step " << i << "/" << max_steps << std::endl;
}
std::cout << "\n\n Stable Time Step is : " << time_step << "\n \n"
<< std::endl;
return 0;
}
diff --git a/examples/heat_transfer/heat_transfer_dynamics_3d.cc b/examples/heat_transfer/heat_transfer_dynamics_3d.cc
index 05b9ade59..155eda419 100644
--- a/examples/heat_transfer/heat_transfer_dynamics_3d.cc
+++ b/examples/heat_transfer/heat_transfer_dynamics_3d.cc
@@ -1,98 +1,98 @@
/**
* Copyright (©) 2011-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include
#include
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "heat_transfer_model.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
Int spatial_dimension = 3;
ElementType type = _tetrahedron_4;
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
Mesh mesh(spatial_dimension);
mesh.read("cube.msh");
HeatTransferModel model(mesh);
// initialize everything
model.initFull();
// get and set stable time step
Real time_step = model.getStableTimeStep() * 0.8;
std::cout << "Stable Time Step is : " << time_step / .8 << std::endl;
std::cout << "time step is:" << time_step << std::endl;
model.setTimeStep(time_step);
/// boundary conditions
const Array & nodes = mesh.getNodes();
Array & boundary = model.getBlockedDOFs();
Array & temperature = model.getTemperature();
auto nb_nodes = mesh.getNbNodes();
double length;
length = 1.;
for (Int i = 0; i < nb_nodes; ++i) {
temperature(i) = 100.;
// to insert a heat source
Real dx = nodes(i, 0) - length / 2.;
Real dy = nodes(i, 1) - length / 2.;
Real dz = nodes(i, 2) - length / 2.;
Real d = sqrt(dx * dx + dy * dy + dz * dz);
if (d < 0.1) {
boundary(i) = true;
temperature(i) = 300.;
}
}
model.setBaseName("heat_transfer_cube3d");
model.addDumpField("temperature");
model.addDumpField("temperature_rate");
model.addDumpField("internal_heat_rate");
// //for testing
int max_steps = 1000;
for (int i = 0; i < max_steps; i++) {
model.solveStep();
if (i % 100 == 0)
model.dump();
if (i % 10 == 0) {
std::cout << "Step " << i << "/" << max_steps << std::endl;
}
}
return 0;
}
diff --git a/examples/heat_transfer/heat_transfer_static_2d.cc b/examples/heat_transfer/heat_transfer_static_2d.cc
index 0c14c6f80..edf161db6 100644
--- a/examples/heat_transfer/heat_transfer_static_2d.cc
+++ b/examples/heat_transfer/heat_transfer_static_2d.cc
@@ -1,84 +1,84 @@
/**
* Copyright (©) 2011-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "heat_transfer_model.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
/* -------------------------------------------------------------------------- */
Int spatial_dimension = 2;
std::string base_name;
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
// create mesh
Mesh mesh(spatial_dimension);
mesh.read("square.msh");
HeatTransferModel model(mesh);
// initialize everything
model.initFull(_analysis_method = _static);
// boundary conditions
const Array & nodes = mesh.getNodes();
Array & blocked_dofs = model.getBlockedDOFs();
Array & temperature = model.getTemperature();
double length = 1.;
Int nb_nodes = nodes.size();
for (Int i = 0; i < nb_nodes; ++i) {
temperature(i) = 100.;
Real dx = nodes(i, 0);
Real dy = nodes(i, 1);
Vector dX = {dx, dy};
dX.array() -= length / 4.;
Real d = dX.norm();
if (d < 0.1) {
blocked_dofs(i) = true;
temperature(i) = 300.;
}
if (std::abs(dx) < 1e-4 || std::abs(dy) < 1e-4)
blocked_dofs(i) = true;
if (std::abs(dx - length) < 1e-4 || std::abs(dy - length) < 1e-4)
blocked_dofs(i) = true;
}
model.setBaseName("heat_transfer_static_2d");
model.addDumpField("temperature");
model.addDumpField("internal_heat_rate");
model.addDumpField("conductivity");
model.addDumpField("blocked_dofs");
model.dump();
model.solveStep();
model.dump();
return 0;
}
diff --git a/examples/implicit/implicit_dynamic.cc b/examples/implicit/implicit_dynamic.cc
index d15b1dda6..31ea8fbfb 100644
--- a/examples/implicit/implicit_dynamic.cc
+++ b/examples/implicit/implicit_dynamic.cc
@@ -1,140 +1,140 @@
/**
* Copyright (©) 2014-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "communicator.hh"
#include "non_linear_solver.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
/* -------------------------------------------------------------------------- */
const Real bar_length = 10.;
const Real bar_height = 1.;
const Real bar_depth = 1.;
const Real F = 5e3;
const Real L = bar_length;
const Real I = bar_depth * bar_height * bar_height * bar_height / 12.;
const Real E = 12e7;
const Real rho = 1000;
const Real m = rho * bar_height * bar_depth;
static Real w(UInt n) {
return n * n * M_PI * M_PI / (L * L) * sqrt(E * I / m);
}
static Real analytical_solution(Real time) {
return 2 * F * L * L * L / (pow(M_PI, 4) * E * I) *
((1. - cos(w(1) * time)) + (1. - cos(w(3) * time)) / 81. +
(1. - cos(w(5) * time)) / 625.);
}
const Int spatial_dimension = 2;
const Real time_step = 1e-4;
const Real max_time = 0.62;
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material_dynamic.dat", argc, argv);
Mesh mesh(spatial_dimension);
const auto & comm = Communicator::getStaticCommunicator();
Int prank = comm.whoAmI();
if (prank == 0)
mesh.read("beam.msh");
mesh.distribute();
SolidMechanicsModel model(mesh);
/// model initialization
model.initFull(_analysis_method = _implicit_dynamic);
Material & mat = model.getMaterial(0);
mat.setParam("E", E);
mat.setParam("rho", rho);
Array & force = model.getExternalForce();
Array & displacment = model.getDisplacement();
// boundary conditions
model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "blocked");
model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "blocked");
model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "roller");
const Array & trac_nodes =
mesh.getElementGroup("traction").getNodeGroup().getNodes();
bool dump_node = false;
if (trac_nodes.size() > 0 && mesh.isLocalOrMasterNode(trac_nodes(0))) {
force(trac_nodes(0), 1) = F;
dump_node = true;
}
// output setup
std::ofstream pos;
pos.open("position.csv");
if (!pos.good())
AKANTU_ERROR("Cannot open file \"position.csv\"");
pos << "id,time,position,solution" << std::endl;
model.setBaseName("dynamic");
model.addDumpFieldVector("displacement");
model.addDumpField("velocity");
model.addDumpField("acceleration");
model.addDumpField("external_force");
model.addDumpField("internal_force");
model.dump();
model.setTimeStep(time_step);
auto & solver = model.getNonLinearSolver();
solver.set("max_iterations", 100);
solver.set("threshold", 1e-12);
solver.set("convergence_type", SolveConvergenceCriteria::_solution);
/// time loop
Real time = 0.;
for (Int s = 1; time < max_time; ++s, time += time_step) {
if (prank == 0)
std::cout << s << "\r" << std::flush;
model.solveStep();
if (dump_node)
pos << s << "," << time << "," << displacment(trac_nodes(0), 1) << ","
<< analytical_solution(s * time_step) << std::endl;
if (s % 100 == 0)
model.dump();
}
std::cout << std::endl;
pos.close();
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/io/dumper/dumpable_interface.cc b/examples/io/dumper/dumpable_interface.cc
index a6df28ef3..c6dc21984 100644
--- a/examples/io/dumper/dumpable_interface.cc
+++ b/examples/io/dumper/dumpable_interface.cc
@@ -1,180 +1,180 @@
/**
* Copyright (©) 2015-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "element_group.hh"
#include "group_manager_inline_impl.hh"
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
#include "dumpable_inline_impl.hh"
#include "dumper_iohelper_paraview.hh"
/* -------------------------------------------------------------------------- */
#include "locomotive_tools.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
/*
In this example, we present dumpers::Dumpable which is an interface
for other classes who want to dump themselves.
Several classes of Akantu inheritate from Dumpable (Model, Mesh, ...).
In this example we reproduce the same tasks as example_dumper_low_level.cc
using this time Dumpable interface inherted by Mesh, NodeGroup and
ElementGroup.
It is then advised to read first example_dumper_low_level.cc.
*/
initialize(argc, argv);
// To start let us load the swiss train mesh and its mesh data information.
Int spatial_dimension = 2;
Mesh mesh(spatial_dimension);
mesh.read("swiss_train.msh");
/*
swiss_train.msh has the following physical groups that can be viewed with
GMSH:
"$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
6
2 1 "red"
2 2 "white"
2 3 "lwheel_1"
2 4 "lwheel_2"
2 5 "rwheel_2"
2 6 "rwheel_1"
$EndPhysicalNames
..."
*/
// Grouping nodes and elements belonging to train wheels (=four mesh data).
ElementGroup & wheels_elements =
mesh.createElementGroup("wheels", spatial_dimension);
wheels_elements.append(mesh.getElementGroup("lwheel_1"));
wheels_elements.append(mesh.getElementGroup("lwheel_2"));
wheels_elements.append(mesh.getElementGroup("rwheel_1"));
wheels_elements.append(mesh.getElementGroup("rwheel_2"));
const Array & lnode_1 =
(mesh.getElementGroup("lwheel_1")).getNodeGroup().getNodes();
const Array & lnode_2 =
(mesh.getElementGroup("lwheel_2")).getNodeGroup().getNodes();
const Array & rnode_1 =
(mesh.getElementGroup("rwheel_1")).getNodeGroup().getNodes();
const Array & rnode_2 =
(mesh.getElementGroup("rwheel_2")).getNodeGroup().getNodes();
Array & node = mesh.getNodes();
Int nb_nodes = mesh.getNbNodes();
// This time a 2D Array is created and a padding size of 3 is passed to
// NodalField in order to warp train deformation on Paraview.
Array displacement(nb_nodes, spatial_dimension);
// Create an ElementTypeMapArray for the colour
ElementTypeMapArray colour("colour");
colour.initialize(mesh, _with_nb_element = true);
/* ------------------------------------------------------------------------ */
/* Creating dumpers */
/* ------------------------------------------------------------------------ */
// Create dumper for the complete mesh and register it as default dumper.
auto && dumper =
std::make_shared("train", "./paraview/dumpable", false);
mesh.registerExternalDumper(dumper, "train", true);
mesh.addDumpMesh(mesh);
// The dumper for the filtered mesh can be directly taken from the
// ElementGroup and then registered as "wheels_elements" dumper.
auto && wheels = mesh.getGroupDumper("paraview_wheels", "wheels");
mesh.registerExternalDumper(wheels.shared_from_this(), "wheels");
mesh.setDirectoryToDumper("wheels", "./paraview/dumpable");
// Arrays and ElementTypeMapArrays can be added as external fields directly
mesh.addDumpFieldExternal("displacement", displacement);
ElementTypeMapArrayFilter filtered_colour(colour,
wheels_elements.getElements());
auto colour_field_wheel =
std::make_shared, true>>(
filtered_colour);
mesh.addDumpFieldExternal("color", colour_field_wheel);
mesh.addDumpFieldExternalToDumper("wheels", "displacement", displacement);
mesh.addDumpFieldExternalToDumper("wheels", "colour", colour);
// For some specific cases the Fields should be created, as when you want to
// pad an array
auto displacement_vector_field =
mesh.createNodalField(&displacement, "all", 3);
mesh.addDumpFieldExternal("displacement_as_paraview_vector",
displacement_vector_field);
mesh.addDumpFieldExternalToDumper("wheels", "displacement_as_paraview_vector",
displacement_vector_field);
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
// Fill the ElementTypeMapArray colour.
fillColour(mesh, colour);
/// Apply displacement and wheels rotation.
Real tot_displacement = 50.;
Real radius = 1.;
auto nb_steps = 100;
Real theta = tot_displacement / radius;
Vector l_center(spatial_dimension);
Vector r_center(spatial_dimension);
for (Int i = 0; i < spatial_dimension; ++i) {
l_center(i) = node(14, i);
r_center(i) = node(2, i);
}
for (Int i = 0; i < nb_steps; ++i) {
displacement.zero();
Real step_ratio = Real(i) / Real(nb_steps);
Real angle = step_ratio * theta;
applyRotation(l_center, angle, node, displacement, lnode_1);
applyRotation(l_center, angle, node, displacement, lnode_2);
applyRotation(r_center, angle, node, displacement, rnode_1);
applyRotation(r_center, angle, node, displacement, rnode_2);
for (Int j = 0; j < nb_nodes; ++j) {
displacement(j, _x) += step_ratio * tot_displacement;
}
/// Dump call is finally made through Dumpable interface.
mesh.dump();
mesh.dump("wheels");
}
finalize();
return 0;
}
diff --git a/examples/io/dumper/dumper_low_level.cc b/examples/io/dumper/dumper_low_level.cc
index 0ed354efd..44f3f8719 100644
--- a/examples/io/dumper/dumper_low_level.cc
+++ b/examples/io/dumper/dumper_low_level.cc
@@ -1,188 +1,188 @@
/**
* Copyright (©) 2015-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "element_group.hh"
#include "group_manager.hh"
#include "mesh.hh"
#include "dumper_elemental_field.hh"
#include "dumper_nodal_field.hh"
#include "dumper_iohelper_paraview.hh"
#include "locomotive_tools.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
/* This example aims at illustrating how to manipulate low-level methods of
DumperIOHelper.
The aims is to visualize a colorized moving train with Paraview */
initialize(argc, argv);
// To start let us load the swiss train mesh and its mesh data information.
// We aknowledge here a weel-known swiss industry for mesh donation.
Int spatial_dimension = 2;
Mesh mesh(spatial_dimension);
mesh.read("swiss_train.msh");
Array & nodes = mesh.getNodes();
Int nb_nodes = mesh.getNbNodes();
/* swiss_train.msh has the following physical groups that can be viewed with
GMSH:
"$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
6
2 1 "red"
2 2 "white"
2 3 "lwheel_1"
2 4 "lwheel_2"
2 5 "rwheel_2"
2 6 "rwheel_1"
$EndPhysicalNames
..."
*/
// Grouping nodes and elements belonging to train wheels (=four mesh data)
ElementGroup & wheels_elements =
mesh.createElementGroup("wheels", spatial_dimension);
wheels_elements.append(mesh.getElementGroup("lwheel_1"));
wheels_elements.append(mesh.getElementGroup("lwheel_2"));
wheels_elements.append(mesh.getElementGroup("rwheel_1"));
wheels_elements.append(mesh.getElementGroup("rwheel_2"));
const Array & lnode_1 =
(mesh.getElementGroup("lwheel_1")).getNodeGroup().getNodes();
const Array & lnode_2 =
(mesh.getElementGroup("lwheel_2")).getNodeGroup().getNodes();
const Array & rnode_1 =
(mesh.getElementGroup("rwheel_1")).getNodeGroup().getNodes();
const Array & rnode_2 =
(mesh.getElementGroup("rwheel_2")).getNodeGroup().getNodes();
/* Note this Array is constructed with three components in order to warp train
deformation on Paraview. A more appropriate way to do this is to set a
padding in the NodalField (See example_dumpable_interface.cc.) */
Array displacement(nb_nodes, 3);
// ElementalField are constructed with an ElementTypeMapArray.
ElementTypeMapArray colour;
colour.initialize(mesh, _with_nb_element = true);
/* ------------------------------------------------------------------------ */
/* Dumper creation */
/* ------------------------------------------------------------------------ */
// Creation of two DumperParaview. One for full mesh, one for a filtered
// mesh.
DumperParaview dumper("train", "./paraview/dumper", false);
DumperParaview wheels("wheels", "./paraview/dumper", false);
// Register the full mesh
dumper.registerMesh(mesh);
// Register a filtered mesh limited to nodes and elements from wheels groups
wheels.registerFilteredMesh(mesh, wheels_elements.getElements(),
wheels_elements.getNodeGroup().getNodes());
// Generate an output file of the two mesh registered.
dumper.dump();
wheels.dump();
/* At this stage no fields are attached to the two dumpers. To do so, a
dumpers::Field object has to be created. Several types of dumpers::Field
exist. In this example we present two of them.
NodalField to describe nodal displacements of our train.
ElementalField handling the color of our different part.
*/
// NodalField are constructed with an Array.
auto displ_field = std::make_shared>(displacement);
auto colour_field = std::make_shared>(colour);
// Register the freshly created fields to our dumper.
dumper.registerField("displacement", displ_field);
dumper.registerField("colour", colour_field);
// For the dumper wheels, fields have to be filtered at registration.
// Filtered NodalField can be simply registered by adding an Array
// listing the nodes.
auto displ_field_wheel = std::make_shared>(
displacement, 0, 0, &(wheels_elements.getNodeGroup().getNodes()));
wheels.registerField("displacement", displ_field_wheel);
// For the ElementalField, an ElementTypeMapArrayFilter has to be created.
ElementTypeMapArrayFilter filtered_colour(colour,
wheels_elements.getElements());
auto colour_field_wheel =
std::make_shared, true>>(
filtered_colour);
wheels.registerField("colour", colour_field_wheel);
/* ------------------------------------------------------------------------ */
// Now that the dumpers are created and the fields are associated, let's
// paint and move the train!
// Fill the ElementTypeMapArray colour according to mesh data information.
fillColour(mesh, colour);
// Apply displacement and wheels rotation.
Real tot_displacement = 50.;
Real radius = 1.;
Int nb_steps = 100;
Real theta = tot_displacement / radius;
Vector l_center(3);
Vector r_center(3);
for (Int i = 0; i < spatial_dimension; ++i) {
l_center(i) = nodes(14, i);
r_center(i) = nodes(2, i);
}
for (Int i = 0; i < nb_steps; ++i) {
displacement.zero();
Real angle = (Real)i / (Real)nb_steps * theta;
applyRotation(l_center, angle, nodes, displacement, lnode_1);
applyRotation(l_center, angle, nodes, displacement, lnode_2);
applyRotation(r_center, angle, nodes, displacement, rnode_1);
applyRotation(r_center, angle, nodes, displacement, rnode_2);
for (Int j = 0; j < nb_nodes; ++j) {
displacement(j, 0) += (Real)i / (Real)nb_steps * tot_displacement;
}
// Output results after each moving steps for main and wheel dumpers.
dumper.dump();
wheels.dump();
}
finalize();
return 0;
}
diff --git a/examples/io/dumper/locomotive_tools.cc b/examples/io/dumper/locomotive_tools.cc
index 556a0b3d0..986e83153 100644
--- a/examples/io/dumper/locomotive_tools.cc
+++ b/examples/io/dumper/locomotive_tools.cc
@@ -1,81 +1,81 @@
/**
* Copyright (©) 2015-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
#include "locomotive_tools.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
/* -------------------------------------------------------------------------- */
void applyRotation(const Vector & center, Real angle,
const Array & nodes, Array & displacement,
const Array & node_group) {
auto nodes_it = nodes.begin(nodes.getNbComponent());
auto disp_it = displacement.begin(center.size());
auto node_num_it = node_group.begin();
auto node_num_end = node_group.end();
Vector pos_rel(center.size());
for (; node_num_it != node_num_end; ++node_num_it) {
const Vector pos = nodes_it[*node_num_it];
for (Int i = 0; i < pos.size(); ++i)
pos_rel(i) = pos(i);
Vector dis = disp_it[*node_num_it];
pos_rel -= center;
Real radius = pos_rel.norm();
if (std::abs(radius) < Math::getTolerance())
continue;
Real phi_i = std::acos(pos_rel(_x) / radius);
if (pos_rel(_y) < 0)
phi_i *= -1;
dis(_x) = std::cos(phi_i - angle) * radius - pos_rel(_x);
dis(_y) = std::sin(phi_i - angle) * radius - pos_rel(_y);
}
}
/* -------------------------------------------------------------------------- */
void fillColour(const Mesh & mesh, ElementTypeMapArray & colour) {
const ElementTypeMapArray & phys_data =
mesh.getData("physical_names");
const Array & txt_colour = phys_data(_triangle_3);
auto & id_colour = colour(_triangle_3);
for (Int i = 0; i < txt_colour.size(); ++i) {
std::string phy_name = txt_colour(i);
if (phy_name == "red")
id_colour(i) = 3;
else if (phy_name == "white" || phy_name == "lwheel_1" ||
phy_name == "rwheel_1")
id_colour(i) = 2;
else
id_colour(i) = 1;
}
}
diff --git a/examples/io/dumper/locomotive_tools.hh b/examples/io/dumper/locomotive_tools.hh
index 336635670..0ed5f84f1 100644
--- a/examples/io/dumper/locomotive_tools.hh
+++ b/examples/io/dumper/locomotive_tools.hh
@@ -1,29 +1,29 @@
/**
* Copyright (©) 2012-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
void applyRotation(const ::akantu::Vector<::akantu::Real> ¢er,
::akantu::Real angle,
const ::akantu::Array<::akantu::Real> &nodes,
::akantu::Array<::akantu::Real> &displacement,
const ::akantu::Array<::akantu::Idx> &node_group);
void fillColour(const ::akantu::Mesh &mesh,
::akantu::ElementTypeMapArray<::akantu::Int> &colour);
diff --git a/examples/io/parser/example_parser.cc b/examples/io/parser/example_parser.cc
index fbcff4ae6..5c81c6402 100644
--- a/examples/io/parser/example_parser.cc
+++ b/examples/io/parser/example_parser.cc
@@ -1,77 +1,77 @@
/**
* Copyright (©) 2015-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "non_linear_solver.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
// Precise in initialize the name of the text input file to parse.
initialize("input_file.dat", argc, argv);
// Get the user ParserSection.
const ParserSection & usersect = getUserParser();
// getParameterValue() allows to extract data associated to a given parameter
// name
// and cast it in the desired type set as template paramter.
Mesh mesh(usersect.getParameterValue("spatial_dimension"));
mesh.read(usersect.getParameterValue("mesh_file"));
// getParameter() can be used with variable declaration (destination type is
// explicitly known).
Int max_iter = usersect.getParameter("max_nb_iterations");
Real precision = usersect.getParameter("precision");
// Following NumPy convention, data can be interpreted as Vector or Matrix
// structures.
Matrix eigen_stress = usersect.getParameter("stress");
SolidMechanicsModel model(mesh);
model.initFull(SolidMechanicsModelOptions(_static));
model.applyBC(BC::Dirichlet::FixedValue(0.0, _x),
usersect.getParameterValue("outter_crust"));
model.applyBC(BC::Dirichlet::FixedValue(0.0, _y),
usersect.getParameterValue("outter_crust"));
model.applyBC(BC::Neumann::FromStress(eigen_stress),
usersect.getParameterValue("inner_holes"));
model.setDirectory("./paraview");
model.setBaseName("swiss_cheese");
model.addDumpFieldVector("displacement");
auto & solver = model.getNonLinearSolver();
solver.set("max_iterations", max_iter);
solver.set("threshold", precision);
model.solveStep();
model.dump();
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/new_material/local_material_damage.cc b/examples/new_material/local_material_damage.cc
index 59ad26b41..4b133faa2 100644
--- a/examples/new_material/local_material_damage.cc
+++ b/examples/new_material/local_material_damage.cc
@@ -1,98 +1,98 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "local_material_damage.hh"
#include "solid_mechanics_model.hh"
namespace akantu {
/* -------------------------------------------------------------------------- */
LocalMaterialDamage::LocalMaterialDamage(SolidMechanicsModel & model,
const ID & id)
: Material(model, id), damage("damage", *this) {
AKANTU_DEBUG_IN();
this->registerParam("E", E, 0., _pat_parsable, "Young's modulus");
this->registerParam("nu", nu, 0.5, _pat_parsable, "Poisson's ratio");
this->registerParam("lambda", lambda, _pat_readable,
"First Lamé coefficient");
this->registerParam("mu", mu, _pat_readable, "Second Lamé coefficient");
this->registerParam("kapa", kpa, _pat_readable, "Bulk coefficient");
this->registerParam("Yd", Yd, 50., _pat_parsmod);
this->registerParam("Sd", Sd, 5000., _pat_parsmod);
damage.initialize(1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void LocalMaterialDamage::initMaterial() {
AKANTU_DEBUG_IN();
Material::initMaterial();
lambda = nu * E / ((1 + nu) * (1 - 2 * nu));
mu = E / (2 * (1 + nu));
kpa = lambda + 2. / 3. * mu;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void LocalMaterialDamage::computeStress(ElementType el_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
Real * dam = damage(el_type, ghost_type).data();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
computeStressOnQuad(grad_u, sigma, *dam);
++dam;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void LocalMaterialDamage::computePotentialEnergy(ElementType el_type) {
AKANTU_DEBUG_IN();
Real * epot = potential_energy(el_type).data();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost);
computePotentialEnergyOnQuad(grad_u, sigma, *epot);
epot++;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
static bool material_is_alocated_local_damage [[gnu::unused]] =
MaterialFactory::getInstance().registerAllocator(
"local_damage",
[](UInt, const ID &, SolidMechanicsModel & model,
const ID & id) -> std::unique_ptr {
return std::make_unique(model, id);
});
} // namespace akantu
diff --git a/examples/new_material/local_material_damage.hh b/examples/new_material/local_material_damage.hh
index ab3df80cf..a7943f1ff 100644
--- a/examples/new_material/local_material_damage.hh
+++ b/examples/new_material/local_material_damage.hh
@@ -1,109 +1,109 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "material.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_LOCAL_MATERIAL_DAMAGE_HH_
#define AKANTU_LOCAL_MATERIAL_DAMAGE_HH_
namespace akantu {
class LocalMaterialDamage : public Material {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
LocalMaterialDamage(SolidMechanicsModel &model, const ID &id = "");
virtual ~LocalMaterialDamage(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
void initMaterial() override;
/// constitutive law for all element of a type
void computeStress(ElementType el_type,
GhostType ghost_type = _not_ghost) override;
/// compute the potential energy for all elements
void computePotentialEnergy(ElementType el_type) override;
protected:
/// constitutive law for a given quadrature point
template
inline void computeStressOnQuad(Eigen::MatrixBase &grad_u,
Eigen::MatrixBase &sigma, Real &damage);
/// compute the potential energy for on element
template
inline void computePotentialEnergyOnQuad(Eigen::MatrixBase &grad_u,
Eigen::MatrixBase &sigma,
Real &epot);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// compute the celerity of the fastest wave in the material
inline Real getCelerity(const Element &element) const override;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real);
private:
/// the young modulus
Real E;
/// Poisson coefficient
Real nu;
/// First Lamé coefficient
Real lambda;
/// Second Lamé coefficient (shear modulus)
Real mu;
/// resistance to damage
Real Yd;
/// damage threshold
Real Sd;
/// Bulk modulus
Real kpa;
/// damage internal variable
InternalField damage;
};
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "local_material_damage_inline_impl.hh"
#endif /* AKANTU_LOCAL_MATERIAL_DAMAGE_HH_ */
diff --git a/examples/new_material/local_material_damage_inline_impl.hh b/examples/new_material/local_material_damage_inline_impl.hh
index ddda748cb..251295014 100644
--- a/examples/new_material/local_material_damage_inline_impl.hh
+++ b/examples/new_material/local_material_damage_inline_impl.hh
@@ -1,80 +1,80 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_LOCAL_MATERIAL_DAMAGE_INLINE_IMPL_HH_
#define AKANTU_LOCAL_MATERIAL_DAMAGE_INLINE_IMPL_HH_
namespace akantu {
/* -------------------------------------------------------------------------- */
template
inline void LocalMaterialDamage::computeStressOnQuad(
Eigen::MatrixBase &grad_u, Eigen::MatrixBase &sigma, Real &dam) {
Real trace = grad_u.trace();
/// \sigma_{ij} = \lambda * (\nabla u)_{kk} * \delta_{ij} + \mu * (\nabla
/// u_{ij} + \nabla u_{ji})
for (Int i = 0; i < spatial_dimension; ++i) {
for (Int j = 0; j < spatial_dimension; ++j) {
sigma(i, j) =
(i == j) * lambda * trace + mu * (grad_u(i, j) + grad_u(j, i));
}
}
Real Y = 0;
for (Int i = 0; i < spatial_dimension; ++i) {
for (Int j = 0; j < spatial_dimension; ++j) {
Y += sigma(i, j) * grad_u(i, j);
}
}
Y *= 0.5;
Real Fd = Y - Yd - Sd * dam;
if (Fd > 0)
dam = (Y - Yd) / Sd;
dam = std::min(dam, 1.);
sigma *= 1 - dam;
}
/* -------------------------------------------------------------------------- */
template
inline void LocalMaterialDamage::computePotentialEnergyOnQuad(
Eigen::MatrixBase &grad_u, Eigen::MatrixBase &sigma, Real &epot) {
epot = 0.;
for (Int i = 0, t = 0; i < spatial_dimension; ++i)
for (Int j = 0; j < spatial_dimension; ++j, ++t)
epot += sigma(i, j) * (grad_u(i, j) - (i == j));
epot *= .5;
}
/* -------------------------------------------------------------------------- */
inline Real LocalMaterialDamage::getCelerity(__attribute__((unused))
const Element &element) const {
// Here the fastest celerity is the push wave speed
return (std::sqrt((2 * mu + lambda) / rho));
}
} // namespace akantu
#endif /* AKANTU_LOCAL_MATERIAL_DAMAGE_INLINE_IMPL_HH_ */
diff --git a/examples/new_material/new_local_material.cc b/examples/new_material/new_local_material.cc
index 4fb84d413..a10970ccf 100644
--- a/examples/new_material/new_local_material.cc
+++ b/examples/new_material/new_local_material.cc
@@ -1,91 +1,91 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "local_material_damage.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
#define bar_length 10.
#define bar_height 4.
akantu::Real eps = 1e-10;
int main(int argc, char * argv[]) {
akantu::initialize("material.dat", argc, argv);
Int max_steps = 10000;
Real epot, ekin;
const Int spatial_dimension = 2;
Mesh mesh(spatial_dimension);
mesh.read("barre_trou.msh");
/// model creation
SolidMechanicsModel model(mesh);
/// model initialization
model.initFull(_analysis_method = _explicit_lumped_mass);
std::cout << model.getMaterial(0) << std::endl;
Real time_step = model.getStableTimeStep();
model.setTimeStep(time_step / 10.);
/// Dirichlet boundary conditions
model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "Fixed_x");
model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "Fixed_y");
// Neumann boundary condition
Matrix stress(2, 2);
stress.eye(3e2);
model.applyBC(BC::Neumann::FromStress(stress), "Traction");
model.setBaseName("local_material");
model.addDumpField("displacement");
model.addDumpField("velocity");
model.addDumpField("acceleration");
model.addDumpField("external_force");
model.addDumpField("internal_force");
model.addDumpField("grad_u");
model.addDumpField("stress");
model.addDumpField("damage");
model.dump();
for (Int s = 0; s < max_steps; ++s) {
model.solveStep();
epot = model.getEnergy("potential");
ekin = model.getEnergy("kinetic");
if (s % 100 == 0)
std::cout << s << " " << epot << " " << ekin << " " << epot + ekin
<< std::endl;
if (s % 1000 == 0)
model.dump();
}
akantu::finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc b/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc
index a8e4607e8..1873a7631 100644
--- a/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc
+++ b/examples/new_material/viscoelastic_maxwell/material_viscoelastic_maxwell_energies.cc
@@ -1,157 +1,157 @@
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
#include "material_viscoelastic_maxwell.hh"
#include "non_linear_solver.hh"
#include "solid_mechanics_model.hh"
#include "sparse_matrix.hh"
using namespace akantu;
/* -------------------------------------------------------------------------- */
/* Main */
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
akantu::initialize("material_viscoelastic_maxwell.dat", argc, argv);
// sim data
Real eps = 0.1;
const Int dim = 2;
Real sim_time = 100.;
Real T = 10.;
Mesh mesh(dim);
mesh.read("material_viscoelastic_maxwell_mesh.msh");
SolidMechanicsModel model(mesh);
/* ------------------------------------------------------------------------ */
/* Initialization */
/* ------------------------------------------------------------------------ */
model.initFull(_analysis_method = _static);
std::cout << model.getMaterial(0) << std::endl;
std::stringstream filename_sstr;
filename_sstr << "material_viscoelastic_maxwell_output.out";
std::ofstream output_data;
output_data.open(filename_sstr.str().c_str());
Material & mat = model.getMaterial(0);
Real time_step = 0.1;
const Array & coordinate = mesh.getNodes();
Array & displacement = model.getDisplacement();
Array & blocked = model.getBlockedDOFs();
/// Setting time step
model.setTimeStep(time_step);
model.setBaseName("dynamic");
model.addDumpFieldVector("displacement");
model.addDumpField("blocked_dofs");
model.addDumpField("external_force");
model.addDumpField("internal_force");
model.addDumpField("grad_u");
model.addDumpField("stress");
model.addDumpField("strain");
Int max_steps = sim_time / time_step + 1;
Real time = 0.;
auto & solver = model.getNonLinearSolver();
solver.set("max_iterations", 10);
solver.set("threshold", 1e-7);
solver.set("convergence_type", SolveConvergenceCriteria::_residual);
/* ------------------------------------------------------------------------ */
/* Main loop */
/* ------------------------------------------------------------------------ */
for (Int s = 0; s <= max_steps; ++s) {
std::cout << "Time Step = " << time_step << "s" << std::endl;
std::cout << "Time = " << time << std::endl;
// impose displacement
Real epsilon = 0;
if (time < T) {
epsilon = eps * time / T;
} else {
epsilon = eps;
}
for (auto && [coord, disp, block] :
zip(make_view(coordinate, dim), make_view(displacement, dim),
make_view(blocked, dim))) {
if (Math::are_float_equal(coord(_x), 0.0)) {
disp(_x) = 0;
disp(_y) = epsilon * coord(_y);
block.set(true);
} else if (Math::are_float_equal(coord(_y), 0.0)) {
disp(_x) = epsilon * coord(_x);
disp(_y) = 0;
block.set(true);
} else if (Math::are_float_equal(coord(_x), 0.001)) {
disp = epsilon * coord;
block.set(true);
} else if (Math::are_float_equal(coord(_y), 0.001)) {
disp = epsilon * coord;
block.set(true);
}
}
try {
model.solveStep();
} catch (debug::NLSNotConvergedException & e) {
std::cout << "Didn't converge after " << e.niter
<< " iterations. Error is " << e.error << std::endl;
return EXIT_FAILURE;
}
// for debugging
// auto int_force = model.getInternalForce();
// auto &K = model.getDOFManager().getMatrix("K");
// K.saveMatrix("K.mtx");
Int nb_iter = solver.get("nb_iterations");
std::cout << "Converged in " << nb_iter << " iterations" << std::endl;
model.dump();
Real epot = mat.getEnergy("potential");
Real edis = mat.getEnergy("dissipated");
Real work = mat.getEnergy("work");
// data output
output_data << s * time_step << " " << epsilon << " " << epot << " " << edis
<< " " << work << std::endl;
time += time_step;
}
output_data.close();
finalize();
}
diff --git a/examples/parallel/parallel_2d.cc b/examples/parallel/parallel_2d.cc
index c8e77f0b5..288abdb98 100644
--- a/examples/parallel/parallel_2d.cc
+++ b/examples/parallel/parallel_2d.cc
@@ -1,97 +1,97 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "communicator.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
Int spatial_dimension = 2;
Int max_steps = 10000;
Real time_factor = 0.8;
Real max_disp = 1e-6;
Mesh mesh(spatial_dimension);
const auto & comm = Communicator::getStaticCommunicator();
Int prank = comm.whoAmI();
if (prank == 0) {
// Read the mesh
mesh.read("square_2d.msh");
}
mesh.distribute();
SolidMechanicsModel model(mesh);
model.initFull();
if (prank == 0)
std::cout << model.getMaterial(0) << std::endl;
model.setBaseName("multi");
model.addDumpFieldVector("displacement");
model.addDumpFieldVector("velocity");
model.addDumpFieldVector("acceleration");
model.addDumpFieldTensor("stress");
model.addDumpFieldTensor("grad_u");
/// boundary conditions
Real eps = 1e-16;
const Array & pos = mesh.getNodes();
Array & disp = model.getDisplacement();
Array & boun = model.getBlockedDOFs();
Real left_side = mesh.getLowerBounds()(0);
Real right_side = mesh.getUpperBounds()(0);
for (Int i = 0; i < mesh.getNbNodes(); ++i) {
if (std::abs(pos(i, 0) - left_side) < eps) {
disp(i, 0) = max_disp;
boun(i, 0) = true;
}
if (std::abs(pos(i, 0) - right_side) < eps) {
disp(i, 0) = -max_disp;
boun(i, 0) = true;
}
}
Real time_step = model.getStableTimeStep() * time_factor;
std::cout << "Time Step = " << time_step << "s" << std::endl;
model.setTimeStep(time_step);
model.dump();
for (Int s = 1; s <= max_steps; ++s) {
model.solveStep();
if (s % 200 == 0)
model.dump();
if (prank == 0 && s % 100 == 0)
std::cout << "passing step " << s << "/" << max_steps << std::endl;
}
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/phase_field/phase_field_notch.cc b/examples/phase_field/phase_field_notch.cc
index 73d07c836..de2019ab3 100644
--- a/examples/phase_field/phase_field_notch.cc
+++ b/examples/phase_field/phase_field_notch.cc
@@ -1,109 +1,109 @@
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "coupler_solid_phasefield.hh"
#include "non_linear_solver.hh"
#include "phase_field_model.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
using clk = std::chrono::high_resolution_clock;
using second = std::chrono::duration;
using millisecond = std::chrono::duration;
const Int spatial_dimension = 2;
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material_notch.dat", argc, argv);
// create mesh
Mesh mesh(spatial_dimension);
mesh.read("square_notch.msh");
CouplerSolidPhaseField coupler(mesh);
auto & model = coupler.getSolidMechanicsModel();
auto & phase = coupler.getPhaseFieldModel();
model.initFull(_analysis_method = _static);
auto && mat_selector =
std::make_shared>("physical_names",
model);
model.setMaterialSelector(mat_selector);
auto && selector = std::make_shared>(
"physical_names", phase);
phase.setPhaseFieldSelector(selector);
phase.initFull(_analysis_method = _static);
model.applyBC(BC::Dirichlet::FixedValue(0., _y), "bottom");
model.applyBC(BC::Dirichlet::FixedValue(0., _x), "left");
model.setBaseName("phase_notch");
model.addDumpField("stress");
model.addDumpField("grad_u");
model.addDumpFieldVector("displacement");
model.addDumpField("damage");
model.dump();
auto nbSteps = 1500;
Real increment = 1e-5;
auto start_time = clk::now();
for (Int s = 1; s < nbSteps; ++s) {
if (s >= 500) {
increment = 1.e-6;
}
if (s % 10 == 0) {
constexpr char wheel[] = "/-\\|";
auto elapsed = clk::now() - start_time;
auto time_per_step = elapsed / s;
std::cout << "\r[" << wheel[(s / 10) % 4] << "] " << std::setw(5) << s
<< "/" << nbSteps << " (" << std::setprecision(2) << std::fixed
<< std::setw(8) << millisecond(time_per_step).count()
<< "ms/step - elapsed: " << std::setw(8)
<< second(elapsed).count() << "s - ETA: " << std::setw(8)
<< second((nbSteps - s) * time_per_step).count() << "s)"
<< std::string(' ', 20) << std::flush;
}
model.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "top");
coupler.solve();
if (s % 100 == 0) {
model.dump();
}
}
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/python/cohesive/custom_material.py b/examples/python/cohesive/custom_material.py
index 25c8ad4ae..92677463d 100755
--- a/examples/python/cohesive/custom_material.py
+++ b/examples/python/cohesive/custom_material.py
@@ -1,196 +1,194 @@
-__copyright__ = "Copyright (©) 2022-2023 EPFL (Ecole Polytechnique Fédérale" \
- " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \
- " en Mécanique des Solides)"
-__license__ = "LGPLv3"
-
-
#!/usr/bin/env python3
# pylint: disable=missing-module-docstring
# pylint: disable=missing-function-docstring
"""Example of custom cohesive material written in python."""
+__copyright__ = "Copyright (©) 2022-2023 EPFL (Ecole Polytechnique Fédérale" \
+ " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \
+ " en Mécanique des Solides)"
+__license__ = "LGPLv3"
+
import numpy as np
import akantu as aka
spatial_dimension = 2
class LinearCohesive(aka.MaterialCohesive):
"""Material Cohesive Linear."""
def __init__(self, model, _id):
"""Construct the material and register the parameters to the parser."""
super().__init__(model, _id)
super().registerParamReal(
"G_c", aka._pat_readable | aka._pat_parsable, "Fracture energy"
)
super().registerParamReal("beta",
aka._pat_readable | aka._pat_parsable)
self.registerInternalReal("delta_max", 1)
self.beta = 0.
self.sigma_c = 0.
self.delta_c = 0.
self.G_c = 0.
def initMaterial(self):
"""Initialize the parameters for the material."""
super().initMaterial()
self.sigma_c = self.getReal("sigma_c")
self.G_c = self.getReal("G_c")
self.beta = self.getReal("beta")
self.delta_c = 2 * self.G_c / self.sigma_c
def checkInsertion(self, check_only):
"""Check if need to insert a cohesive element."""
model = self.getModel()
facets = self.getFacetFilter()
inserter = model.getElementInserter()
for type_facet in facets.elementTypes(dim=(spatial_dimension - 1)):
facet_filter = facets(type_facet)
nb_facet = facet_filter.shape[0]
if nb_facet == 0:
continue
fe_engine = model.getFEEngine("FacetsFEEngine")
facets_check = inserter.getCheckFacets(type_facet)
insertion = inserter.getInsertionFacets(type_facet)
nb_quad_facet = fe_engine.getNbIntegrationPoints(type_facet)
normals = fe_engine.getNormalsOnIntegrationPoints(type_facet)
facets_stresses = model.getStressOnFacets(type_facet).reshape(
normals.shape[0] // nb_quad_facet,
nb_quad_facet,
2,
spatial_dimension,
spatial_dimension,
)
tangents = model.getTangents(type_facet)
for facet, facet_stresses in zip(facet_filter, facets_stresses):
if not facets_check[facet]:
continue
ref_stress = 0
for q, quad_stresses in enumerate(facet_stresses):
current_quad = facet * nb_quad_facet + q
normal = normals[current_quad, :].ravel()
tangent = tangents[current_quad, :].ravel()
stress_1 = quad_stresses.T[0]
stress_2 = quad_stresses.T[1]
avg_stress = stress_1 + stress_2 / 2.0
traction = avg_stress.dot(normal)
n = traction.dot(normal)
n = max(0, n)
t = traction.dot(tangent)
ref_stress = max(
ref_stress, np.sqrt(n * n + t * t / (self.beta**2))
)
if ref_stress > self.sigma_c:
insertion[facet] = True
# constitutive law
def computeTraction(self, normals, el_type, ghost_type):
"""Compute the traction for a given opening."""
openings = self.getOpening(el_type, ghost_type)
tractions = self.getTraction(el_type, ghost_type)
delta_max = self.getInternalReal("delta_max")(el_type)
for el in range(normals.shape[0]):
normal = normals[el].ravel()
opening = openings[el].ravel()
delta_n = opening.dot(normal) * normal
delta_s = opening - delta_n
delta = (
self.beta * np.linalg.norm(delta_s) ** 2 +
np.linalg.norm(delta_n) ** 2
)
delta_max[el] = max(delta_max[el], delta)
tractions[el, :] = (
(delta * delta_s + delta_n)
* self.sigma_c
/ delta
* (1 - delta / self.delta_c)
)
def allocator(_dim, unused, model, _id):
"""Register the material to the material factory."""
return LinearCohesive(model, _id)
mat_factory = aka.MaterialFactory.getInstance()
mat_factory.registerAllocator("local_cohesive", allocator)
# -------------------------------------------------------------------------
# Initialization
# -------------------------------------------------------------------------
aka.parseInput("local_material.dat")
mesh = aka.Mesh(spatial_dimension)
mesh.read("plate.msh")
model = aka.SolidMechanicsModelCohesive(mesh)
model.initFull(_analysis_method=aka._static, _is_extrinsic=True)
model.initNewSolver(aka._explicit_lumped_mass)
model.setBaseName("plate")
model.addDumpFieldVector("displacement")
model.addDumpFieldVector("external_force")
model.addDumpField("strain")
model.addDumpField("stress")
model.addDumpField("blocked_dofs")
model.setBaseNameToDumper("cohesive elements", "cohesive")
model.addDumpFieldVectorToDumper("cohesive elements", "displacement")
model.addDumpFieldToDumper("cohesive elements", "damage")
model.addDumpFieldVectorToDumper("cohesive elements", "tractions")
model.addDumpFieldVectorToDumper("cohesive elements", "opening")
# -------------------------------------------------------------------------
# Boundary conditions
# -------------------------------------------------------------------------
model.applyBC(aka.FixedValue(0.0, aka._x), "XBlocked")
model.applyBC(aka.FixedValue(0.0, aka._y), "YBlocked")
trac = np.zeros(spatial_dimension)
traction = 0.095
trac[int(aka._y)] = traction
model.getExternalForce()[:] = 0
model.applyBC(aka.FromTraction(trac), "Traction")
print("Solve for traction ", traction)
solver = model.getNonLinearSolver("static")
solver.set("max_iterations", 100)
solver.set("threshold", 1e-10)
solver.set("convergence_type", aka.SolveConvergenceCriteria.residual)
model.solveStep("static")
model.dump()
model.dump("cohesive elements")
model.setTimeStep(model.getStableTimeStep() * 0.1)
maxsteps = 100
for i in range(0, maxsteps):
print("{0}/{1}".format(i, maxsteps))
model.checkCohesiveStress()
model.solveStep("explicit_lumped")
if i % 10 == 0:
model.dump()
model.dump("cohesive elements")
-
diff --git a/examples/python/cohesive/plate.py b/examples/python/cohesive/plate.py
index a0cba5249..aa86acb45 100644
--- a/examples/python/cohesive/plate.py
+++ b/examples/python/cohesive/plate.py
@@ -1,106 +1,92 @@
#!/usr/bin/env python3
__copyright__ = "Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale" \
" de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \
" en Mécanique des Solides)"
__license__ = "LGPLv3"
-
-"""plate.py: Python example: plate with a hole breaking with cohesive
-elements"""
-
-__author__ = "Guillaume Anciaux"
-__credits__ = [
- "Guillaume Anciaux ",
-]
-__copyright__ = (
- "Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale"
- " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation"
- " en Mécanique des Solides)"
-)
-__license__ = "LGPLv3"
-
import akantu as aka
import numpy as np
+
def set_dumpers(model):
model.setBaseName("plate")
model.addDumpFieldVector("displacement")
model.addDumpFieldVector("external_force")
model.addDumpField("strain")
model.addDumpField("stress")
model.addDumpField("blocked_dofs")
model.setBaseNameToDumper("cohesive elements", "cohesive")
model.addDumpFieldVectorToDumper("cohesive elements", "displacement")
model.addDumpFieldToDumper("cohesive elements", "damage")
model.addDumpFieldVectorToDumper("cohesive elements", "tractions")
model.addDumpFieldVectorToDumper("cohesive elements", "opening")
def solve(material_file, mesh_file, traction):
aka.parseInput(material_file)
spatial_dimension = 2
# -------------------------------------------------------------------------
# Initialization
# -------------------------------------------------------------------------
mesh = aka.Mesh(spatial_dimension)
mesh.read(mesh_file)
model = aka.SolidMechanicsModelCohesive(mesh)
model.initFull(_analysis_method=aka._static, _is_extrinsic=True)
model.initNewSolver(aka._explicit_lumped_mass)
set_dumpers(model)
# -------------------------------------------------------------------------
# Boundary conditions
# -------------------------------------------------------------------------
model.applyBC(aka.FixedValue(0.0, aka._x), "XBlocked")
model.applyBC(aka.FixedValue(0.0, aka._y), "YBlocked")
trac = np.zeros(spatial_dimension)
trac[int(aka._y)] = traction
print("Solve for traction ", traction)
model.getExternalForce()[:] = 0
model.applyBC(aka.FromTraction(trac), "Traction")
solver = model.getNonLinearSolver("static")
solver.set("max_iterations", 100)
solver.set("threshold", 1e-10)
solver.set("convergence_type", aka.SolveConvergenceCriteria.residual)
model.solveStep("static")
model.dump()
model.dump("cohesive elements")
model.setTimeStep(model.getStableTimeStep() * 0.1)
maxsteps = 100
for i in range(0, maxsteps):
print("{0}/{1}".format(i, maxsteps))
model.checkCohesiveStress()
model.solveStep("explicit_lumped")
if i % 10 == 0:
model.dump()
model.dump("cohesive elements")
# -----------------------------------------------------------------------------
# main
# -----------------------------------------------------------------------------
def main():
mesh_file = "plate.msh"
material_file = "material.dat"
traction = 0.095
solve(material_file, mesh_file, traction)
# -----------------------------------------------------------------------------
if __name__ == "__main__":
main()
diff --git a/examples/python/custom-material/custom-material.py b/examples/python/custom-material/custom-material.py
index 4385443b2..2764f2afb 100644
--- a/examples/python/custom-material/custom-material.py
+++ b/examples/python/custom-material/custom-material.py
@@ -1,193 +1,180 @@
#!/usr/bin/env python3
__copyright__ = "Copyright (©) 2016-2023 EPFL (Ecole Polytechnique Fédérale" \
" de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \
" en Mécanique des Solides)"
__license__ = "LGPLv3"
-""" custom-material.py: Custom material example"""
-
-__author__ = "Guillaume Anciaux"
-__credits__ = [
- "Guillaume Anciaux ",
-]
-__copyright__ = (
- "Copyright (©) 2016-2021 EPFL (Ecole Polytechnique Fédérale"
- " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation"
- " en Mécanique des Solides)"
-)
-__license__ = "LGPLv3"
-
import numpy as np
import akantu as aka
# ------------------------------------------------------------------------------
class LocalElastic(aka.Material):
def __init__(self, model, _id):
super().__init__(model, _id)
super().registerParamReal(
"E", aka._pat_readable | aka._pat_parsable, "Youngs modulus"
)
super().registerParamReal(
"nu", aka._pat_readable | aka._pat_parsable, "Poisson ratio"
)
def initMaterial(self):
nu = self.getReal("nu")
E = self.getReal("E")
self.mu = E / (2 * (1 + nu))
self.lame_lambda = nu * E / ((1.0 + nu) * (1.0 - 2.0 * nu))
# Second Lame coefficient (shear modulus)
self.lame_mu = E / (2.0 * (1.0 + nu))
super().initMaterial()
# declares all the parameters that are needed
def getPushWaveSpeed(self, element):
rho = self.getReal("rho")
return np.sqrt((self.lame_lambda + 2 * self.lame_mu) / rho)
# compute small deformation tensor
@staticmethod
def computeEpsilon(grad_u):
return 0.5 * (grad_u + np.einsum("aij->aji", grad_u))
# constitutive law
def computeStress(self, el_type, ghost_type):
grad_u = self.getGradU(el_type, ghost_type)
sigma = self.getStress(el_type, ghost_type)
n_quads = grad_u.shape[0]
grad_u = grad_u.reshape((n_quads, 2, 2))
epsilon = self.computeEpsilon(grad_u)
sigma = sigma.reshape((n_quads, 2, 2))
trace = np.einsum("aii->a", grad_u)
sigma[:, :, :] = (
np.einsum("a,ij->aij", trace, self.lame_lambda * np.eye(2))
+ 2.0 * self.lame_mu * epsilon
)
# constitutive law tangent modulii
def computeTangentModuli(self, el_type, tangent_matrix, ghost_type):
n_quads = tangent_matrix.shape[0]
tangent = tangent_matrix.reshape(n_quads, 3, 3)
Miiii = self.lame_lambda + 2 * self.lame_mu
Miijj = self.lame_lambda
Mijij = self.lame_mu
tangent[:, 0, 0] = Miiii
tangent[:, 1, 1] = Miiii
tangent[:, 0, 1] = Miijj
tangent[:, 1, 0] = Miijj
tangent[:, 2, 2] = Mijij
# computes the energy density
def computePotentialEnergy(self, el_type):
sigma = self.getStress(el_type)
grad_u = self.getGradU(el_type)
nquads = sigma.shape[0]
stress = sigma.reshape(nquads, 2, 2)
grad_u = grad_u.reshape((nquads, 2, 2))
epsilon = self.computeEpsilon(grad_u)
energy_density = self.getPotentialEnergy(el_type)
energy_density[:, 0] = 0.5 * np.einsum("aij,aij->a", stress, epsilon)
# register material to the MaterialFactory
def allocator(_dim, unused, model, _id):
return LocalElastic(model, _id)
mat_factory = aka.MaterialFactory.getInstance()
mat_factory.registerAllocator("local_elastic", allocator)
# ------------------------------------------------------------------------------
# main
# ------------------------------------------------------------------------------
spatial_dimension = 2
aka.parseInput("material.dat")
mesh_file = "bar.msh"
max_steps = 250
time_step = 1e-3
# ------------------------------------------------------------------------------
# Initialization
# ------------------------------------------------------------------------------
mesh = aka.Mesh(spatial_dimension)
mesh.read(mesh_file)
# parse input file
aka.parseInput("material.dat")
model = aka.SolidMechanicsModel(mesh)
model.initFull(_analysis_method=aka._explicit_lumped_mass)
model.setBaseName("waves")
model.addDumpFieldVector("displacement")
model.addDumpFieldVector("acceleration")
model.addDumpFieldVector("velocity")
model.addDumpFieldVector("internal_force")
model.addDumpFieldVector("external_force")
model.addDumpField("strain")
model.addDumpField("stress")
model.addDumpField("blocked_dofs")
# ------------------------------------------------------------------------------
# boundary conditions
# ------------------------------------------------------------------------------
model.applyBC(aka.FixedValue(0, aka._x), "XBlocked")
model.applyBC(aka.FixedValue(0, aka._y), "YBlocked")
# ------------------------------------------------------------------------------
# initial conditions
# ------------------------------------------------------------------------------
displacement = model.getDisplacement()
nb_nodes = mesh.getNbNodes()
position = mesh.getNodes()
pulse_width = 1
A = 0.01
for i in range(0, nb_nodes):
# Sinus * Gaussian
x = position[i, 0] - 5.0
L = pulse_width
k = 0.1 * 2 * np.pi * 3 / L
displacement[i, 0] = A * np.sin(k * x) * np.exp(-(k * x) * (k * x) / (L * L))
# ------------------------------------------------------------------------------
# timestep value computation
# ------------------------------------------------------------------------------
time_factor = 0.8
stable_time_step = model.getStableTimeStep() * time_factor
print("Stable Time Step = {0}".format(stable_time_step))
print("Required Time Step = {0}".format(time_step))
time_step = stable_time_step * time_factor
model.setTimeStep(time_step)
# ------------------------------------------------------------------------------
# loop for evolution of motion dynamics
# ------------------------------------------------------------------------------
print("step,step * time_step,epot,ekin,epot + ekin")
for step in range(0, max_steps + 1):
model.solveStep()
if step % 10 == 0:
model.dump()
epot = model.getEnergy("potential")
ekin = model.getEnergy("kinetic")
# output energy calculation to screen
print(
"{0},{1},{2},{3},{4}".format(step, step * time_step, epot, ekin, (epot + ekin))
)
diff --git a/examples/python/phase-field/phasefield-static.py b/examples/python/phase-field/phasefield-static.py
index 5644c1f21..dec93de04 100644
--- a/examples/python/phase-field/phasefield-static.py
+++ b/examples/python/phase-field/phasefield-static.py
@@ -1,138 +1,124 @@
#!/usr/bin/env python
# coding: utf-8
__copyright__ = "Copyright (©) 2021-2023 EPFL (Ecole Polytechnique Fédérale" \
" de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \
" en Mécanique des Solides)"
__license__ = "LGPLv3"
-""" phasefield-static.py: Static phase field example"""
-
-__author__ = "Mohit Pundir"
-__credits__ = [
- "Mohit Pundir ",
-]
-__copyright__ = (
- "Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale"
- " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation"
- " en Mécanique des Solides)"
-)
-__license__ = "LGPLv3"
-
-
import numpy as np
import akantu as aka
aka.parseInput("material_static.dat")
dim = 2
mesh = aka.Mesh(dim)
mesh.read("plate_static.msh")
model = aka.CouplerSolidPhaseField(mesh)
solid = model.getSolidMechanicsModel()
phase = model.getPhaseFieldModel()
solid.initFull(_analysis_method=aka._static)
solver = solid.getNonLinearSolver("static")
solver.set("max_iterations", 100)
solver.set("threshold", 1e-8)
solver.set("convergence_type", aka.SolveConvergenceCriteria.solution)
solid.getNewSolver(
"linear_static", aka.TimeStepSolverType.static, aka.NonLinearSolverType.linear
)
solid.setIntegrationScheme(
"linear_static", "displacement", aka.IntegrationSchemeType.pseudo_time
)
phase.initFull(_analysis_method=aka._static)
phase.getNewSolver(
"nonlinear_static",
aka.TimeStepSolverType.static,
aka.NonLinearSolverType.newton_raphson,
)
phase.setIntegrationScheme(
"nonlinear_static", "damage", aka.IntegrationSchemeType.pseudo_time
)
solver = phase.getNonLinearSolver("nonlinear_static")
solver.set("max_iterations", 100)
solver.set("threshold", 1e-4)
solver.set("convergence_type", aka.SolveConvergenceCriteria.solution)
solid.applyBC(aka.FixedValue(0, aka._y), "bottom")
solid.applyBC(aka.FixedValue(0, aka._x), "left")
# Initialization for bulk vizualisation
solid.setBaseName("phasefield-static")
solid.addDumpFieldVector("displacement")
solid.addDumpFieldVector("external_force")
solid.addDumpField("strain")
solid.addDumpField("stress")
solid.addDumpField("damage")
solid.addDumpField("blocked_dofs")
nb_dofs = solid.getMesh().getNbNodes() * dim
increment = solid.getIncrement()
displacement = solid.getDisplacement()
displacement = displacement.reshape(nb_dofs)
blocked_dofs = solid.getBlockedDOFs()
blocked_dofs = blocked_dofs.reshape(nb_dofs)
damage = phase.getDamage()
tolerance = 1e-5
steps = 1000
increment = 5e-6
for n in range(steps):
print("Computing iteration " + str(n + 1) + "/" + str(steps))
solid.applyBC(aka.IncrementValue(increment, aka._y), "top")
mask = blocked_dofs == False # NOQA: E712
iiter = 0
error_disp = 1
error_dam = 1
displacement_prev = displacement[mask].copy()
damage_prev = damage.copy()
damage_prev = damage_prev
# solve using staggered scheme
while error_disp > tolerance or error_dam > tolerance:
model.solve("linear_static", "")
displacement_new = displacement[mask]
damage_new = damage
delta_disp = displacement_new - displacement_prev
delta_dam = damage_new - damage_prev
error_disp = np.linalg.norm(delta_disp)
error_dam = np.linalg.norm(delta_dam)
iiter += 1
displacement_prev = displacement_new.copy()
damage_prev = damage_new.copy()
print(error_dam, error_disp)
if iiter > 500:
raise Exception("Convergence not reached")
if n % 50 == 0:
solid.dump()
solid.dump()
diff --git a/examples/static/static.cc b/examples/static/static.cc
index df5a193bb..ddc9059a0 100644
--- a/examples/static/static.cc
+++ b/examples/static/static.cc
@@ -1,69 +1,69 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "non_linear_solver.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
#define bar_length 0.01
#define bar_height 0.01
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material.dat", argc, argv);
const Int spatial_dimension = 2;
Mesh mesh(spatial_dimension);
mesh.read("square.msh");
SolidMechanicsModel model(mesh);
/// model initialization
model.initFull(_analysis_method = _static);
model.setBaseName("static");
model.addDumpFieldVector("displacement");
model.addDumpField("external_force");
model.addDumpField("internal_force");
model.addDumpField("grad_u");
model.addDumpField("stress");
/// Dirichlet boundary conditions
model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "Fixed_x");
model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "Fixed_y");
model.applyBC(BC::Dirichlet::FixedValue(0.0001, _y), "Traction");
model.dump();
auto & solver = model.getNonLinearSolver();
solver.set("max_iterations", 2);
solver.set("threshold", 2e-4);
solver.set("convergence_type", SolveConvergenceCriteria::_solution);
model.solveStep();
model.dump();
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/structural_mechanics/bernoulli_beam_2_example.cc b/examples/structural_mechanics/bernoulli_beam_2_example.cc
index 27e92545a..0762f4ce7 100644
--- a/examples/structural_mechanics/bernoulli_beam_2_example.cc
+++ b/examples/structural_mechanics/bernoulli_beam_2_example.cc
@@ -1,131 +1,131 @@
/**
* Copyright (©) 2011-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "mesh_accessor.hh"
#include "structural_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
#define TYPE _bernoulli_beam_2
using namespace akantu;
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize(argc, argv);
// Defining the mesh
Mesh beams(2);
const auto q = 6000.;
const auto L = 10.;
const auto M = -3600.; // Momentum at 3
auto nb_nodes = 3;
auto nb_element = nb_nodes - 1;
MeshAccessor mesh_accessor(beams);
Array & nodes = mesh_accessor.getNodes();
nodes.resize(nb_nodes);
beams.addConnectivityType(_bernoulli_beam_2);
auto & connectivity = mesh_accessor.getConnectivity(_bernoulli_beam_2);
connectivity.resize(nb_element);
nodes.zero();
nodes(1, 0) = 10;
nodes(2, 0) = 18;
for (int i = 0; i < nb_element; ++i) {
connectivity(i, 0) = i;
connectivity(i, 1) = i + 1;
}
mesh_accessor.makeReady();
// Defining the materials
StructuralMechanicsModel model(beams);
StructuralMaterial mat1;
mat1.E = 3e10;
mat1.I = 0.0025;
mat1.A = 0.01;
model.addMaterial(mat1);
StructuralMaterial mat2;
mat2.E = 3e10;
mat2.I = 0.00128;
mat2.A = 0.01;
model.addMaterial(mat2);
// Defining the forces
model.initFull();
auto & forces = model.getExternalForce();
auto & displacement = model.getDisplacement();
auto & boundary = model.getBlockedDOFs();
const auto & N_M = model.getStress(_bernoulli_beam_2);
auto & element_material = model.getElementMaterial(_bernoulli_beam_2);
boundary.set(false);
forces.zero();
displacement.zero();
element_material(1) = 1;
forces(0, 1) = -q * L / 2.;
forces(0, 2) = -q * L * L / 12.;
forces(1, 1) = -q * L / 2.;
forces(1, 2) = q * L * L / 12.;
forces(2, 2) = M;
forces(2, 0) = mat2.E * mat2.A / 18;
// Defining the boundary conditions
boundary(0, 0) = true;
boundary(0, 1) = true;
boundary(0, 2) = true;
boundary(1, 1) = true;
boundary(2, 1) = true;
model.addDumpFieldVector("displacement");
model.addDumpField("rotation");
model.addDumpFieldVector("force");
model.addDumpField("momentum");
model.solveStep();
model.assembleResidual();
// Post-Processing
std::cout << " d1 = " << displacement(1, 2) << std::endl;
std::cout << " d2 = " << displacement(2, 2) << std::endl;
std::cout << " d3 = " << displacement(1, 0) << std::endl;
std::cout << " M1 = " << N_M(0, 1) << std::endl;
std::cout << " M2 = " << N_M(2 * (nb_nodes - 2), 1) << std::endl;
model.dump();
finalize();
}
diff --git a/extra_packages/extra-materials/src/material_FE2/material_FE2.cc b/extra_packages/extra-materials/src/material_FE2/material_FE2.cc
index cf785b64f..caffc1c7b 100644
--- a/extra_packages/extra-materials/src/material_FE2/material_FE2.cc
+++ b/extra_packages/extra-materials/src/material_FE2/material_FE2.cc
@@ -1,204 +1,204 @@
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "material_FE2.hh"
#include "communicator.hh"
#include "solid_mechanics_model_RVE.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
template
MaterialFE2::MaterialFE2(SolidMechanicsModel & model,
const ID & id)
: Parent(model, id), C("material_stiffness", *this) {
AKANTU_DEBUG_IN();
this->C.initialize(voigt_h::size * voigt_h::size);
this->initialize();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
MaterialFE2::~MaterialFE2() = default;
/* -------------------------------------------------------------------------- */
template void MaterialFE2::initialize() {
this->registerParam("element_type", el_type, _triangle_3,
_pat_parsable | _pat_modifiable,
"element type in RVE mesh");
this->registerParam("mesh_file", mesh_file, _pat_parsable | _pat_modifiable,
"the mesh file for the RVE");
this->registerParam("nb_gel_pockets", nb_gel_pockets,
_pat_parsable | _pat_modifiable,
"the number of gel pockets in each RVE");
}
/* -------------------------------------------------------------------------- */
template
void MaterialFE2::initMaterial() {
AKANTU_DEBUG_IN();
Parent::initMaterial();
/// create a Mesh and SolidMechanicsModel on each integration point of the
/// material
auto C_it = this->C(this->el_type).begin(voigt_h::size, voigt_h::size);
for (auto && data :
enumerate(make_view(C(this->el_type), voigt_h::size, voigt_h::size))) {
auto q = std::get<0>(data);
auto & C = std::get<1>(data);
meshes.emplace_back(std::make_unique(
spatial_dimension, "RVE_mesh_" + std::to_string(q), q + 1));
auto & mesh = *meshes.back();
mesh.read(mesh_file);
RVEs.emplace_back(std::make_unique(
mesh, true, this->nb_gel_pockets, _all_dimensions,
"SMM_RVE_" + std::to_string(q), q + 1));
auto & RVE = *RVEs.back();
RVE.initFull(_analysis_method = _static);
/// compute intial stiffness of the RVE
RVE.homogenizeStiffness(C);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialFE2::computeStress(ElementType el_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
// Compute thermal stresses first
Parent::computeStress(el_type, ghost_type);
Array::const_scalar_iterator sigma_th_it =
this->sigma_th(el_type, ghost_type).begin();
// Wikipedia convention:
// 2*eps_ij (i!=j) = voigt_eps_I
// http://en.wikipedia.org/wiki/Voigt_notation
Array::const_matrix_iterator C_it =
this->C(el_type, ghost_type).begin(voigt_h::size, voigt_h::size);
// create vectors to store stress and strain in Voigt notation
// for efficient computation of stress
Vector voigt_strain(voigt_h::size);
Vector voigt_stress(voigt_h::size);
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
const Matrix & C_mat = *C_it;
const Real & sigma_th = *sigma_th_it;
/// copy strains in Voigt notation
for (Int I = 0; I < voigt_h::size; ++I) {
/// copy stress in
Real voigt_factor = voigt_h::factors[I];
UInt i = voigt_h::vec[I][0];
UInt j = voigt_h::vec[I][1];
voigt_strain(I) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.;
}
// compute stresses in Voigt notation
voigt_stress.mul(C_mat, voigt_strain);
/// copy stresses back in full vectorised notation
for (Int I = 0; I < voigt_h::size; ++I) {
UInt i = voigt_h::vec[I][0];
UInt j = voigt_h::vec[I][1];
sigma(i, j) = sigma(j, i) = voigt_stress(I) + (i == j) * sigma_th;
}
++C_it;
++sigma_th_it;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialFE2::computeTangentModuli(
ElementType el_type, Array & tangent_matrix, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Array::const_matrix_iterator C_it =
this->C(el_type, ghost_type).begin(voigt_h::size, voigt_h::size);
MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix);
tangent.copy(*C_it);
++C_it;
MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template
void MaterialFE2::advanceASR(
const Matrix & prestrain) {
AKANTU_DEBUG_IN();
for (auto && data :
zip(RVEs,
make_view(this->gradu(this->el_type), spatial_dimension,
spatial_dimension),
make_view(this->eigengradu(this->el_type), spatial_dimension,
spatial_dimension),
make_view(this->C(this->el_type), voigt_h::size, voigt_h::size),
this->delta_T(this->el_type))) {
auto & RVE = *(std::get<0>(data));
/// apply boundary conditions based on the current macroscopic displ.
/// gradient
RVE.applyBoundaryConditions(std::get<1>(data));
/// apply homogeneous temperature field to each RVE to obtain thermoelastic
/// effect
RVE.applyHomogeneousTemperature(std::get<4>(data));
/// advance the ASR in every RVE
RVE.advanceASR(prestrain);
/// compute the average eigen_grad_u
RVE.homogenizeEigenGradU(std::get<2>(data));
/// compute the new effective stiffness of the RVE
RVE.homogenizeStiffness(std::get<3>(data));
}
AKANTU_DEBUG_OUT();
}
INSTANTIATE_MATERIAL(material_FE2, MaterialFE2);
} // namespace akantu
diff --git a/extra_packages/extra-materials/src/material_FE2/material_FE2.hh b/extra_packages/extra-materials/src/material_FE2/material_FE2.hh
index 8ab7018e0..976bcc871 100644
--- a/extra_packages/extra-materials/src/material_FE2/material_FE2.hh
+++ b/extra_packages/extra-materials/src/material_FE2/material_FE2.hh
@@ -1,115 +1,115 @@
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "material.hh"
#include "material_thermal.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_MATERIAL_FE_2_HH_
#define AKANTU_MATERIAL_FE_2_HH_
namespace akantu {
class SolidMechanicsModelRVE;
}
namespace akantu {
/* -------------------------------------------------------------------------- */
/// /!\ This material works ONLY for meshes with a single element type!!!!!
/* -------------------------------------------------------------------------- */
/**
* MaterialFE2
*
* parameters in the material files :
* - mesh_file
*/
template class MaterialFE2 : public MaterialThermal {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
private:
typedef MaterialThermal Parent;
public:
MaterialFE2(SolidMechanicsModel & model, const ID & id = "");
virtual ~MaterialFE2();
typedef VoigtHelper voigt_h;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
virtual void initMaterial();
/// constitutive law for all element of a type
virtual void computeStress(ElementType el_type,
GhostType ghost_type = _not_ghost);
/// compute the tangent stiffness matrix for an element type
void computeTangentModuli(ElementType el_type, Array & tangent_matrix,
GhostType ghost_type = _not_ghost);
/// advance alkali-silica reaction
void advanceASR(const Matrix & prestrain);
private:
void initialize();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// Underlying RVE at each integration point
std::vector> RVEs;
/// Meshes for all RVEs
std::vector> meshes;
/// the element type of the associated mesh (this material handles only one
/// type!!)
ElementType el_type;
/// the name of RVE mesh file
ID mesh_file;
/// Elastic stiffness tensor at each Gauss point (in voigt notation)
InternalField C;
/// number of gel pockets in each underlying RVE
UInt nb_gel_pockets;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "material_FE2_inline_impl.hh"
} // namespace akantu
#endif /* AKANTU_MATERIAL_FE_2_HH_ */
diff --git a/extra_packages/extra-materials/src/material_FE2/material_FE2_inline_impl.hh b/extra_packages/extra-materials/src/material_FE2/material_FE2_inline_impl.hh
index 9d9896118..0709c7cd5 100644
--- a/extra_packages/extra-materials/src/material_FE2/material_FE2_inline_impl.hh
+++ b/extra_packages/extra-materials/src/material_FE2/material_FE2_inline_impl.hh
@@ -1,20 +1,20 @@
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
diff --git a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc
index 04eef1a38..06144882a 100644
--- a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc
+++ b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc
@@ -1,596 +1,596 @@
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model_RVE.hh"
#include "element_group.hh"
#include "material_damage_iterative.hh"
#include "node_group.hh"
#include "non_linear_solver.hh"
#include "non_local_manager.hh"
#include "parser.hh"
#include "sparse_matrix.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh,
bool use_RVE_mat_selector,
UInt nb_gel_pockets, Int dim,
const ID & id)
: SolidMechanicsModel(mesh, dim, id), volume(0.),
use_RVE_mat_selector(use_RVE_mat_selector),
nb_gel_pockets(nb_gel_pockets), nb_dumps(0) {
AKANTU_DEBUG_IN();
/// find the four corner nodes of the RVE
findCornerNodes();
/// remove the corner nodes from the surface node groups:
/// This most be done because corner nodes a not periodic
mesh.getElementGroup("top").removeNode(corner_nodes(2));
mesh.getElementGroup("top").removeNode(corner_nodes(3));
mesh.getElementGroup("left").removeNode(corner_nodes(3));
mesh.getElementGroup("left").removeNode(corner_nodes(0));
mesh.getElementGroup("bottom").removeNode(corner_nodes(1));
mesh.getElementGroup("bottom").removeNode(corner_nodes(0));
mesh.getElementGroup("right").removeNode(corner_nodes(2));
mesh.getElementGroup("right").removeNode(corner_nodes(1));
const auto & bottom = mesh.getElementGroup("bottom").getNodeGroup();
bottom_nodes.insert(bottom.begin(), bottom.end());
const auto & left = mesh.getElementGroup("left").getNodeGroup();
left_nodes.insert(left.begin(), left.end());
// /// enforce periodicity on the displacement fluctuations
// auto surface_pair_1 = std::make_pair("top", "bottom");
// auto surface_pair_2 = std::make_pair("right", "left");
// SurfacePairList surface_pairs_list;
// surface_pairs_list.push_back(surface_pair_1);
// surface_pairs_list.push_back(surface_pair_2);
// TODO: To Nicolas correct the PBCs
// this->setPBC(surface_pairs_list);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModelRVE::~SolidMechanicsModelRVE() = default;
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::initFullImpl(const ModelOptions & options) {
AKANTU_DEBUG_IN();
auto options_cp(options);
options_cp.analysis_method = AnalysisMethod::_static;
SolidMechanicsModel::initFullImpl(options_cp);
// this->initMaterials();
auto & fem = this->getFEEngine("SolidMechanicsFEEngine");
/// compute the volume of the RVE
GhostType gt = _not_ghost;
for (auto element_type :
this->mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) {
Array Volume(this->mesh.getNbElement(element_type) *
fem.getNbIntegrationPoints(element_type),
1, 1.);
this->volume = fem.integrate(Volume, element_type);
}
std::cout << "The volume of the RVE is " << this->volume << std::endl;
/// dumping
std::stringstream base_name;
base_name << this->id;
this->setBaseName(base_name.str());
this->addDumpFieldVector("displacement");
this->addDumpField("stress");
this->addDumpField("grad_u");
this->addDumpField("eigen_grad_u");
this->addDumpField("blocked_dofs");
this->addDumpField("material_index");
this->addDumpField("damage");
this->addDumpField("Sc");
this->addDumpField("external_force");
this->addDumpField("equivalent_stress");
this->addDumpField("internal_force");
this->addDumpField("delta_T");
this->dump();
this->nb_dumps += 1;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::applyBoundaryConditions(
const Matrix & displacement_gradient) {
AKANTU_DEBUG_IN();
/// get the position of the nodes
const Array & pos = mesh.getNodes();
/// storage for the coordinates of a given node and the displacement that will
/// be applied
Vector x(spatial_dimension);
Vector appl_disp(spatial_dimension);
/// fix top right node
UInt node = this->corner_nodes(2);
x(0) = pos(node, 0);
x(1) = pos(node, 1);
appl_disp.mul(displacement_gradient, x);
(*this->blocked_dofs)(node, 0) = true;
(*this->displacement)(node, 0) = appl_disp(0);
(*this->blocked_dofs)(node, 1) = true;
(*this->displacement)(node, 1) = appl_disp(1);
// (*this->blocked_dofs)(node,0) = true; (*this->displacement)(node,0) = 0.;
// (*this->blocked_dofs)(node,1) = true; (*this->displacement)(node,1) = 0.;
/// apply Hx at all the other corner nodes; H: displ. gradient
node = this->corner_nodes(0);
x(0) = pos(node, 0);
x(1) = pos(node, 1);
appl_disp.mul(displacement_gradient, x);
(*this->blocked_dofs)(node, 0) = true;
(*this->displacement)(node, 0) = appl_disp(0);
(*this->blocked_dofs)(node, 1) = true;
(*this->displacement)(node, 1) = appl_disp(1);
node = this->corner_nodes(1);
x(0) = pos(node, 0);
x(1) = pos(node, 1);
appl_disp.mul(displacement_gradient, x);
(*this->blocked_dofs)(node, 0) = true;
(*this->displacement)(node, 0) = appl_disp(0);
(*this->blocked_dofs)(node, 1) = true;
(*this->displacement)(node, 1) = appl_disp(1);
node = this->corner_nodes(3);
x(0) = pos(node, 0);
x(1) = pos(node, 1);
appl_disp.mul(displacement_gradient, x);
(*this->blocked_dofs)(node, 0) = true;
(*this->displacement)(node, 0) = appl_disp(0);
(*this->blocked_dofs)(node, 1) = true;
(*this->displacement)(node, 1) = appl_disp(1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::applyHomogeneousTemperature(
const Real & temperature) {
for (UInt m = 0; m < this->getNbMaterials(); ++m) {
Material & mat = this->getMaterial(m);
const ElementTypeMapArray & filter_map = mat.getElementFilter();
// Loop over all element types
for (auto && type : filter_map.elementTypes(spatial_dimension)) {
const Array & filter = filter_map(type);
if (filter.size() == 0)
continue;
Array & delta_T = mat.getArray("delta_T", type);
Array::scalar_iterator delta_T_it = delta_T.begin();
Array::scalar_iterator delta_T_end = delta_T.end();
for (; delta_T_it != delta_T_end; ++delta_T_it) {
*delta_T_it = temperature;
}
}
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::findCornerNodes() {
AKANTU_DEBUG_IN();
// find corner nodes
const auto & position = mesh.getNodes();
const auto & lower_bounds = mesh.getLowerBounds();
const auto & upper_bounds = mesh.getUpperBounds();
AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!");
corner_nodes.resize(4);
corner_nodes.set(UInt(-1));
for (auto && data : enumerate(make_view(position, spatial_dimension))) {
auto node = std::get<0>(data);
const auto & X = std::get<1>(data);
auto distance = X.distance(lower_bounds);
// node 1
if (Math::are_float_equal(distance, 0)) {
corner_nodes(0) = node;
}
// node 2
else if (Math::are_float_equal(X(_x), upper_bounds(_x)) &&
Math::are_float_equal(X(_y), lower_bounds(_y))) {
corner_nodes(1) = node;
}
// node 3
else if (Math::are_float_equal(X(_x), upper_bounds(_x)) &&
Math::are_float_equal(X(_y), upper_bounds(_y))) {
corner_nodes(2) = node;
}
// node 4
else if (Math::are_float_equal(X(_x), lower_bounds(_x)) &&
Math::are_float_equal(X(_y), upper_bounds(_y))) {
corner_nodes(3) = node;
}
}
for (UInt i = 0; i < corner_nodes.size(); ++i) {
if (corner_nodes(i) == UInt(-1))
AKANTU_ERROR("The corner node " << i + 1 << " wasn't found");
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::advanceASR(const Matrix & prestrain) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!");
/// apply the new eigenstrain
for (auto element_type :
mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) {
Array & prestrain_vect =
const_cast &>(this->getMaterial("gel").getInternal(
"eigen_grad_u")(element_type));
auto prestrain_it =
prestrain_vect.begin(spatial_dimension, spatial_dimension);
auto prestrain_end =
prestrain_vect.end(spatial_dimension, spatial_dimension);
for (; prestrain_it != prestrain_end; ++prestrain_it)
(*prestrain_it) = prestrain;
}
/// advance the damage
MaterialDamageIterative<2> & mat_paste =
dynamic_cast &>(*this->materials[1]);
MaterialDamageIterative<2> & mat_aggregate =
dynamic_cast &>(*this->materials[0]);
UInt nb_damaged_elements = 0;
Real max_eq_stress_aggregate = 0;
Real max_eq_stress_paste = 0;
auto & solver = this->getNonLinearSolver();
solver.set("max_iterations", 2);
solver.set("threshold", 1e-6);
solver.set("convergence_type", SolveConvergenceCriteria::_solution);
do {
this->solveStep();
/// compute damage
max_eq_stress_aggregate = mat_aggregate.getNormMaxEquivalentStress();
max_eq_stress_paste = mat_paste.getNormMaxEquivalentStress();
nb_damaged_elements = 0;
if (max_eq_stress_aggregate > max_eq_stress_paste)
nb_damaged_elements = mat_aggregate.updateDamage();
else if (max_eq_stress_aggregate < max_eq_stress_paste)
nb_damaged_elements = mat_paste.updateDamage();
else
nb_damaged_elements =
(mat_paste.updateDamage() + mat_aggregate.updateDamage());
std::cout << "the number of damaged elements is " << nb_damaged_elements
<< std::endl;
} while (nb_damaged_elements);
if (this->nb_dumps % 10 == 0) {
this->dump();
}
this->nb_dumps += 1;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index,
const ID & field_type) {
AKANTU_DEBUG_IN();
auto & fem = this->getFEEngine("SolidMechanicsFEEngine");
Real average = 0;
GhostType gt = _not_ghost;
for (auto element_type :
mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) {
if (field_type == "stress") {
for (UInt m = 0; m < this->materials.size(); ++m) {
const auto & stress_vec = this->materials[m]->getStress(element_type);
const auto & elem_filter =
this->materials[m]->getElementFilter(element_type);
Array int_stress_vec(elem_filter.size(),
spatial_dimension * spatial_dimension,
"int_of_stress");
fem.integrate(stress_vec, int_stress_vec,
spatial_dimension * spatial_dimension, element_type,
_not_ghost, elem_filter);
for (UInt k = 0; k < elem_filter.size(); ++k)
average += int_stress_vec(k, row_index * spatial_dimension +
col_index); // 3 is the value
// for the yy (in
// 3D, the value is
// 4)
}
} else if (field_type == "strain") {
for (UInt m = 0; m < this->materials.size(); ++m) {
const auto & gradu_vec = this->materials[m]->getGradU(element_type);
const auto & elem_filter =
this->materials[m]->getElementFilter(element_type);
Array int_gradu_vec(elem_filter.size(),
spatial_dimension * spatial_dimension,
"int_of_gradu");
fem.integrate(gradu_vec, int_gradu_vec,
spatial_dimension * spatial_dimension, element_type,
_not_ghost, elem_filter);
for (UInt k = 0; k < elem_filter.size(); ++k)
/// averaging is done only for normal components, so stress and strain
/// are equal
average +=
0.5 *
(int_gradu_vec(k, row_index * spatial_dimension + col_index) +
int_gradu_vec(k, col_index * spatial_dimension + row_index));
}
} else if (field_type == "eigen_grad_u") {
for (UInt m = 0; m < this->materials.size(); ++m) {
const auto & eigen_gradu_vec =
this->materials[m]->getInternal("eigen_grad_u")(element_type);
const auto & elem_filter =
this->materials[m]->getElementFilter(element_type);
Array int_eigen_gradu_vec(elem_filter.size(),
spatial_dimension * spatial_dimension,
"int_of_gradu");
fem.integrate(eigen_gradu_vec, int_eigen_gradu_vec,
spatial_dimension * spatial_dimension, element_type,
_not_ghost, elem_filter);
for (UInt k = 0; k < elem_filter.size(); ++k)
/// averaging is done only for normal components, so stress and strain
/// are equal
average +=
int_eigen_gradu_vec(k, row_index * spatial_dimension + col_index);
}
} else {
AKANTU_ERROR("Averaging not implemented for this field!!!");
}
}
return average / this->volume;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::homogenizeStiffness(Matrix & C_macro) {
AKANTU_DEBUG_IN();
const Int dim = 2;
AKANTU_DEBUG_ASSERT(this->spatial_dimension == dim,
"Is only implemented for 2D!!!");
/// apply three independent loading states to determine C
/// 1. eps_el = (1;0;0) 2. eps_el = (0,1,0) 3. eps_el = (0,0,0.5)
/// clear the eigenstrain
Matrix zero_eigengradu(dim, dim, 0.);
GhostType gt = _not_ghost;
for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) {
auto & prestrain_vect =
const_cast &>(this->getMaterial("gel").getInternal(
"eigen_grad_u")(element_type));
auto prestrain_it =
prestrain_vect.begin(spatial_dimension, spatial_dimension);
auto prestrain_end =
prestrain_vect.end(spatial_dimension, spatial_dimension);
for (; prestrain_it != prestrain_end; ++prestrain_it)
(*prestrain_it) = zero_eigengradu;
}
/// storage for results of 3 different loading states
UInt voigt_size = VoigtHelper::size;
Matrix stresses(voigt_size, voigt_size, 0.);
Matrix strains(voigt_size, voigt_size, 0.);
Matrix H(dim, dim, 0.);
/// save the damage state before filling up cracks
// ElementTypeMapReal saved_damage("saved_damage");
// saved_damage.initialize(getFEEngine(), _nb_component = 1, _default_value =
// 0);
// this->fillCracks(saved_damage);
/// virtual test 1:
H(0, 0) = 0.01;
this->performVirtualTesting(H, stresses, strains, 0);
/// virtual test 2:
H.zero();
H(1, 1) = 0.01;
this->performVirtualTesting(H, stresses, strains, 1);
/// virtual test 3:
H.zero();
H(0, 1) = 0.01;
this->performVirtualTesting(H, stresses, strains, 2);
/// drain cracks
// this->drainCracks(saved_damage);
/// compute effective stiffness
Matrix eps_inverse(voigt_size, voigt_size);
eps_inverse.inverse(strains);
C_macro.mul(stresses, eps_inverse);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::performVirtualTesting(const Matrix & H,
Matrix & eff_stresses,
Matrix & eff_strains,
const UInt test_no) {
AKANTU_DEBUG_IN();
this->applyBoundaryConditions(H);
auto & solver = this->getNonLinearSolver();
solver.set("max_iterations", 2);
solver.set("threshold", 1e-6);
solver.set("convergence_type", SolveConvergenceCriteria::_solution);
this->solveStep();
/// get average stress and strain
eff_stresses(0, test_no) = this->averageTensorField(0, 0, "stress");
eff_strains(0, test_no) = this->averageTensorField(0, 0, "strain");
eff_stresses(1, test_no) = this->averageTensorField(1, 1, "stress");
eff_strains(1, test_no) = this->averageTensorField(1, 1, "strain");
eff_stresses(2, test_no) = this->averageTensorField(1, 0, "stress");
eff_strains(2, test_no) = 2. * this->averageTensorField(1, 0, "strain");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::homogenizeEigenGradU(
Matrix & eigen_gradu_macro) {
AKANTU_DEBUG_IN();
eigen_gradu_macro(0, 0) = this->averageTensorField(0, 0, "eigen_grad_u");
eigen_gradu_macro(1, 1) = this->averageTensorField(1, 1, "eigen_grad_u");
eigen_gradu_macro(0, 1) = this->averageTensorField(0, 1, "eigen_grad_u");
eigen_gradu_macro(1, 0) = this->averageTensorField(1, 0, "eigen_grad_u");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::initMaterials() {
AKANTU_DEBUG_IN();
// make sure the material are instantiated
if (!are_materials_instantiated)
instantiateMaterials();
if (use_RVE_mat_selector) {
const Vector & lowerBounds = mesh.getLowerBounds();
const Vector & upperBounds = mesh.getUpperBounds();
Real bottom = lowerBounds(1);
Real top = upperBounds(1);
Real box_size = std::abs(top - bottom);
Real eps = box_size * 1e-6;
auto tmp = std::make_shared(*this, box_size, "gel",
this->nb_gel_pockets, eps);
tmp->setFallback(material_selector);
material_selector = tmp;
}
this->assignMaterialToElements();
// synchronize the element material arrays
this->synchronize(SynchronizationTag::_material_id);
for (auto & material : materials) {
/// init internals properties
const auto type = material->getID();
if (type.find("material_FE2") != std::string::npos)
continue;
material->initMaterial();
}
this->synchronize(SynchronizationTag::_smm_init_mat);
if (this->non_local_manager) {
this->non_local_manager->initialize();
}
// SolidMechanicsModel::initMaterials();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::fillCracks(ElementTypeMapReal & saved_damage) {
const auto & mat_gel = this->getMaterial("gel");
Real E_gel = mat_gel.get("E");
Real E_homogenized = 0.;
for (auto && mat : materials) {
if (mat->getName() == "gel" || mat->getName() == "FE2_mat")
continue;
Real E = mat->get("E");
auto & damage = mat->getInternal("damage");
for (auto && type : damage.elementTypes()) {
const auto & elem_filter = mat->getElementFilter(type);
auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type);
auto sav_dam_it =
make_view(saved_damage(type), nb_integration_point).begin();
for (auto && data :
zip(elem_filter, make_view(damage(type), nb_integration_point))) {
auto el = std::get<0>(data);
auto & dam = std::get<1>(data);
Vector sav_dam = sav_dam_it[el];
sav_dam = dam;
for (auto q : arange(dam.size())) {
E_homogenized = (E_gel - E) * dam(q) + E;
dam(q) = 1. - (E_homogenized / E);
}
}
}
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::drainCracks(
const ElementTypeMapReal & saved_damage) {
for (auto && mat : materials) {
if (mat->getName() == "gel" || mat->getName() == "FE2_mat")
continue;
auto & damage = mat->getInternal("damage");
for (auto && type : damage.elementTypes()) {
const auto & elem_filter = mat->getElementFilter(type);
auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type);
auto sav_dam_it =
make_view(saved_damage(type), nb_integration_point).begin();
for (auto && data :
zip(elem_filter, make_view(damage(type), nb_integration_point))) {
auto el = std::get<0>(data);
auto & dam = std::get<1>(data);
Vector sav_dam = sav_dam_it[el];
dam = sav_dam;
}
}
}
}
} // namespace akantu
diff --git a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh
index 5c2fcc243..217f66618 100644
--- a/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh
+++ b/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.hh
@@ -1,224 +1,224 @@
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
- *
+ *
* Akantu 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.
- *
+ *
* Akantu 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 Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_SOLID_MECHANICS_MODEL_RVE_HH_
#define AKANTU_SOLID_MECHANICS_MODEL_RVE_HH_
/* -------------------------------------------------------------------------- */
#include "aka_grid_dynamic.hh"
#include "solid_mechanics_model.hh"
#include
/* -------------------------------------------------------------------------- */
namespace akantu {
class SolidMechanicsModelRVE : public SolidMechanicsModel {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector = true,
UInt nb_gel_pockets = 400,
Int spatial_dimension = _all_dimensions,
const ID & id = "solid_mechanics_model");
virtual ~SolidMechanicsModelRVE();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
void initFullImpl(const ModelOptions & option) override;
/// initialize the materials
void initMaterials() override;
public:
/// apply boundary contions based on macroscopic deformation gradient
virtual void
applyBoundaryConditions(const Matrix & displacement_gradient);
/// apply homogeneous temperature field from the macroscale level to the RVEs
virtual void applyHomogeneousTemperature(const Real & temperature);
/// advance the reactions -> grow gel and apply homogenized properties
void advanceASR(const Matrix & prestrain);
/// compute average stress or strain in the model
Real averageTensorField(UInt row_index, UInt col_index,
const ID & field_type);
/// compute effective stiffness of the RVE
void homogenizeStiffness(Matrix