initial commit

.vscode/settings.json

+{
+    "files.associations": {
+        "array": "cpp",
+        "bit": "cpp",
+        "*.tcc": "cpp",
+        "cctype": "cpp",
+        "clocale": "cpp",
+        "cmath": "cpp",
+        "complex": "cpp",
+        "cstddef": "cpp",
+        "cstdint": "cpp",
+        "cstdio": "cpp",
+        "cstdlib": "cpp",
+        "cstring": "cpp",
+        "cwchar": "cpp",
+        "cwctype": "cpp",
+        "deque": "cpp",
+        "set": "cpp",
+        "unordered_map": "cpp",
+        "unordered_set": "cpp",
+        "vector": "cpp",
+        "exception": "cpp",
+        "fstream": "cpp",
+        "functional": "cpp",
+        "initializer_list": "cpp",
+        "iosfwd": "cpp",
+        "iostream": "cpp",
+        "istream": "cpp",
+        "limits": "cpp",
+        "memory": "cpp",
+        "new": "cpp",
+        "numeric": "cpp",
+        "optional": "cpp",
+        "ostream": "cpp",
+        "sstream": "cpp",
+        "stdexcept": "cpp",
+        "streambuf": "cpp",
+        "string": "cpp",
+        "string_view": "cpp",
+        "system_error": "cpp",
+        "type_traits": "cpp",
+        "tuple": "cpp",
+        "typeinfo": "cpp",
+        "utility": "cpp",
+        "cstdarg": "cpp",
+        "ctime": "cpp",
+        "atomic": "cpp",
+        "strstream": "cpp",
+        "bitset": "cpp",
+        "chrono": "cpp",
+        "forward_list": "cpp",
+        "list": "cpp",
+        "map": "cpp",
+        "algorithm": "cpp",
+        "iterator": "cpp",
+        "memory_resource": "cpp",
+        "random": "cpp",
+        "ratio": "cpp",
+        "iomanip": "cpp",
+        "thread": "cpp",
+        "cfenv": "cpp",
+        "cinttypes": "cpp",
+        "typeindex": "cpp",
+        "variant": "cpp"
+    },
+    "editor.formatOnSave": true
+}
\ No newline at end of file

CMakeLists.txt

+cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
+project(ML4WaveDispersion LANGUAGES C CXX)
+
+add_subdirectory(ae108)
+
+add_subdirectory(validation)
+
+add_executable(DispersionSolver DispersionSolver.cc)
+target_link_libraries(DispersionSolver
+                      PRIVATE ae108::solve
+                      PRIVATE ae108::assembly
+                      PRIVATE ae108::elements
+                      PRIVATE ae108::cppslepc
+                      PRIVATE ae108::cmdline
+) 
\ No newline at end of file

