To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 5ec890e0 authored by Bastian Telgen's avatar Bastian Telgen
Browse files

add output

parent c5026908
......@@ -24,9 +24,13 @@
#include "ae108/cpppetsc/Mesh.h"
#include "ae108/cpppetsc/ParallelComputePolicy.h"
#include "ae108/cpppetsc/Vector.h"
#include "ae108/cpppetsc/Viewer.h"
#include "ae108/cpppetsc/asTransformedMatrix.h"
#include "ae108/cpppetsc/asTransposedMatrix.h"
#include "ae108/cpppetsc/computeElementsOfMatrix.h"
#include "ae108/cpppetsc/createVectorFromSource.h"
#include "ae108/cpppetsc/setName.h"
#include "ae108/cpppetsc/writeToViewer.h"
#include "ae108/cppslepc/Context.h"
#include "ae108/elements/TimoshenkoBeamElementWithMass.h"
#include "ae108/elements/mesh/read_mesh_from_file.h"
......@@ -38,8 +42,8 @@
#include "utilities/compute_eigenvalues.h"
#include "utilities/read_config.h"
#include <Eigen/Dense>
#include <fstream>
#include <iostream>
#include <sstream>
using namespace ae108;
......@@ -49,6 +53,7 @@ using Context = cppslepc::Context<Policy>;
using Mesh = cpppetsc::Mesh<Policy>;
using Vector = cpppetsc::Vector<Policy>;
using BoundaryCondition = cpppetsc::GeneralizedMeshBoundaryCondition<Mesh>;
using Viewer = cpppetsc::Viewer<Policy>;
constexpr auto topological_dimension = Mesh::TopologicalDimension{1};
constexpr auto coordinate_dimension = Mesh::CoordinateDimension{3};
......@@ -75,14 +80,12 @@ using elements::tensor::as_vector;
int main(int argc, char **argv) {
const auto context = Context(&argc, &argv);
PetscPrintf(Policy::communicator(), "Read config\n");
auto index = Mesh::size_type(1);
std::string config_file = "/mnt/io/validation/data/cube.config";
std::string input_folder = "/mnt/io/validation/data/";
std::string output_folder = "/mnt/io/";
std::string config_file = "/mnt/io/tools/temp/lattice_vectors.csv";
std::string input_folder = "/mnt/io/tools/temp/";
std::string output_folder = "/mnt/io/tools/temp/";
auto number_of_eigenfrequencies = Mesh::size_type(18);
auto max_element_length = Mesh::real_type(0.004);
auto max_element_length = Mesh::real_type(10);
cmdline::CommandLineOptionParser(std::cerr)
.withOption("config_file,c", "Config file.", &config_file)
.withOption("index,n", "Index of configuration.", &index)
......@@ -95,6 +98,7 @@ int main(int argc, char **argv) {
&max_element_length)
.parse(argc, argv);
PetscPrintf(Policy::communicator(), "Reading config\n");
std::ifstream fConfig(config_file);
const auto config = read_config<coordinate_dimension>(fConfig, index);
......@@ -112,25 +116,44 @@ int main(int argc, char **argv) {
".dat"),
max_element_length);
PetscPrintf(Policy::communicator(), "Name: %s\n", config.name.c_str());
PetscPrintf(Policy::communicator(), " Name: %s\n", config.name.c_str());
PetscPrintf(Policy::communicator(),
"The mesh consists of %d elements and %d nodes\n",
" Number of elements: %d \n Number of nodes: %d\n",
geometry.connectivity().size(), geometry.number_of_positions());
PetscPrintf(Policy::communicator(), " Beam radius: %f \n", config.radius);
const auto mesh = Mesh::fromConnectivity(
topological_dimension, coordinate_dimension, geometry.connectivity(),
geometry.number_of_positions(), Element::degrees_of_freedom(), 0);
using DataSource = std::function<void(Mesh::size_type, Mesh::value_type *)>;
auto coordinates = ae108::cpppetsc::createVectorFromSource(
mesh, coordinate_dimension,
DataSource([&](const Mesh::size_type index, Mesh::value_type *const out) {
const auto &position = geometry.position_of_vertex(index);
std::copy(position.begin(), position.end(), out);
}));
ae108::cpppetsc::setName("coordinates", &coordinates);
auto viewer = Viewer::fromHdf5FilePath((config.name + ".ae108").c_str(),
Viewer::Mode::write);
ae108::cpppetsc::writeToViewer(mesh, coordinates, &viewer);
auto assembler = Assembler();
PetscPrintf(Policy::communicator(), "Find periodic couples\n");
PetscPrintf(Policy::communicator(), "Finding periodic couples\n");
auto element_source = ae108::homogenization::utilities::assign_lattice_to(
geometry.connectivity(), geometry.positions(), lattice_vectors);
const auto source_of = ae108::homogenization::utilities::assign_lattice_to(
geometry.positions(), lattice_vectors);
PetscPrintf(Policy::communicator(), "Assemble mass and stiffness matrix\n");
std::set<std::size_t> unique_sources(source_of.begin(), source_of.end());
PetscPrintf(Policy::communicator(),
" Number of sources: %d \n Number of eliminated targets: %d\n",
unique_sources.size(),
geometry.number_of_positions() - unique_sources.size());
double solid_volume = 0;
PetscPrintf(Policy::communicator(), "Assembling mass and stiffness matrix\n");
for (const auto &element : mesh.localElements()) {
elements::tensor::Tensor<Mesh::real_type, coordinate_dimension>
element_axis;
......@@ -143,15 +166,22 @@ int main(int argc, char **argv) {
const auto weight =
1. / std::count(element_source.begin(), element_source.end(),
element_source[element.index()]);
solid_volume = solid_volume +
weight * properties.area * as_vector(&element_axis).norm();
assembler.emplaceElement(
element,
Element::Element(weight * timoshenko_beam_stiffness_matrix(element_axis,
properties)),
weight *
timoshenko_beam_consistent_mass_matrix(element_axis, properties));
weight * timoshenko_beam_lumped_mass_matrix(element_axis, properties));
}
const auto lat = ae108::elements::tensor::as_matrix_of_rows(&lattice_vectors);
double domain = lat.row(0).dot(lat.row(1).cross(lat.row(2)));
PetscPrintf(Policy::communicator(), " Domain volume: %f \n", domain);
PetscPrintf(Policy::communicator(), " Volume fraction: %f \n",
solid_volume / domain);
const auto K = [&]() {
auto K = Mesh::matrix_type::fromMesh(mesh);
assembler.assembleStiffnessMatrix(Mesh::vector_type::fromLocalMesh(mesh),
......@@ -160,7 +190,6 @@ int main(int argc, char **argv) {
return K;
}();
PetscPrintf(Policy::communicator(), "Solve Bloch Eigenvalue Problems\n");
const auto M = [&]() {
auto M = Mesh::matrix_type::fromMesh(mesh);
assembler.assembleMassMatrix(&M);
......@@ -168,6 +197,7 @@ int main(int argc, char **argv) {
return M;
}();
PetscPrintf(Policy::communicator(), "Solving Bloch Eigenvalue Problems\n");
std::string line;
while (std::getline(fKpath, line)) {
std::stringstream linestream(line);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment