- /**
- * Implement setting dof indices on a cell. TO CHECK ZHAOWEI + LH
- */
- template <int dim, int spacedim, bool level_dof_access>
- static void
- set_dof_indices(
- const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
- const std::vector<types::global_dof_index> &local_dof_indices)
- {
- Assert(accessor.has_children() == false, ExcInternalError());
-
- const unsigned int dofs_per_vertex =
- accessor.get_fe().n_dofs_per_vertex(),
- dofs_per_line = accessor.get_fe().n_dofs_per_line(),
- dofs_per_quad = accessor.get_fe().n_dofs_per_quad(),
- dofs_per_hex = accessor.get_fe().n_dofs_per_hex();
-
- Assert(local_dof_indices.size() == accessor.get_fe().dofs_per_cell,
- ExcInternalError());
-
- unsigned int index = 0;
-
- for (unsigned int vertex = 0;
- vertex < GeometryInfo<dim>::vertices_per_cell;
- ++vertex)
- for (unsigned int d = 0; d < dofs_per_vertex; ++d, ++index)
- accessor.set_vertex_dof_index(vertex,
- d,
- local_dof_indices[index],
- accessor.active_fe_index());
- // now copy dof numbers into the line. for lines in 3d with the
- // wrong orientation, we have already made sure that we're ok
- // by picking the correct vertices (this happens automatically
- // in the vertex() function). however, if the line is in wrong
- // orientation, we look at it in flipped orientation and we
- // will have to adjust the shape function indices that we see
- // to correspond to the correct (cell-local) ordering.
- //
- // of course, if dim<3, then there is nothing to adjust
- for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_cell;
- ++line)
- for (unsigned int d = 0; d < dofs_per_line; ++d, ++index)
- accessor.line(line)->set_dof_index(
- dim < 3 ?
- d :
- accessor.get_fe().adjust_line_dof_index_for_line_orientation(
- d, accessor.line_orientation(line)),
- local_dof_indices[index],
- accessor.active_fe_index());
- // now copy dof numbers into the face. for faces in 3d with the
- // wrong orientation, we have already made sure that we're ok
- // by picking the correct lines and vertices (this happens
- // automatically in the line() and vertex()
- // functions). however, if the face is in wrong orientation,
- // we look at it in flipped orientation and we will have to
- // adjust the shape function indices that we see to correspond
- // to the correct (cell-local) ordering. The same applies, if
- // the face_rotation or face_orientation is non-standard
- //
- // again, if dim<3, then there is nothing to adjust
- for (unsigned int quad = 0; quad < GeometryInfo<dim>::quads_per_cell;
- ++quad)
- for (unsigned int d = 0; d < dofs_per_quad; ++d, ++index)
- accessor.quad(quad)->set_dof_index(
- dim < 3 ?
- d :
- accessor.get_fe().adjust_quad_dof_index_for_face_orientation(
- d,
- accessor.face_orientation(quad),
- accessor.face_flip(quad),
- accessor.face_rotation(quad)),
- local_dof_indices[index],
- accessor.active_fe_index());
- for (unsigned int d = 0; d < dofs_per_hex; ++d, ++index)
- accessor.set_dof_index(d,
- local_dof_indices[index],
- accessor.active_fe_index());
- Assert(index == accessor.get_fe().dofs_per_cell, ExcInternalError());
- }
-
-
-
/**
* Implement setting non-local dof indices on a cell.
*/
{
Assert(accessor.has_children() == false, ExcInternalError());
- const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell;
+ const unsigned int dofs_per_cell = accessor.get_fe().n_dofs_per_cell();
+ const unsigned int n_non_local_dofs =
+ accessor.get_fe().n_non_local_dofs_per_cell();
+ const unsigned int n_local_dofs = dofs_per_cell - n_non_local_dofs;
+
+ AssertDimension(local_non_local_dof_indices.size(), n_non_local_dofs);
+
+ // First the easy case.
+ if (n_non_local_dofs == 0)
+ return;
types::global_dof_index *next_dof_index =
const_cast<types::global_dof_index *>(
get_cache_ptr(accessor.dof_handler,
accessor.present_level,
accessor.present_index,
- dofs_per_cell));
-
- const unsigned int non_local_dofs =
- accessor.get_fe().n_non_local_dofs_per_cell(),
- n_dofs = accessor.get_fe().n_dofs_per_cell();
-
- unsigned int index = 0;
+ dofs_per_cell)) +
+ n_local_dofs;
- for (unsigned int d = 0; d < n_dofs; ++d, ++next_dof_index)
- if (d >= n_dofs - non_local_dofs)
- {
- *next_dof_index = local_non_local_dof_indices[index];
- ++index;
- }
- Assert(index == accessor.get_fe().n_non_local_dofs_per_cell(),
- ExcInternalError());
+ for (unsigned int d = 0; d < n_non_local_dofs; ++d, ++next_dof_index)
+ *next_dof_index = local_non_local_dof_indices[d];
}
void
distribute_mg_dofs();
+ /**
+ * Distribute non local degrees of freedom. The local DoFs need to be
+ * distributed using distribute_dofs() before calling this function.
+ */
+ void
+ distribute_non_local_dofs();
+
/**
* This function returns whether this DoFHandler has DoFs distributed on
* each multigrid level or in other words if distribute_mg_dofs() has been
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_non_local_dof_handler_h
+#define dealii_non_local_dof_handler_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/subscriptor.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+// Forward declarations
+#ifndef DOXYGEN
+template <int dim, int spacedim, bool>
+class DoFCellAccessor;
+#endif
+
+/**
+ * This class is used to enumerate non local degrees of freedom. Its default
+ * implementation does nothing, since in general FiniteElement spaces only
+ * define degrees of freedom on vertices, edges, faces, or cells. There are
+ * cases, however, in which a FiniteElement space may define some non-local
+ * basis functions which are non-zero on a given cell, even if the basis
+ * function cannot be attributed to local subobjects of the given cell.
+ *
+ * In these cases, the DoFHandler class needs to know how to distribute these
+ * extra degrees of freedom, and since it cannot do so on its own, it asks
+ * the FiniteElement to return a NonLocalDoFHandler class that knows how to
+ * enumerate and distribute these non local degrees of freedom.
+ *
+ * This class is intended as a base class for all FiniteElement spaces that
+ * need to enumerate non local degrees of freedom.
+ *
+ * @ingroup dofs
+ */
+template <int dim, int spacedim = dim>
+class NonLocalDoFHandler : public Subscriptor
+{
+public:
+ /**
+ * Make the dimension available in function templates.
+ */
+ static const unsigned int dimension = dim;
+
+ /**
+ * Make the space dimension available in function templates.
+ */
+ static const unsigned int space_dimension = spacedim;
+
+ /**
+ * Standard constructor, not initializing any data.
+ */
+ NonLocalDoFHandler() = default;
+
+ /**
+ * Destructor.
+ */
+ virtual ~NonLocalDoFHandler() = default;
+
+ /**
+ * Copy operator. NonLocalDoFHandler objects may be large and expensive.
+ * They should not be copied, in particular not by accident, but
+ * rather deliberately constructed. As a consequence, this operator
+ * is explicitly removed from the interface of this class.
+ */
+ NonLocalDoFHandler &
+ operator=(const NonLocalDoFHandler &) = delete;
+
+ /**
+ * Return the non local dof indices associated to the current cell, for
+ * active cell iterators.
+ */
+ virtual std::vector<types::global_dof_index>
+ get_non_local_dof_indices(
+ const DoFCellAccessor<dim, spacedim, false> &) const;
+
+ /**
+ * Return the non local dof indices associated to the current cell, for
+ * level cell iterators.
+ */
+ virtual std::vector<types::global_dof_index>
+ get_non_local_dof_indices(const DoFCellAccessor<dim, spacedim, true> &) const;
+
+ /**
+ * Return the number of non local dof indices that are required in addition to
+ * the local ones.
+ */
+ virtual types::global_dof_index
+ n_additional_non_local_dofs() const;
+
+ /**
+ * At any given moment, only one NonLocalDoFHandler can be used. Throw an
+ * exception if two cells define an active FiniteElement for which
+ * get_non_local_dof_handler() return a different object.
+ *
+ * @ingroup Exceptions
+ */
+ DeclException0(ExcDifferentNonLocalDoFHandlers);
+};
+
+
+
+// ----------------------------------------------------------------------------
+template <int dim, int spacedim>
+inline std::vector<types::global_dof_index>
+NonLocalDoFHandler<dim, spacedim>::get_non_local_dof_indices(
+ const DoFCellAccessor<dim, spacedim, true> &) const
+{
+ // By default, we return an empty vector.
+ return std::vector<types::global_dof_index>();
+}
+
+
+
+template <int dim, int spacedim>
+inline std::vector<types::global_dof_index>
+NonLocalDoFHandler<dim, spacedim>::get_non_local_dof_indices(
+ const DoFCellAccessor<dim, spacedim, false> &) const
+{
+ // By default, we return an empty vector.
+ return std::vector<types::global_dof_index>();
+}
+
+
+
+template <int dim, int spacedim>
+inline types::global_dof_index
+NonLocalDoFHandler<dim, spacedim>::n_additional_non_local_dofs() const
+{
+ return 0;
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
#include <deal.II/base/config.h>
+#include <deal.II/dofs/non_local_dof_handler.h>
+
#include <deal.II/fe/block_mask.h>
#include <deal.II/fe/component_mask.h>
#include <deal.II/fe/fe_base.h>
//@}
+ /**
+ * @name Non local dofs support
+ * @{
+ */
+ /**
+ * Return an object that knows how to handle non local dof indices.
+ */
+ virtual std::shared_ptr<const NonLocalDoFHandler<dim, spacedim>>
+ get_non_local_dof_handler() const;
+ //@}
+
/**
* Determine an estimate for the memory consumption (in bytes) of this
* object.
+template <int dim, int spacedim>
+inline std::shared_ptr<const NonLocalDoFHandler<dim, spacedim>>
+FiniteElement<dim, spacedim>::get_non_local_dof_handler() const
+{
+ return std::make_shared<const NonLocalDoFHandler<dim, spacedim>>();
+}
+
+
DEAL_II_NAMESPACE_CLOSE
#endif
&*this->tria) == nullptr)
this->block_info_object.initialize(*this, false, true);
}
+ // now distribute non_local_dofs.
+ distribute_non_local_dofs();
}
+template <int dim, int spacedim>
+void
+DoFHandler<dim, spacedim>::distribute_non_local_dofs()
+{
+ Assert(
+ this->object_dof_indices.size() > 0,
+ ExcMessage(
+ "Distribute active DoFs using distribute_dofs() before calling distribute_non_local_dofs()."));
+
+ auto non_local_dh = begin()->get_fe().get_non_local_dof_handler();
+ for (auto cell : active_cell_iterators())
+ {
+ // Make sure all cells use the same NonLocalDoFHandler.
+ Assert(non_local_dh == cell->get_fe().get_non_local_dof_handler(),
+ ExcMessage(
+ "You are trying to use different NonLocalDoFHandler objects"));
+ cell->set_non_local_dof_indices(
+ non_local_dh->get_non_local_dof_indices(*cell));
+ }
+ this->number_cache.n_global_dofs +=
+ non_local_dh->n_additional_non_local_dofs();
+}
+
+
+
template <int dim, int spacedim>
void
DoFHandler<dim, spacedim>::initialize_local_block_info()
return dpo;
}
+
+/**
+ * Associates non local dofs according to the vertex index.
+ */
+template <int dim>
+class NonLocalQ1DoFHandler : public NonLocalDoFHandler<dim>
+{
+public:
+ virtual std::vector<types::global_dof_index>
+ get_non_local_dof_indices(
+ const DoFCellAccessor<dim, dim, false> &accessor) const override
+ {
+ if (!tria)
+ tria = &(accessor.get_triangulation());
+ std::vector<types::global_dof_index> dofs(accessor.n_vertices());
+ for (unsigned int i = 0; i < dofs.size(); ++i)
+ dofs[i] = accessor.vertex_index(i);
+ return dofs;
+ }
+
+ virtual types::global_dof_index
+ n_additional_non_local_dofs() const override
+ {
+ Assert(tria, ExcInternalError());
+ return tria->n_vertices();
+ }
+
+private:
+ mutable SmartPointer<const Triangulation<dim>> tria = nullptr;
+};
+
+
+
template <int dim>
class FE_Q1_Nonlocal : public FE_Q_Base<TensorProductPolynomials<dim>, dim, dim>
{
public:
- FE_Q1_Nonlocal()
+ FE_Q1_Nonlocal(const std::shared_ptr<NonLocalQ1DoFHandler<dim>> &ptr =
+ std::shared_ptr<NonLocalQ1DoFHandler<dim>>())
: FE_Q_Base<TensorProductPolynomials<dim>, dim, dim>(
TensorProductPolynomials<dim>(
Polynomials::generate_complete_Lagrange_basis(
1,
FiniteElementData<dim>::H1),
std::vector<bool>(1, false))
+ , non_local_dh(ptr ? ptr : std::make_shared<NonLocalQ1DoFHandler<dim>>())
{
this->unit_support_points = QTrapez<dim>().get_points();
}
virtual std::unique_ptr<FiniteElement<dim>>
clone() const override
{
- return std::make_unique<FE_Q1_Nonlocal<dim>>();
+ return std::make_unique<FE_Q1_Nonlocal<dim>>(non_local_dh);
}
virtual std::string
{
return "FE_Q_Nonlocal<dim>";
}
+
+ virtual std::shared_ptr<const NonLocalDoFHandler<dim>>
+ get_non_local_dof_handler() const override
+ {
+ return non_local_dh;
+ }
+
+ virtual void
+ convert_generalized_support_point_values_to_dof_values(
+ const std::vector<Vector<double>> &support_point_values,
+ std::vector<double> & nodal_values) const override
+ {
+ AssertDimension(support_point_values.size(),
+ this->get_unit_support_points().size());
+ AssertDimension(support_point_values.size(), nodal_values.size());
+ AssertDimension(this->n_dofs_per_cell(), nodal_values.size());
+
+ for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
+ {
+ AssertDimension(support_point_values[i].size(), 1);
+ nodal_values[i] = support_point_values[i](0);
+ }
+ }
+
+private:
+ std::shared_ptr<NonLocalQ1DoFHandler<dim>> non_local_dh;
};
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// build an FE_Q1_Nonlocal finite element, distribute dofs on a simple
+// Triangulation, and interpolate a scalar function.
+
+#include <deal.II/base/function_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../non_local/fe_q1_nonlocal.h"
+
+template <int dim>
+void
+test()
+{
+ FE_Q1_Nonlocal<dim> fe;
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(1);
+ DoFHandler<dim> dh(tria);
+ dh.distribute_dofs(fe);
+
+ Tensor<1, dim> ones;
+ for (unsigned int d = 0; d < dim; ++d)
+ ones[d] = 1;
+
+ Vector<double> solution(dh.n_dofs());
+
+ VectorTools::interpolate(dh, Functions::Monomial<dim>(ones), solution);
+
+ deallog << solution << std::endl;
+}
+
+int
+main()
+{
+ initlog();
+
+ test<1>();
+ test<2>();
+ test<3>();
+}
--- /dev/null
+
+DEAL::0.00000 1.00000 0.500000
+DEAL::0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.500000 0.500000 0.250000
+DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.00000 0.500000 0.500000 0.00000 0.00000 0.00000 0.250000 0.250000 0.250000 0.125000
projects / dealii.git / commitdiff