27#include <boost/container/small_vector.hpp>
38template <
int dim,
int spacedim>
50template <
int dim,
int spacedim>
56 const std::array<Point<spacedim>, 2>
vertices{{p1, p2}};
62template <
int dim,
int spacedim>
68 const double tol = 1
e-10;
69 const unsigned int n_points = surrounding_points.size();
75 "There should be as many surrounding points as weights given."));
79 ExcMessage(
"The weights for the individual points should sum to 1!"));
84 boost::container::small_vector<unsigned int, 100> permutation(n_points);
85 std::iota(permutation.begin(), permutation.end(), 0u);
86 std::sort(permutation.begin(),
88 [&weights](
const std::size_t a,
const std::size_t b) {
89 return weights[a] < weights[b];
94 double w = weights[permutation[0]];
96 for (
unsigned int i = 1; i < n_points; ++i)
99 if (
std::abs(weights[permutation[i]] + w) < tol)
102 weight =
w / (weights[permutation[i]] +
w);
106 p = get_intermediate_point(p,
107 surrounding_points[permutation[i]],
112 p = surrounding_points[permutation[i]];
114 w += weights[permutation[i]];
122template <
int dim,
int spacedim>
131 for (
unsigned int row = 0; row < weights.size(0); ++row)
145 const int spacedim = 2;
150 ((p - face->vertex(0)).norm_square() > (p - face->vertex(1)).norm_square() ?
151 -get_tangent_vector(p, face->vertex(0)) :
152 get_tangent_vector(p, face->vertex(1)));
156 return normal / normal.
norm();
166 const int spacedim = 3;
168 const std::array<Point<spacedim>, 4>
vertices{
169 {face->vertex(0), face->vertex(1), face->vertex(2), face->vertex(3)}};
175 std::max(distances[2], distances[3]));
183 double abs_cos_angle = std::numeric_limits<double>::max();
186 for (
unsigned int i = 0; i < 3; ++i)
187 if (distances[i] > 1e-8 * max_distance)
188 for (
unsigned int j = i + 1; j < 4; ++j)
189 if (distances[j] > 1e-8 * max_distance)
192 (distances[i] * distances[j]);
195 if (
std::abs(new_angle) < 0.999 * abs_cos_angle)
197 abs_cos_angle =
std::abs(new_angle);
203 ExcMessage(
"The search for possible directions did not succeed."));
211 std::vector<bool> tested_vertices(
vertices.size(),
false);
212 tested_vertices[first_index] =
true;
213 tested_vertices[second_index] =
true;
218 t1 = get_tangent_vector(p,
vertices[first_index]);
219 t2 = get_tangent_vector(p,
vertices[second_index]);
220 normal = cross_product_3d(t1, t2);
223 if (normal.
norm_square() < 1e4 * std::numeric_limits<double>::epsilon() *
228 std::find(tested_vertices.begin(), tested_vertices.end(),
false);
229 if (first_false == tested_vertices.end())
236 second_index = first_false - tested_vertices.begin();
247 normal.
norm_square() > 1e4 * std::numeric_limits<double>::epsilon() *
250 "Manifold::normal_vector was unable to find a suitable combination "
251 "of vertices to compute a normal on this face. We chose the secants "
252 "that are as orthogonal as possible, but tangents appear to be "
253 "linearly dependent. Check for distorted faces in your triangulation."));
262 if (reference_normal * normal < 0.0)
265 return normal / normal.
norm();
270template <
int dim,
int spacedim>
286 FaceVertexNormals &n)
const
288 n[0] = cross_product_2d(get_tangent_vector(face->vertex(0), face->vertex(1)));
290 -cross_product_2d(get_tangent_vector(face->vertex(1), face->vertex(0)));
292 for (
unsigned int i = 0; i < 2; ++i)
296 "computed normals have "
308 FaceVertexNormals &n)
const
310 n[0] = cross_product_3d(get_tangent_vector(face->vertex(0), face->vertex(1)),
311 get_tangent_vector(face->vertex(0), face->vertex(2)));
313 n[1] = cross_product_3d(get_tangent_vector(face->vertex(1), face->vertex(3)),
314 get_tangent_vector(face->vertex(1), face->vertex(0)));
316 n[2] = cross_product_3d(get_tangent_vector(face->vertex(2), face->vertex(0)),
317 get_tangent_vector(face->vertex(2), face->vertex(3)));
319 n[3] = cross_product_3d(get_tangent_vector(face->vertex(3), face->vertex(2)),
320 get_tangent_vector(face->vertex(3), face->vertex(1)));
322 for (
unsigned int i = 0; i < 4; ++i)
326 "computed normals have "
334template <
int dim,
int spacedim>
338 FaceVertexNormals &n)
const
340 for (
unsigned int v = 0; v < face->reference_cell().n_vertices(); ++v)
342 n[v] = normal_vector(face, face->vertex(v));
349template <
int dim,
int spacedim>
361template <
int dim,
int spacedim>
373template <
int dim,
int spacedim>
383 return get_new_point_on_line(face);
385 return get_new_point_on_quad(face);
393template <
int dim,
int spacedim>
401 return get_new_point_on_line(cell);
403 return get_new_point_on_quad(cell);
405 return get_new_point_on_hex(cell);
479template <
int dim,
int spacedim>
495 const auto points_weights =
503template <
int dim,
int spacedim>
510 const std::array<Point<spacedim>, 2> points{{x1, x2}};
514 return (neighbor_point - x1) /
epsilon;
524 normalized_alternating_product(
const Tensor<1, 3> (&)[1])
536 normalized_alternating_product(
const Tensor<1, 3> (&basis_vectors)[2])
538 Tensor<1, 3> tmp = cross_product_3d(basis_vectors[0], basis_vectors[1]);
539 return tmp / tmp.
norm();
547template <
int dim,
int spacedim>
550 const double tolerance)
551 : periodicity(periodicity)
552 , tolerance(tolerance)
557template <
int dim,
int spacedim>
558std::unique_ptr<Manifold<dim, spacedim>>
561 return std::make_unique<FlatManifold<dim, spacedim>>(periodicity, tolerance);
566template <
int dim,
int spacedim>
574 ExcMessage(
"The weights for the new point should sum to 1!"));
581 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
582 p += surrounding_points[i] * weights[i];
588 for (
unsigned int d = 0;
d < spacedim; ++
d)
589 if (periodicity[d] > 0)
590 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
592 minP[
d] =
std::min(minP[d], surrounding_points[i][d]);
593 Assert((surrounding_points[i][d] <
594 periodicity[d] + tolerance * periodicity[d]) ||
595 (surrounding_points[i][d] >=
596 -tolerance * periodicity[d]),
597 ExcPeriodicBox(d, surrounding_points[i], periodicity[d]));
602 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
605 for (
unsigned int d = 0;
d < spacedim; ++
d)
606 if (periodicity[d] > 0)
608 ((surrounding_points[i][
d] - minP[
d]) > periodicity[d] / 2.0 ?
612 p += (surrounding_points[i] + dp) * weights[i];
616 for (
unsigned int d = 0;
d < spacedim; ++
d)
617 if (periodicity[d] > 0)
619 p[
d] += periodicity[
d];
622 return project_to_manifold(surrounding_points, p);
627template <
int dim,
int spacedim>
635 if (weights.size(0) == 0)
639 const std::size_t n_points = surrounding_points.size();
646 for (
unsigned int row = 0; row < weights.size(0); ++row)
648 &weights(row, 0) + n_points,
651 ExcMessage(
"The weights for each of the points should sum to "
654 constexpr std::size_t n_lanes =
657 const std::size_t n_regular_cols = (n_points / n_lanes) * n_lanes;
658 for (
unsigned int row = 0; row < weights.size(0); row += n_lanes)
660 std::array<unsigned int, n_lanes> offsets;
663 for (std::size_t i = 0; i < n_lanes; ++i)
665 std::min<unsigned int>((row + i) * n_points,
666 (weights.size(0) - 1) * n_points);
668 for (std::size_t col = 0; col < n_regular_cols; col += n_lanes)
670 std::array<VectorizedArrayType, n_lanes> vectorized_weights;
672 &weights(0, 0) + col,
674 vectorized_weights.data());
675 for (std::size_t i = 0; i < n_lanes; ++i)
676 point += vectorized_weights[i] * surrounding_points[col + i];
678 for (std::size_t col = n_regular_cols; col < n_points; ++col)
680 VectorizedArrayType vectorized_weights;
681 vectorized_weights.gather(&weights(0, 0) + col, offsets.data());
682 point += vectorized_weights * surrounding_points[col];
684 for (
unsigned int r = row;
685 r < std::min<unsigned int>(weights.size(0), row + n_lanes);
689 for (
unsigned int d = 0;
d < spacedim; ++
d)
690 new_points[r][d] = point[d][r - row];
692 project_to_manifold(surrounding_points, new_points[r]);
697 for (
unsigned int row = 0; row < weights.size(0); ++row)
699 get_new_point(surrounding_points,
705template <
int dim,
int spacedim>
716template <
int dim,
int spacedim>
725template <
int dim,
int spacedim>
736 for (
unsigned int d = 0;
d < spacedim; ++
d)
737 if (periodicity[d] > tolerance)
739 if (direction[d] < -periodicity[d] / 2)
740 direction[
d] += periodicity[
d];
741 else if (direction[d] > periodicity[d] / 2)
742 direction[
d] -= periodicity[
d];
789 const Tensor<1, 2> tangent = face->vertex(1) - face->vertex(0);
793 face_vertex_normals[vertex] =
Tensor<1, 2>({tangent[1], -tangent[0]});
817 static const unsigned int neighboring_vertices[4][2] = {{1, 2},
821 for (
unsigned int vertex = 0; vertex < vertices_per_face; ++vertex)
826 face->vertex(neighboring_vertices[vertex][0]) - face->vertex(vertex),
827 face->vertex(neighboring_vertices[vertex][1]) - face->vertex(vertex)};
831 face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
883template <
int dim,
int spacedim>
921 const unsigned int facedim = dim - 1;
925 const auto face_reference_cell = face->reference_cell();
928 (face_reference_cell == ReferenceCells::get_hypercube<facedim>()))
930 for (
unsigned int i = 0; i < facedim; ++i)
935 for (
unsigned int i = 0; i < facedim; ++i)
939 const double eps = 1
e-12;
941 unsigned int iteration = 0;
945 for (
const unsigned int v : face->vertex_indices())
947 face->vertex(v) * face_reference_cell.d_linear_shape_function(xi, v);
949 for (
unsigned int i = 0; i < facedim; ++i)
955 face_reference_cell.d_linear_shape_function_gradient(xi, v)[i];
959 for (
unsigned int i = 0; i < facedim; ++i)
960 for (
unsigned int j = 0; j < spacedim; ++j)
961 J[i] += grad_F[i][j] * (F - p)[j];
964 for (
unsigned int i = 0; i < facedim; ++i)
965 for (
unsigned int j = 0; j < facedim; ++j)
966 for (
unsigned int k = 0; k < spacedim; ++k)
967 H[i][j] += grad_F[i][k] * grad_F[j][k];
975 "The Newton iteration to find the reference point "
976 "did not converge in 10 iterations. Do you have a "
977 "deformed cell? (See the glossary for a definition "
978 "of what a deformed cell is. You may want to output "
979 "the vertices of your cell."));
988 const double normalized_delta_world = (F - p).norm() / face->diameter();
990 if (delta_xi.
norm() < eps || normalized_delta_world < eps)
998 return internal::normalized_alternating_product(grad_F);
1003template <
int dim,
int spacedim,
int chartdim>
1006 : sub_manifold(periodicity)
1011template <
int dim,
int spacedim,
int chartdim>
1016 const double w)
const
1018 const std::array<Point<spacedim>, 2> points{{p1, p2}};
1019 const std::array<double, 2> weights{{1. - w, w}};
1025template <
int dim,
int spacedim,
int chartdim>
1031 const std::size_t n_points = surrounding_points.size();
1033 boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1035 for (
unsigned int i = 0; i < n_points; ++i)
1036 chart_points[i] = pull_back(surrounding_points[i]);
1039 sub_manifold.get_new_point(chart_points, weights);
1041 return push_forward(p_chart);
1046template <
int dim,
int spacedim,
int chartdim>
1056 const std::size_t n_points = surrounding_points.size();
1058 boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1059 for (std::size_t i = 0; i < n_points; ++i)
1060 chart_points[i] = pull_back(surrounding_points[i]);
1062 boost::container::small_vector<Point<chartdim>, 200> new_points_on_chart(
1064 sub_manifold.get_new_points(chart_points, weights, new_points_on_chart);
1066 for (std::size_t row = 0; row < weights.size(0); ++row)
1067 new_points[row] = push_forward(new_points_on_chart[row]);
1072template <
int dim,
int spacedim,
int chartdim>
1085template <
int dim,
int spacedim,
int chartdim>
1092 push_forward_gradient(pull_back(x1));
1101 1e-12 * F_prime.
norm(),
1103 "The derivative of a chart function must not be singular."));
1106 sub_manifold.get_tangent_vector(pull_back(x1), pull_back(x2));
1109 for (
unsigned int i = 0; i < spacedim; ++i)
1110 result[i] += F_prime[i] * delta;
1117template <
int dim,
int spacedim,
int chartdim>
1121 return sub_manifold.get_periodicity();
1125#include "manifold.inst"
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
ChartManifold(const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >())
virtual void get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim > > new_points) const override
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const
const Tensor< 1, chartdim > & get_periodicity() const
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Manifold< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const override
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim > > &points, const Point< spacedim > &candidate) const override
const Tensor< 1, spacedim > & get_periodicity() const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
FlatManifold(const Tensor< 1, spacedim > &periodicity=Tensor< 1, spacedim >(), const double tolerance=1e-10)
virtual void get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim > > new_points) const override
virtual Point< spacedim > get_new_point_on_hex(const typename Triangulation< dim, spacedim >::hex_iterator &hex) const
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim > > &surrounding_points, const Point< spacedim > &candidate) const
virtual void get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim > > new_points) const
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Point< spacedim > get_new_point_on_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
numbers::NumberTraits< Number >::real_type norm() const
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcEmptyObject()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcPureFunctionCalled()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename IteratorSelector::hex_iterator hex_iterator
typename IteratorSelector::quad_iterator quad_iterator
typename IteratorSelector::line_iterator line_iterator
#define DEAL_II_NOT_IMPLEMENTED()
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
std::pair< std::array< Point< MeshIteratorType::AccessorType::space_dimension >, n_default_points_per_cell< MeshIteratorType >()>, std::array< double, n_default_points_per_cell< MeshIteratorType >()> > get_default_points_and_weights(const MeshIteratorType &iterator, const bool with_interpolation=false)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
constexpr const ReferenceCell Line
static const unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
void vectorized_load_and_transpose(const unsigned int n_entries, const Number *in, const unsigned int *offsets, VectorizedArray< Number, width > *out)