From 55b331e622d3d30a190a77fc3c61532153fc9d69 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Mon, 15 Jan 2018 02:09:12 +0100
Subject: [PATCH] [PL] SD: Set integration point initial conditions.
---
ProcessLib/Output/IntegrationPointWriter.cpp | 48 +++++++++++++++++
ProcessLib/Output/IntegrationPointWriter.h | 6 +++
.../LocalAssemblerInterface.h | 4 ++
.../SmallDeformation/SmallDeformationFEM.h | 46 ++++++++++++++++
.../SmallDeformationProcess-impl.h | 53 +++++++++++++++++++
5 files changed, 157 insertions(+)
diff --git a/ProcessLib/Output/IntegrationPointWriter.cpp b/ProcessLib/Output/IntegrationPointWriter.cpp
index fc4249e864e..5fca5b46382 100644
--- a/ProcessLib/Output/IntegrationPointWriter.cpp
+++ b/ProcessLib/Output/IntegrationPointWriter.cpp
@@ -71,6 +71,22 @@ static void addIntegrationPointMetaData(
std::back_inserter(dictionary));
}
+/// For the given json object and the name extract integration point meta data,
+/// or fail if no meta data was found for the given name.
+static ProcessLib::IntegrationPointMetaData extractIntegrationPointMetaData(
+ json const& meta_data, std::string const& name)
+{
+ for (auto const& md : meta_data["integration_point_arrays"])
+ {
+ if (md["name"] == name)
+ {
+ return {name, md["number_of_components"], md["integration_order"]};
+ }
+ }
+ OGS_FATAL("No integration point meta data with name \"%s\" found.",
+ name.c_str());
+}
+
namespace ProcessLib
{
void addIntegrationPointWriter(
@@ -88,4 +104,36 @@ void addIntegrationPointWriter(
addIntegrationPointMetaData(mesh, meta_data);
}
}
+
+IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Mesh const& mesh,
+ std::string const& name)
+{
+ if (!mesh.getProperties().existsPropertyVector<char>(
+ "IntegrationPointMetaData"))
+ {
+ OGS_FATAL(
+ "Integration point data '%s' is present in the vtk field "
+ "data but the required \"IntegrationPointMetaData\" array "
+ "is not available.",
+ name.c_str());
+ }
+ auto const& mesh_property_ip_meta_data =
+ *mesh.getProperties().template getPropertyVector<char>(
+ "IntegrationPointMetaData");
+
+ if (mesh_property_ip_meta_data.getMeshItemType() !=
+ MeshLib::MeshItemType::IntegrationPoint)
+ {
+ OGS_FATAL("IntegrationPointMetaData array must be field data.");
+ }
+
+ // Find the current integration point data entry and extract the
+ // meta data.
+ auto const ip_meta_data = extractIntegrationPointMetaData(
+ json::parse(mesh_property_ip_meta_data.begin(),
+ mesh_property_ip_meta_data.end()),
+ name);
+
+ return ip_meta_data;
+}
} // namespace ProcessLib
diff --git a/ProcessLib/Output/IntegrationPointWriter.h b/ProcessLib/Output/IntegrationPointWriter.h
index da772e81186..fc7d7274e65 100644
--- a/ProcessLib/Output/IntegrationPointWriter.h
+++ b/ProcessLib/Output/IntegrationPointWriter.h
@@ -9,6 +9,7 @@
*
*/
#include <memory>
+#include <nlohmann/json_fwd.hpp>
#include <vector>
#pragma once
@@ -51,4 +52,9 @@ struct IntegrationPointMetaData
int const integration_order;
};
+/// Returns integration point meta data for the given field name.
+///
+/// The data is read from a JSON encoded string stored in field data array.
+IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Mesh const& mesh,
+ std::string const& name);
} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
index 7166973c445..49c28f481c5 100644
--- a/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
+++ b/ProcessLib/SmallDeformation/LocalAssemblerInterface.h
@@ -26,6 +26,10 @@ struct SmallDeformationLocalAssemblerInterface
public MaterialForcesInterface,
public NumLib::ExtrapolatableElement
{
+ virtual std::size_t setIPDataInitialConditions(
+ std::string const& name, double const* values,
+ int const integration_order) = 0;
+
virtual std::vector<double> const& getIntPtFreeEnergyDensity(
const double /*t*/,
GlobalVector const& /*current_solution*/,
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index a93fe56df6c..8fee02c9f8c 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -151,6 +151,29 @@ public:
}
}
+ /// Returns number of read integration points.
+ std::size_t setIPDataInitialConditions(std::string const& name,
+ double const* values,
+ int const integration_order) override
+ {
+ if (integration_order !=
+ static_cast<int>(_integration_method.getIntegrationOrder()))
+ {
+ OGS_FATAL(
+ "Setting integration point initial conditions; The integration "
+ "order of the local assembler for element %d is different from "
+ "the integration order in the initial condition.",
+ _element.getID());
+ }
+
+ if (name == "sigma_ip")
+ {
+ return setSigma(values);
+ }
+
+ return 0;
+ }
+
void assemble(double const /*t*/, std::vector<double> const& /*local_x*/,
std::vector<double>& /*local_M_data*/,
std::vector<double>& /*local_K_data*/,
@@ -314,6 +337,29 @@ public:
return cache;
}
+ std::size_t setSigma(double const* values)
+ {
+ auto const kelvin_vector_size =
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value;
+ auto const n_integration_points = _ip_data.size();
+
+ std::vector<double> ip_sigma_values;
+ auto sigma_values =
+ Eigen::Map<Eigen::Matrix<double, kelvin_vector_size, Eigen::Dynamic,
+ Eigen::ColMajor> const>(
+ values, kelvin_vector_size, n_integration_points);
+
+ for (unsigned ip = 0; ip < n_integration_points; ++ip)
+ {
+ _ip_data[ip].sigma =
+ MathLib::KelvinVector::symmetricTensorToKelvinVector(
+ sigma_values.col(ip));
+ }
+
+ return n_integration_points;
+ }
+
// TODO (naumov) This method is same as getIntPtSigma but for arguments and
// the ordering of the cache_mat.
// There should be only one.
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h b/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
index 107bd51ec67..28900e122be 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess-impl.h
@@ -10,6 +10,7 @@
#pragma once
#include <cassert>
+#include <nlohmann/json.hpp>
#include "BaseLib/Functional.h"
#include "ProcessLib/Process.h"
@@ -72,6 +73,8 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
MeshLib::Mesh const& mesh,
unsigned const integration_order)
{
+ using nlohmann::json;
+
ProcessLib::SmallDeformation::createLocalAssemblers<
DisplacementDim, SmallDeformationLocalAssembler>(
mesh.getElements(), dof_table, _local_assemblers,
@@ -157,6 +160,56 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
makeExtrapolator(num_components, getExtrapolator(),
_local_assemblers, std::move(getIntPtValues)));
}
+
+ // Set initial conditions for integration point data.
+ for (auto const& ip_writer : _integration_point_writer)
+ {
+ // Find the mesh property with integration point writer's name.
+ auto const& name = ip_writer->name();
+ if (!mesh.getProperties().existsPropertyVector<double>(name))
+ {
+ continue;
+ }
+ auto const& mesh_property =
+ *mesh.getProperties().template getPropertyVector<double>(name);
+
+ // The mesh property must be defined on integration points.
+ if (mesh_property.getMeshItemType() !=
+ MeshLib::MeshItemType::IntegrationPoint)
+ {
+ continue;
+ }
+
+ auto const ip_meta_data = getIntegrationPointMetaData(mesh, name);
+
+ // Check the number of components.
+ if (ip_meta_data.n_components != mesh_property.getNumberOfComponents())
+ {
+ OGS_FATAL(
+ "Different number of components in meta data (%d) than in "
+ "the integration point field data for \"%s\": %d.",
+ ip_meta_data.n_components, name.c_str(),
+ mesh_property.getNumberOfComponents());
+ }
+
+ // Now we have a properly named vtk's field data array and the
+ // corresponding meta data.
+ std::size_t position = 0;
+ for (auto& local_asm : _local_assemblers)
+ {
+ std::size_t const integration_points_read =
+ local_asm->setIPDataInitialConditions(
+ name, &mesh_property[position],
+ ip_meta_data.integration_order);
+ if (integration_points_read == 0)
+ {
+ OGS_FATAL(
+ "No integration points read in the integration point "
+ "initial conditions set function.");
+ }
+ position += integration_points_read * ip_meta_data.n_components;
+ }
+ }
}
template <int DisplacementDim>
--
GitLab