DispersionSolver.cc

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this fKpath except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+/////////////////////////////////////////////////////////////////////
+// compare to https://doi.org/10.1038/s41598-019-46089-9
+/////////////////////////////////////////////////////////////////////
+
+#include "ae108/assembly//plugins/AssembleMassMatrixPlugin.h"
+#include "ae108/assembly/Assembler.h"
+#include "ae108/cmdline/CommandLineOptionParser.h"
+#include "ae108/cpppetsc/Context.h"
+#include "ae108/cpppetsc/GeneralizedMeshBoundaryCondition.h"
+#include "ae108/cpppetsc/Mesh.h"
+#include "ae108/cpppetsc/ParallelComputePolicy.h"
+#include "ae108/cpppetsc/Vector.h"
+#include "ae108/cpppetsc/asTransformedMatrix.h"
+#include "ae108/cpppetsc/asTransposedMatrix.h"
+#include "ae108/cpppetsc/computeElementsOfMatrix.h"
+#include "ae108/cppslepc/Context.h"
+#include "ae108/elements/TimoshenkoBeamElementWithMass.h"
+#include "ae108/elements/mesh/read_mesh_from_file.h"
+#include "ae108/elements/mesh/refine_segment_mesh.h"
+#include "ae108/elements/tensor/as_vector.h"
+#include "ae108/solve/boundaryConditionsToTransform.h"
+#include "utilities/assign_lattice_to.h"
+#include "utilities/circular_beam_properties.h"
+#include "utilities/compute_eigenvalues.h"
+#include "utilities/read_config.h"
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+using namespace ae108;
+
+using Policy = cpppetsc::ParallelComputePolicy;
+using Context = cppslepc::Context<Policy>;
+using Mesh = cpppetsc::Mesh<Policy>;
+using Vector = cpppetsc::Vector<Policy>;
+using BoundaryCondition = cpppetsc::GeneralizedMeshBoundaryCondition<Mesh>;
+
+constexpr auto topological_dimension = Mesh::TopologicalDimension{1};
+constexpr auto coordinate_dimension = Mesh::CoordinateDimension{3};
+
+using Point = std::array<Mesh::real_type, coordinate_dimension>;
+
+using Element = elements::TimoshenkoBeamElementWithMass<
+    coordinate_dimension, Vector::value_type, Vector::real_type>;
+
+using Properties =
+    elements::TimoshenkoBeamWithMassProperties<Mesh::real_type,
+                                               coordinate_dimension>;
+
+using Plugins =
+    assembly::FeaturePlugins<assembly::plugins::AssembleEnergyPlugin,
+                             assembly::plugins::AssembleForceVectorPlugin,
+                             assembly::plugins::AssembleStiffnessMatrixPlugin,
+                             assembly::plugins::AssembleMassMatrixPlugin>;
+
+using Assembler = assembly::Assembler<Element, Plugins, Policy>;
+
+using elements::tensor::as_vector;
+
+int main(int argc, char **argv) {
+  const auto context = Context(&argc, &argv);
+
+  auto index = Mesh::size_type(0);
+  std::string config_file = "/mnt/io/validation/data/cube.config";
+  std::string input_folder = "/mnt/io/validation/data/";
+  std::string output_folder = "/mnt/io/";
+  auto number_of_eigenfrequencies = Mesh::size_type(18);
+  auto max_element_length = Mesh::real_type(0.004);
+  cmdline::CommandLineOptionParser(std::cerr)
+      .withOption("config_file,c", "Config file.", &config_file)
+      .withOption("index,n", "Index of configuration.", &index)
+      .withOption("input_folder,i", "Input folder.", &input_folder)
+      .withOption("output_folder,o", "Output folder.", &output_folder)
+      .withOption("number_of_eigenfrequencies,e",
+                  "Requested number of eigenfrequencies.",
+                  &number_of_eigenfrequencies)
+      .withOption("max_element_length,m", "Maximum length of a beam element.",
+                  &max_element_length)
+      .parse(argc, argv);
+
+  std::ifstream fConfig(config_file);
+  const auto config = read_config<coordinate_dimension>(fConfig, index);
+
+  std::ifstream fKpath(input_folder + config.name + ".kpath");
+  FILE *out = fopen((output_folder + config.name + ".dispersion").c_str(), "w");
+
+  auto properties = circular_beam_properties<Properties>(
+      config.young_modulus, config.poisson_ratio, config.radius,
+      config.density);
+
+  const auto &lattice_vectors = config.lattice_vectors;
+
+  const auto geometry = elements::mesh::refine_segment_mesh(
+      elements::mesh::read_mesh_from_file<Point>(input_folder + config.name +
+                                                 ".dat"),
+      max_element_length);
+
+  const auto mesh = Mesh::fromConnectivity(
+      topological_dimension, coordinate_dimension, geometry.connectivity(),
+      geometry.number_of_positions(), Element::degrees_of_freedom(), 0);
+
+  auto assembler = Assembler();
+
+  auto element_source = ae108::homogenization::utilities::assign_lattice_to(
+      geometry.connectivity(), geometry.positions(), lattice_vectors);
+
+  for (const auto &element : mesh.localElements()) {
+    elements::tensor::Tensor<Mesh::real_type, coordinate_dimension>
+        element_axis;
+    as_vector(&element_axis) =
+        as_vector(&geometry.position_of_vertex(
+            geometry.connectivity().at(element.index()).at(1))) -
+        as_vector(&geometry.position_of_vertex(
+            geometry.connectivity().at(element.index()).at(0)));
+
+    const auto weight =
+        1. / std::count(element_source.begin(), element_source.end(),
+                        element_source[element.index()]);
+
+    assembler.emplaceElement(
+        element,
+        Element::Element(weight * timoshenko_beam_stiffness_matrix(element_axis,
+                                                                   properties)),
+        weight *
+            timoshenko_beam_consistent_mass_matrix(element_axis, properties));
+  }
+
+  const auto K = [&]() {
+    auto K = Mesh::matrix_type::fromMesh(mesh);
+    assembler.assembleStiffnessMatrix(Mesh::vector_type::fromLocalMesh(mesh),
+                                      Element::Time{0.}, &K);
+    K.finalize();
+    return K;
+  }();
+
+  const auto M = [&]() {
+    auto M = Mesh::matrix_type::fromMesh(mesh);
+    assembler.assembleMassMatrix(&M);
+    M.finalize();
+    return M;
+  }();
+
+  const auto source_of = ae108::homogenization::utilities::assign_lattice_to(
+      geometry.positions(), lattice_vectors);
+
+  std::string line;
+  while (std::getline(fKpath, line)) {
+    std::stringstream linestream(line);
+    std::string value;
+    Point wave_vector({0, 0, 0});
+    std::size_t i = 0;
+    for (std::size_t i = 0; i < coordinate_dimension; i++) {
+      std::getline(linestream, value, ',');
+      wave_vector[i] = std::atof(value.c_str());
+    }
+    const auto T = [&](const Point &wave_vector) {
+      std::vector<BoundaryCondition> boundary_conditions;
+      for (const auto &vertex : mesh.localVertices()) {
+        const Mesh::size_type &target = vertex.index();
+        const Mesh::size_type &source = source_of[target];
+        if (source != target) {
+          const auto k = [&] {
+            Point translation;
+            as_vector(&translation) =
+                as_vector(&geometry.position_of_vertex(target)) -
+                as_vector(&geometry.position_of_vertex(source));
+            return PETSC_i *
+                   as_vector(&wave_vector).dot(as_vector(&translation));
+          }();
+          for (Mesh::size_type dof = 0; dof < vertex.numberOfDofs(); dof++) {
+            boundary_conditions.push_back(
+                {{target, dof}, {{std::exp(k), {source, dof}}}, 0.});
+          }
+        }
+      }
+      return solve::boundaryConditionsToTransform(boundary_conditions, mesh)
+          .matrix;
+    }(wave_vector);
+
+    auto Tt = cpppetsc::asTransposedMatrix(&T);
+    const auto KT = cpppetsc::asTransformedMatrix(&K, &Tt);
+    const auto MT = cpppetsc::asTransformedMatrix(&M, &Tt);
+
+    const auto result = cppslepc::computeGeneralizedEigenvaluesWithSINVERT(
+        cpppetsc::computeElementsOfMatrix(KT),
+        cpppetsc::computeElementsOfMatrix(MT), number_of_eigenfrequencies);
+
+    PetscFPrintf(Policy::communicator(), out, "%f,%f,%f", wave_vector[0],
+                 wave_vector[1], wave_vector[2]);
+    for (std::size_t n = 0; n < number_of_eigenfrequencies; n++)
+      PetscFPrintf(Policy::communicator(), out, ",%f",
+                   std::sqrt(result[n]).real() / (2 * M_PI));
+    PetscFPrintf(Policy::communicator(), out, "\n");
+  }
+  fKpath.close();
+  fclose(out);
+}
\ No newline at end of file

