LehrFEM++ 1.0.0
A simple Finite Element Library for teaching
Loading...
Searching...
No Matches
fe_tools.h
1#ifndef LF_FETOOLS_H
2#define LF_FETOOLS_H
3/***************************************************************************
4 * LehrFEM++ - A simple C++ finite element libray for teaching
5 * Developed from 2018 at the Seminar of Applied Mathematics of ETH Zurich,
6 * lead developers Dr. R. Casagrande and Prof. R. Hiptmair
7 ***************************************************************************/
8
18#include <lf/assemble/assemble.h>
19#include <lf/mesh/utils/utils.h>
20#include <lf/quad/quad.h>
21
22#include "scalar_fe_space.h"
23
24namespace lf::fe {
25
26static const unsigned int kout_l2_qr = 1;
27static const unsigned int kout_l2_rsfvals = 2;
28
29namespace internal {
30
31// TODO(raffael) convert this method into a lambda function of
32// IntegrateMeshFunction() once
33// https://developercommunity.visualstudio.com/content/problem/200017/if-constexpr-in-lambda.html
34// is resolved.
35template <class MF, class QR_SELECTOR>
36auto LocalIntegral(const mesh::Entity &e, const QR_SELECTOR &qr_selector,
39 auto qr = qr_selector(e);
40 auto values = mf(e, qr.Points());
41 auto weights_ie =
42 (qr.Weights().cwiseProduct(e.Geometry()->IntegrationElement(qr.Points())))
43 .eval();
44 LF_ASSERT_MSG(values.size() == qr.NumPoints(),
45 "mf returns vector with wrong size.");
46 if constexpr (std::is_arithmetic_v<MfType>) { // NOLINT
47 auto value_m = Eigen::Map<Eigen::Matrix<MfType, 1, Eigen::Dynamic>>(
48 &values[0], 1, values.size());
49 return (value_m * weights_ie)(0);
50 }
51
52 if constexpr (base::EigenMatrix<MfType>) { // NOLINT
53 constexpr int size = MfType::SizeAtCompileTime;
54 if constexpr (size != Eigen::Dynamic) {
55 auto value_m = Eigen::Map<
56 Eigen::Matrix<typename MfType::Scalar, size, Eigen::Dynamic>>(
57 &values[0](0, 0), size, values.size());
58 MfType result;
59 auto result_m =
60 Eigen::Map<Eigen::Matrix<typename MfType::Scalar, size, 1>>(
61 &result(0, 0));
62 result_m = value_m * weights_ie;
63 return result;
64 }
65 }
66 // fallback: we cannot make any optimizations:
67 MfType temp = weights_ie(0) * values[0];
68 for (Eigen::Index i = 1; i < qr.NumPoints(); ++i) {
69 temp = temp + weights_ie(i) * values[i];
70 }
71 return temp;
72}
73}; // namespace internal
74
107template <mesh::utils::MeshFunction MF, class QR_SELECTOR,
108 class ENTITY_PREDICATE = base::PredicateTrue>
109 requires std::is_invocable_v<QR_SELECTOR, const mesh::Entity &>
111 const lf::mesh::Mesh &mesh, const MF &mf, const QR_SELECTOR &qr_selector,
112 const ENTITY_PREDICATE &ep = base::PredicateTrue{},
115
116 auto entities = mesh.Entities(codim);
117 auto result = internal::LocalIntegral(**entities.begin(), qr_selector, mf);
118 for (auto i = entities.begin() + 1; i != entities.end(); ++i) {
119 if (!ep(**i)) {
120 continue;
121 }
122 result = result + internal::LocalIntegral(**i, qr_selector, mf);
123 }
124 return result;
125}
126
153template <mesh::utils::MeshFunction MF,
154 class ENTITY_PREDICATE = base::PredicateTrue>
156 const lf::mesh::Mesh &mesh, const MF &mf, int quad_degree,
157 const ENTITY_PREDICATE &ep = base::PredicateTrue{},
159 std::array<quad::QuadRule, 5> qrs;
160 for (auto ref_el :
162 qrs[ref_el.Id()] = quad::make_QuadRule(ref_el, quad_degree);
163 }
165 mesh, mf, [&](const mesh::Entity &e) { return qrs[e.RefEl().Id()]; }, ep,
166 codim);
167}
168
169// ******************************************************************************
170
196template <typename SCALAR, mesh::utils::MeshFunction MF,
197 typename SELECTOR = base::PredicateTrue>
198auto NodalProjection(const lf::fe::ScalarFESpace<SCALAR> &fe_space, MF const &u,
199 SELECTOR &&pred = base::PredicateTrue{})
200 -> Eigen::Vector<decltype(SCALAR{0} * mesh::utils::MeshFunctionReturnType<
201 std::remove_reference_t<MF>>{0}),
202 Eigen::Dynamic> {
203 // choose scalar type so it can hold the scalar type of u as well as
204 // SCALAR
205 using scalarMF_t =
206 mesh::utils::MeshFunctionReturnType<std::remove_reference_t<MF>>;
207 using scalar_t = decltype(SCALAR(0) * scalarMF_t(0)); // NOLINT
208 // Return type, type for FE coefficient vector
209 using dof_vec_t = Eigen::Matrix<scalar_t, Eigen::Dynamic, 1>;
210
211 // Underlying mesh instance
212 const lf::mesh::Mesh &mesh{*fe_space.Mesh()};
213 // Fetch local-to-global index mapping for shape functions
214 const lf::assemble::DofHandler &dofh{fe_space.LocGlobMap()};
215
216 // Vector for returning expansion coefficients
217 dof_vec_t glob_dofvec(dofh.NumDofs());
218 glob_dofvec.setZero();
219
220 // Loop over all cells
221 for (const lf::mesh::Entity *cell : mesh.Entities(0)) {
222 if (!pred(*cell)) {
223 continue;
224 }
225 // Topological type of the cell
226 const lf::base::RefEl ref_el{cell->RefEl()};
227
228 // Information about local shape functions on reference element
229 auto ref_shape_fns = fe_space.ShapeFunctionLayout(*cell);
230 LF_ASSERT_MSG(ref_shape_fns, "reference shape function for "
231 << ref_el << " not available.");
232 // Obtain reference coordinates for evaluation nodes
233 const Eigen::MatrixXd ref_nodes(ref_shape_fns->EvaluationNodes());
234
235 // Collect values of function to be projected in a row vector
236 auto uvalvec = u(*cell, ref_nodes);
237
238 // Compute the resulting local degrees of freedom
239 auto dofvec(ref_shape_fns->NodalValuesToDofs(
240 Eigen::Map<Eigen::Matrix<scalar_t, 1, Eigen::Dynamic>>(
241 &uvalvec[0], 1, uvalvec.size())));
242
243 // Set the corresponing global degrees of freedom
244 // Note: "Setting", not "adding to"
245 // Number of local shape functions
246 const size_type num_loc_dofs = dofh.NumLocalDofs(*cell);
247 LF_ASSERT_MSG(
248 dofvec.size() == num_loc_dofs,
249 "Size mismatch: " << dofvec.size() << " <-> " << num_loc_dofs);
250 // Fetch global numbers of local shape functions
251 const std::span<const lf::assemble::gdof_idx_t> ldof_gidx(
252 dofh.GlobalDofIndices(*cell));
253 // Insert dof values into the global coefficient vector
254 for (int j = 0; j < num_loc_dofs; ++j) {
255 glob_dofvec[ldof_gidx[j]] = dofvec[j];
256 }
257 }
258 return glob_dofvec;
259}
260
301template <typename SCALAR, typename EDGESELECTOR,
302 mesh::utils::MeshFunction FUNCTION>
303std::vector<std::pair<bool, SCALAR>> InitEssentialConditionFromFunction(
304 const lf::fe::ScalarFESpace<SCALAR> &fes, EDGESELECTOR &&esscondflag,
305 FUNCTION const &g) {
306 // static_assert(isMeshFunction<FUNCTION>, "g must by a MeshFunction object");
307
308 // *** I: Preprocessing ***
309
310 // Underlying mesh
311 const lf::mesh::Mesh &mesh{*fes.Mesh()};
312
313 // Vector for returning flags and fixed values for dofs
314 const size_type num_coeffs = fes.LocGlobMap().NumDofs();
315 std::vector<std::pair<bool, SCALAR>> flag_val_vec(num_coeffs,
316 {false, SCALAR{}});
317
318 // *** II: Local computations ****
319 // Visit all edges of the mesh (codim-1 entities)
320 for (const lf::mesh::Entity *edge : mesh.Entities(1)) {
321 // Check whether the current edge carries dofs to be imposed by the
322 // function g. The decision relies on the predicate `esscondflag`
323 if (esscondflag(*edge)) {
324 // retrieve reference finite element:
325 auto fe = fes.ShapeFunctionLayout(*edge);
326 LF_ASSERT_MSG(fe != nullptr,
327 "No ScalarReferenceFiniteElement for this edge!");
328 auto ref_eval_pts = fe->EvaluationNodes();
329
330 // Evaluate mesh function at several points specified by their
331 // reference coordinates.
332 auto g_vals = g(*edge, ref_eval_pts);
333
334 // Compute degrees of freedom from function values in evaluation
335 // points
336 auto dof_vals = fe->NodalValuesToDofs(
337 Eigen::Map<Eigen::Matrix<SCALAR, 1, Eigen::Dynamic>>(&g_vals[0], 1,
338 g_vals.size()));
339 LF_ASSERT_MSG(dof_vals.size() == fe->NumRefShapeFunctions(),
340 "Mismatch " << dof_vals.size() << " <-> "
341 << fe->NumRefShapeFunctions());
342 LF_ASSERT_MSG(fes.LocGlobMap().NumLocalDofs(*edge) == dof_vals.size(),
343 "Mismatch " << fes.LocGlobMap().NumLocalDofs(*edge)
344 << " <-> " << dof_vals.size());
345 // Fetch indices of global shape functions associated with current
346 // edge
347 auto gdof_indices{fes.LocGlobMap().GlobalDofIndices(*edge)};
348 int k = 0;
349 // Set flags and values; setting, no accumulation here!
350 for (const lf::assemble::gdof_idx_t gdof_idx : gdof_indices) {
351 flag_val_vec[gdof_idx] = {true, dof_vals[k++]};
352 }
353 }
354 }
355 return flag_val_vec;
356}
357
358} // namespace lf::fe
359
360#endif
A general (interface) class for DOF handling, see Lecture Document Paragraph 2.7.4....
Definition dofhandler.h:112
virtual std::span< const gdof_idx_t > GlobalDofIndices(const lf::mesh::Entity &entity) const =0
access to indices of global dof's belonging to an entity
virtual size_type NumLocalDofs(const lf::mesh::Entity &entity) const =0
tells the number of degrees of freedom subordinate/_belonging_ to to an entity
virtual size_type NumDofs() const =0
total number of dof's handled by the object
A Function Object that can be invoked with any arguments and that always returns the value true.
Represents a reference element with all its properties.
Definition ref_el.h:109
static constexpr RefEl kSegment()
Returns the (1-dimensional) reference segment.
Definition ref_el.h:153
constexpr RefEl(RefElType type) noexcept
Create a RefEl from a lf::base::RefElType enum.
Definition ref_el.h:175
static constexpr RefEl kTria()
Returns the reference triangle.
Definition ref_el.h:161
static constexpr RefEl kQuad()
Returns the reference quadrilateral.
Definition ref_el.h:169
Space of scalar valued finite element functions on a Mesh.
virtual ScalarReferenceFiniteElement< SCALAR > const * ShapeFunctionLayout(const lf::mesh::Entity &entity) const =0
access to shape function layout for mesh entities
virtual std::shared_ptr< const lf::mesh::Mesh > Mesh() const =0
acess to underlying mesh
virtual const lf::assemble::DofHandler & LocGlobMap() const =0
access to associated local-to-global map
virtual Eigen::VectorXd IntegrationElement(const Eigen::MatrixXd &local) const =0
The integration element (factor appearing in integral transformation formula, see below) at number of...
Interface class representing a topological entity in a cellular complex
Definition entity.h:42
virtual const geometry::Geometry * Geometry() const =0
Describes the geometry of this entity.
Abstract interface for objects representing a single mesh.
virtual std::span< const Entity *const > Entities(unsigned codim) const =0
All entities of a given codimension.
Check if a given type T is an Eigen::Matrix.
Definition eigen_tools.h:70
A MeshFunction is a function object that can be evaluated at any point on the mesh.
Eigen::Index gdof_idx_t
Collects data structures and algorithms designed for scalar finite element methods primarily meant fo...
Definition fe.h:47
auto IntegrateMeshFunction(const lf::mesh::Mesh &mesh, const MF &mf, const QR_SELECTOR &qr_selector, const ENTITY_PREDICATE &ep=base::PredicateTrue{}, int codim=0) -> mesh::utils::MeshFunctionReturnType< MF >
Integrate a MeshFunction over a mesh (with quadrature rules)
Definition fe_tools.h:110
static const unsigned int kout_l2_qr
Definition fe_tools.h:26
lf::assemble::size_type size_type
Definition fe.h:54
std::vector< std::pair< bool, SCALAR > > InitEssentialConditionFromFunction(const lf::fe::ScalarFESpace< SCALAR > &fes, EDGESELECTOR &&esscondflag, FUNCTION const &g)
Initialization of flags/values for dofs of a Scalar finite element space whose values are imposed by ...
Definition fe_tools.h:303
static const unsigned int kout_l2_rsfvals
Definition fe_tools.h:27
auto NodalProjection(const lf::fe::ScalarFESpace< SCALAR > &fe_space, MF const &u, SELECTOR &&pred=base::PredicateTrue{}) -> Eigen::Vector< decltype(SCALAR{0} *mesh::utils::MeshFunctionReturnType< std::remove_reference_t< MF > >{0}), Eigen::Dynamic >
Computes nodal projection of a mesh function and returns the finite element basis expansion coefficie...
Definition fe_tools.h:198
internal::VectorElement_t< internal::MeshFunctionReturnType_t< R > > MeshFunctionReturnType
Determine the type of objects returned by a MeshFunction.
QuadRule make_QuadRule(base::RefEl ref_el, unsigned degree)
Returns a QuadRule object for the given Reference Element and Degree.