From 9b10fde3328d7520721dcf86e76a14c5fcfa93ab Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Wed, 12 Sep 2018 08:34:03 +0200
Subject: [PATCH] [PL] Impl. of volumetric source terms.
First approach use the original bulk mesh as integration domain.
---
.../CreateVolumetricSourceTerm.cpp | 50 +++++++++
.../SourceTerms/CreateVolumetricSourceTerm.h | 43 ++++++++
ProcessLib/SourceTerms/SourceTerm.h | 4 +-
.../SourceTerms/VolumetricSourceTerm.cpp | 61 +++++++++++
ProcessLib/SourceTerms/VolumetricSourceTerm.h | 38 +++++++
.../SourceTerms/VolumetricSourceTermFEM.h | 101 ++++++++++++++++++
6 files changed, 296 insertions(+), 1 deletion(-)
create mode 100644 ProcessLib/SourceTerms/CreateVolumetricSourceTerm.cpp
create mode 100644 ProcessLib/SourceTerms/CreateVolumetricSourceTerm.h
create mode 100644 ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
create mode 100644 ProcessLib/SourceTerms/VolumetricSourceTerm.h
create mode 100644 ProcessLib/SourceTerms/VolumetricSourceTermFEM.h
diff --git a/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.cpp b/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.cpp
new file mode 100644
index 00000000000..c7b424b0d5f
--- /dev/null
+++ b/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.cpp
@@ -0,0 +1,50 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateVolumetricSourceTerm.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/FileTools.h"
+#include "MeshLib/Mesh.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+#include "VolumetricSourceTerm.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<SourceTerm> createVolumetricSourceTerm(
+ BaseLib::ConfigTree const& config,
+ MeshLib::Mesh const& source_term_mesh,
+ NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order, unsigned const shapefunction_order,
+ int const variable_id, int const component_id)
+{
+ //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__type}
+ config.checkConfigParameter("type", "Volumetric");
+
+ DBUG("Constructing VolumetricSourceTerm from config.");
+
+ // source term field name
+ auto const& volumetric_source_term_parameter_name =
+ //! \ogs_file_param{prj__process_variables__process_variable_source_terms__volumetric_source_term__parameter}
+ config.getConfigParameter<std::string>("parameter");
+ auto& volumetric_source_term = findParameter<double>(
+ volumetric_source_term_parameter_name, parameters, 1);
+
+ DBUG("Using '%s` as volumetric source term parameter.",
+ volumetric_source_term.name.c_str());
+
+ return std::make_unique<VolumetricSourceTerm>(
+ source_term_mesh, source_term_dof_table, integration_order,
+ shapefunction_order, variable_id, component_id,
+ volumetric_source_term);
+}
+
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.h b/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.h
new file mode 100644
index 00000000000..2f1fcb62c54
--- /dev/null
+++ b/ProcessLib/SourceTerms/CreateVolumetricSourceTerm.h
@@ -0,0 +1,43 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+}
+
+namespace ProcessLib
+{
+struct ParameterBase;
+class SourceTerm;
+
+std::unique_ptr<SourceTerm> createVolumetricSourceTerm(
+ BaseLib::ConfigTree const& config,
+ MeshLib::Mesh const& source_term_mesh,
+ NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order, unsigned const shapefunction_order,
+ int const variable_id, int const component_id);
+
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/SourceTerm.h b/ProcessLib/SourceTerms/SourceTerm.h
index 89d9e595f13..bbbd878a3e9 100644
--- a/ProcessLib/SourceTerms/SourceTerm.h
+++ b/ProcessLib/SourceTerms/SourceTerm.h
@@ -17,7 +17,8 @@ namespace ProcessLib
class SourceTerm
{
public:
- SourceTerm(const NumLib::LocalToGlobalIndexMap& source_term_dof_table)
+ explicit SourceTerm(
+ const NumLib::LocalToGlobalIndexMap& source_term_dof_table)
: _source_term_dof_table(source_term_dof_table)
{
}
@@ -25,6 +26,7 @@ public:
virtual void integrate(const double t, GlobalVector& b) const = 0;
virtual ~SourceTerm() = default;
+
protected:
NumLib::LocalToGlobalIndexMap const& _source_term_dof_table;
};
diff --git a/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp b/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
new file mode 100644
index 00000000000..166405fa84c
--- /dev/null
+++ b/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
@@ -0,0 +1,61 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "VolumetricSourceTerm.h"
+
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+
+namespace ProcessLib
+{
+VolumetricSourceTerm::VolumetricSourceTerm(
+ MeshLib::Mesh const& source_term_mesh,
+ NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+ unsigned const integration_order, unsigned const shapefunction_order,
+ int const variable_id, int const component_id,
+ Parameter<double> const& volumetric_source_term)
+ : SourceTerm(source_term_dof_table),
+ _volumetric_source_term(volumetric_source_term)
+{
+ // check basic data consistency
+ if (variable_id >=
+ static_cast<int>(source_term_dof_table.getNumberOfVariables()))
+ {
+ OGS_FATAL(
+ "Variable id too high. Actual value: %d, maximum value: %d.",
+ variable_id,
+ source_term_dof_table.getNumberOfVariables());
+ }
+ if (component_id >=
+ source_term_dof_table.getNumberOfVariableComponents(variable_id))
+ {
+ OGS_FATAL(
+ "Component id too high. Actual value: %d, maximum value: %d.",
+ component_id,
+ source_term_dof_table.getNumberOfVariableComponents(variable_id));
+ }
+
+ ProcessLib::createLocalAssemblers<VolumetricSourceTermLocalAssembler>(
+ source_term_mesh.getDimension(), source_term_mesh.getElements(),
+ source_term_dof_table, shapefunction_order, _local_assemblers,
+ source_term_mesh.isAxiallySymmetric(), integration_order,
+ _volumetric_source_term);
+}
+
+void VolumetricSourceTerm::integrate(const double t,
+ GlobalVector& b) const
+{
+ DBUG("Assemble VolumetricSourceTerm.");
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberOnDereferenced(
+ &VolumetricSourceTermLocalAssemblerInterface::integrate,
+ _local_assemblers, _source_term_dof_table, t, b);
+}
+
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/VolumetricSourceTerm.h b/ProcessLib/SourceTerms/VolumetricSourceTerm.h
new file mode 100644
index 00000000000..9f26f93d227
--- /dev/null
+++ b/ProcessLib/SourceTerms/VolumetricSourceTerm.h
@@ -0,0 +1,38 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+#include <vector>
+
+#include "SourceTerm.h"
+#include "VolumetricSourceTermFEM.h"
+
+namespace ProcessLib
+{
+class VolumetricSourceTerm final : public SourceTerm
+{
+public:
+ VolumetricSourceTerm(
+ MeshLib::Mesh const& source_term_mesh,
+ NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+ unsigned const integration_order, unsigned const shapefunction_order,
+ int const variable_id, int const component_id,
+ Parameter<double> const& volumetric_source_term);
+
+ void integrate(const double t, GlobalVector& b) const;
+
+private:
+ Parameter<double> const& _volumetric_source_term;
+ std::vector<std::unique_ptr<VolumetricSourceTermLocalAssemblerInterface>>
+ _local_assemblers;
+};
+
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/VolumetricSourceTermFEM.h b/ProcessLib/SourceTerms/VolumetricSourceTermFEM.h
new file mode 100644
index 00000000000..bb165e4f652
--- /dev/null
+++ b/ProcessLib/SourceTerms/VolumetricSourceTermFEM.h
@@ -0,0 +1,101 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <vector>
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+namespace ProcessLib
+{
+class VolumetricSourceTermLocalAssemblerInterface
+{
+public:
+ virtual void integrate(
+ std::size_t const id,
+ NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+ double const t, GlobalVector& b) = 0;
+ virtual ~VolumetricSourceTermLocalAssemblerInterface() = default;
+};
+
+const unsigned NUM_NODAL_DOF = 1;
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class VolumetricSourceTermLocalAssembler final
+ : public VolumetricSourceTermLocalAssemblerInterface
+{
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+ using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits<
+ ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
+
+ using NodalMatrixType = typename LocalAssemblerTraits::LocalMatrix;
+ using NodalVectorType = typename LocalAssemblerTraits::LocalVector;
+ using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+
+public:
+ VolumetricSourceTermLocalAssembler(
+ MeshLib::Element const& element,
+ std::size_t const local_matrix_size,
+ bool is_axially_symmetric,
+ unsigned const integration_order,
+ Parameter<double> const& volumetric_source_term)
+ : _volumetric_source_term(volumetric_source_term),
+ _integration_method(integration_order),
+ _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ element, is_axially_symmetric, _integration_method)),
+ _local_rhs(local_matrix_size)
+ {
+ }
+
+ void integrate(std::size_t const id,
+ NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+ double const t, GlobalVector& b) override
+ {
+ _local_rhs.setZero();
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition pos;
+ pos.setElementID(id);
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ pos.setIntegrationPoint(ip);
+ auto const& sm = _shape_matrices[ip];
+ auto const& wp = _integration_method.getWeightedPoint(ip);
+ auto const st_val = _volumetric_source_term(t, pos)[0];
+
+ _local_rhs.noalias() +=
+ sm.N * st_val * sm.detJ * sm.integralMeasure * wp.getWeight();
+ }
+ auto const indices = NumLib::getIndices(id, source_term_dof_table);
+ b.add(indices, _local_rhs);
+ }
+
+private:
+ Parameter<double> const& _volumetric_source_term;
+
+ IntegrationMethod const _integration_method;
+ std::vector<ShapeMatrices, Eigen::aligned_allocator<ShapeMatrices>>
+ _shape_matrices;
+ NodalVectorType _local_rhs;
+};
+
+} // namespace ProcessLib
--
GitLab