docker-compose.yml

@@ -15,11 +15,11 @@
 version: "3.3"
 services:
   dev:
+    volumes:
+    - .:/mnt/io
     build:
       args:
         PETSC_SCALAR_TYPE: complex
       context: ./docker
       dockerfile: Dockerfile
-    volumes:
-    - .:/mnt/io
     tty: true

docker/Dockerfile

@@ -40,6 +40,8 @@ RUN pip3 install \
         black==20.8b1 \
         pylint==2.6.0 \
         sympy==1.7.1 \
+        matplotlib \
+        pymatgen \
         gcovr
 
 RUN apt-get update && \

docker/Dockerfile_docu

+# © 2020 ETH Zurich, Mechanics and Materials Lab
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+ARG BASE_IMAGE=debian:buster
+ARG DEBIAN_FRONTEND=noninteractive
+
+FROM $BASE_IMAGE
+
+RUN apt-get update && apt-get install -y \
+      doxygen \
+      git \
+      pandoc \
+      python3 python3-pip \
+      texlive texlive-latex-extra latexmk
+
+RUN pip3 install sphinx breathe

tools/findKPath.py

+from pymatgen.core.structure import Structure
+from pymatgen.symmetry.bandstructure import HighSymmKpath
+
+import numpy as np
+
+lattice = np.loadtxt('lattice.txt')
+
+# scaling factor to satisfy pymatgen
+fac =1
+
+s = Structure(np.nan_to_num(lattice) / fac, [1], [[0, 0, 0]])
+print(s)
+h = HighSymmKpath(s, symprec=0.001, angle_tolerance=1, atol=1e-06)
+
+kpath = h.get_kpoints(line_density=1)
+kpath = ([k / fac for k in kpath[0]], kpath[1])  # back scale
+
+f = open('kpath.txt', 'w')
+for label, k in zip(kpath[1], kpath[0]):
+
+    f.write(str(k[0]) + "," + str(k[1]) + "," + str(k[2]))
+
+    if label:
+        f.write("," + label)
+
+    f.write("\n")
+f.close()

tools/plotting.py

+import numpy as np
+from matplotlib import pyplot as plt
+
+def plot_dispersion(kpath, eigfreq):
+
+    from matplotlib import pyplot as plt
+    import numpy as np
+
+    plt.rc('text', usetex=True)
+    plt.rc('font', family='serif')
+
+    # klabel = np.array(kpath[1])
+    # idx = np.where(klabel <> u'')[0]
+    #
+    il = range(0, len(kpath))
+    # for i in range(len(idx) - 1):
+    #
+    #     if idx[i + 1] == idx[i] + 1:
+    #
+    #         il[idx[i] + 1:] = [l - 1 for l in il[idx[i] + 1:]]
+    #
+    #         if klabel[idx[i]] <> klabel[idx[i + 1]]:
+    #             klabel[idx[i]] = klabel[idx[i]] + klabel[idx[i + 1]]
+    #             klabel[idx[i + 1]] = klabel[idx[i]]
+
+    plt.plot(il, eigfreq)
+
+    #for i in [il[j] for j in idx]:
+     #   plt.axvline(x=i, color='black')
+
+
+    plt.xlim(0, np.max(il))
+
+    #plt.xticks([il[j] for j in idx], klabel[idx])
+    plt.show()
+
+
+name = "cube"
+path = "/mnt/io/"
+#path = "/home/bastian/Workspace/cppPNCOpt/build/Applications/"
+fkpath = path+"data/"+name+'.kpath'
+fdpath = path+name+".dispersion"
+
+plt.rc('text', usetex=True)
+plt.rc('font', family='serif')
+
+dispersion = np.loadtxt(fdpath,delimiter=',')
+
+with open(fkpath, "r") as f:
+    kpath = [x.split(",") for x in f.readlines()]
+
+klabel=[]
+il=[]
+for i, k in enumerate(kpath):
+    if len(k) > 3:
+        klabel.append(k[3].rstrip())
+        il.append(i)
+
+
+
+np.set_printoptions(precision=8, edgeitems=50,  linewidth=2)
+eigfreq = dispersion[:,3:]
+kpath = dispersion[:,:3]
+
+#plt.rc('text', usetex=True)
+
+from matplotlib.patches import Rectangle
+currentAxis = plt.gca()
+i_klabel = []
+for i in range(0,len(il),2):
+
+    idb = il[i]
+    idu = il[i+1]
+
+    plt.axvline(x=idb - i, color='black', linestyle='--',linewidth=0.8)
+    plt.plot(range(idb - i, idu - i), eigfreq[range(idb, idu)], color='black',marker=",",linestyle="None")
+    # plt.gca().set_prop_cycle(None)
+
+    i_klabel.append(idb - i)
+
+    fu = np.nanmin(eigfreq[idb:idu], axis=0)[1:]
+    fl = np.nanmax(eigfreq[idb:idu], axis=0)[:-1]
+
+    pb= fu-fl
+
+    # print pb
+    #
+    for j,b in enumerate(pb):
+         if b>2e-2:
+
+            w = idu - idb -1
+            h = fu[j] - fl[j]
+
+            currentAxis.add_patch(Rectangle((idb-i, fl[j]), w, h, alpha=1,color="#98FB98"))
+
+
+i_klabel.append(idu - i-1)
+fu = np.nanmin(eigfreq, axis=0)[1:]
+fl = np.nanmax(eigfreq, axis=0)[:-1]
+full_bandgap = fu - fl
+for j, b in enumerate(full_bandgap):
+    if b > 2e-2:
+
+        w = i_klabel[-1]
+        h = fu[j] - fl[j]
+
+        currentAxis.add_patch(Rectangle((0, fl[j]), w, h, alpha=1, color="#6B8E23"))
+
+# # k-path labelling
+klabel2 =[]
+double = False
+for i in range(len(klabel)-1):
+
+    idb = i
+    idu = i+1
+    if il[idb] == il[idu]-1:
+        if klabel[idb] != klabel[idu]:
+            klabel2.append("$"+klabel[idb]+r"\vert "+klabel[idb+1]+"$")
+            double = True
+    else:
+        if not double:
+            klabel2.append("$"+klabel[idb]+"$")
+        else:
+            double = False
+klabel2.append("$"+klabel[idu]+"$")
+
+plt.xticks(i_klabel, klabel2)
+plt.xlim(0,len(kpath)-len(klabel))
+plt.ylim(0,7000)
+# plt.title('dispersion curve')
+plt.xlabel('k-space position')
+plt.ylabel('frequency $[Hz]$')
+plt.savefig("test3.svg")
\ No newline at end of file

utilities/Config.h

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#pragma once
+
+template <std::size_t Dimension> struct Config {
+  double radius;
+  double young_modulus;
+  double poisson_ratio;
+  double density;
+  std::array<std::array<double, Dimension>, Dimension> lattice_vectors;
+  std::string name;
+};
\ No newline at end of file

utilities/assign_lattice_to.h

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#pragma once
+
+#include "find_source_of_target.h"
+#include "lattice_vector_multiplier.h"
+#include <vector>
+
+namespace ae108 {
+namespace homogenization {
+namespace utilities {
+
+template <class value_type, std::size_t dimension_>
+std::vector<std::size_t> assign_lattice_to(
+    const std::vector<std::vector<std::size_t>> &connectivities,
+    const std::vector<std::array<value_type, dimension_>> &positions,
+    const std::array<std::array<value_type, dimension_>, dimension_>
+        &lattice_vectors) noexcept {
+
+  const auto multiplier = lattice_vector_multiplier<double, dimension_, 1>();
+
+  std::vector<std::size_t> source(connectivities.size());
+  for (std::size_t target = 0; target < connectivities.size(); target++)
+    source[target] = find_source_of_target(positions, connectivities,
+                                           lattice_vectors, multiplier, target);
+  return source;
+}
+
+template <class value_type, std::size_t dimension_>
+std::vector<std::size_t> assign_lattice_to(
+    const std::vector<std::array<value_type, dimension_>> &positions,
+    const std::array<std::array<value_type, dimension_>, dimension_>
+        &lattice_vectors) noexcept {
+
+  auto connectivities = std::vector<std::vector<std::size_t>>(positions.size());
+  std::size_t i = 0;
+  std::for_each(connectivities.begin(), connectivities.end(),
+                [&i](std::vector<std::size_t> &n) { n.push_back(i++); });
+
+  return assign_lattice_to(connectivities, positions, lattice_vectors);
+}
+
+} // namespace utilities
+} // namespace homogenization
+} // namespace ae108

utilities/cartesian_product.h

+// © 2021 ETH Zurich, Mechanics and Materials Lab
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+