28 #include <boost/container/small_vector.hpp>
37 template <
int dim,
int spacedim>
49 template <
int dim,
int spacedim>
55 const std::array<Point<spacedim>, 2>
vertices{{p1, p2}};
57 w * p2 + (1 -
w) * p1);
62 template <
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."));
77 Assert(std::abs(std::accumulate(weights.
begin(), weights.
end(), 0.0) - 1.0) <
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);
104 if (std::abs(weight) > 1
e-14)
106 p = get_intermediate_point(p,
107 surrounding_points[permutation[i]],
112 p = surrounding_points[permutation[i]];
114 w += weights[permutation[i]];
122 template <
int dim,
int spacedim>
131 for (
unsigned int row = 0; row < weights.size(0); ++row)
135 surrounding_points.end()),
147 const int spacedim = 2;
152 ((p - face->vertex(0)).norm_square() > (p - face->vertex(1)).norm_square() ?
153 -get_tangent_vector(p, face->vertex(0)) :
154 get_tangent_vector(p, face->vertex(1)));
158 return normal / normal.
norm();
168 const int spacedim = 3;
170 const std::array<Point<spacedim>, 4>
vertices{
171 {face->vertex(0), face->vertex(1), face->vertex(2), face->vertex(3)}};
177 std::max(distances[2], distances[3]));
188 for (
unsigned int i = 0; i < 3; ++i)
189 if (distances[i] > 1
e-8 * max_distance)
190 for (
unsigned int j = i + 1; j < 4; ++j)
191 if (distances[j] > 1
e-8 * max_distance)
194 (distances[i] * distances[j]);
197 if (std::abs(new_angle) < 0.999 * abs_cos_angle)
199 abs_cos_angle = std::abs(new_angle);
205 ExcMessage(
"The search for possible directions did not succeed."));
213 std::vector<bool> tested_vertices(
vertices.size(),
false);
214 tested_vertices[first_index] =
true;
215 tested_vertices[second_index] =
true;
220 t1 = get_tangent_vector(p,
vertices[first_index]);
221 t2 = get_tangent_vector(p,
vertices[second_index]);
222 normal = cross_product_3d(t1, t2);
230 std::find(tested_vertices.begin(), tested_vertices.end(),
false);
231 if (first_false == tested_vertices.end())
238 second_index = first_false - tested_vertices.begin();
252 "Manifold::normal_vector was unable to find a suitable combination "
253 "of vertices to compute a normal on this face. We chose the secants "
254 "that are as orthogonal as possible, but tangents appear to be "
255 "linearly dependent. Check for distorted faces in your triangulation."));
264 if (reference_normal * normal < 0.0)
267 return normal / normal.
norm();
272 template <
int dim,
int spacedim>
290 n[0] = cross_product_2d(get_tangent_vector(face->vertex(0), face->vertex(1)));
292 -cross_product_2d(get_tangent_vector(face->vertex(1), face->vertex(0)));
294 for (
unsigned int i = 0; i < 2; ++i)
298 "computed normals have "
312 n[0] = cross_product_3d(get_tangent_vector(face->vertex(0), face->vertex(1)),
313 get_tangent_vector(face->vertex(0), face->vertex(2)));
315 n[1] = cross_product_3d(get_tangent_vector(face->vertex(1), face->vertex(3)),
316 get_tangent_vector(face->vertex(1), face->vertex(0)));
318 n[2] = cross_product_3d(get_tangent_vector(face->vertex(2), face->vertex(0)),
319 get_tangent_vector(face->vertex(2), face->vertex(3)));
321 n[3] = cross_product_3d(get_tangent_vector(face->vertex(3), face->vertex(2)),
322 get_tangent_vector(face->vertex(3), face->vertex(1)));
324 for (
unsigned int i = 0; i < 4; ++i)
328 "computed normals have "
336 template <
int dim,
int spacedim>
342 for (
unsigned int v = 0; v < face->reference_cell().n_vertices(); ++v)
344 n[v] = normal_vector(face, face->vertex(v));
351 template <
int dim,
int spacedim>
358 points_weights.first.end()),
360 points_weights.second.end()));
365 template <
int dim,
int spacedim>
372 points_weights.first.end()),
374 points_weights.second.end()));
379 template <
int dim,
int spacedim>
389 return get_new_point_on_line(face);
391 return get_new_point_on_quad(face);
399 template <
int dim,
int spacedim>
407 return get_new_point_on_line(cell);
409 return get_new_point_on_quad(cell);
411 return get_new_point_on_hex(cell);
485 template <
int dim,
int spacedim>
501 const auto points_weights =
504 points_weights.first.end()),
506 points_weights.second.end()));
511 template <
int dim,
int spacedim>
518 const std::array<Point<spacedim>, 2> points{{x1, x2}};
523 return (neighbor_point - x1) /
epsilon;
533 normalized_alternating_product(
const Tensor<1, 3> (&)[1])
545 normalized_alternating_product(
const Tensor<1, 3> (&basis_vectors)[2])
547 Tensor<1, 3> tmp = cross_product_3d(basis_vectors[0], basis_vectors[1]);
548 return tmp / tmp.
norm();
554 template <
int dim,
int spacedim>
557 const double tolerance)
558 : periodicity(periodicity)
559 , tolerance(tolerance)
564 template <
int dim,
int spacedim>
565 std::unique_ptr<Manifold<dim, spacedim>>
573 template <
int dim,
int spacedim>
579 Assert(std::abs(std::accumulate(weights.
begin(), weights.
end(), 0.0) - 1.0) <
581 ExcMessage(
"The weights for the new point should sum to 1!"));
588 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
589 p += surrounding_points[i] * weights[i];
595 for (
unsigned int d = 0;
d < spacedim; ++
d)
597 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
599 minP[
d] =
std::min(minP[
d], surrounding_points[i][
d]);
600 Assert((surrounding_points[i][
d] <
602 (surrounding_points[i][
d] >=
609 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
612 for (
unsigned int d = 0;
d < spacedim; ++
d)
619 p += (surrounding_points[i] + dp) * weights[i];
623 for (
unsigned int d = 0;
d < spacedim; ++
d)
634 template <
int dim,
int spacedim>
642 if (weights.size(0) == 0)
646 const std::size_t n_points = surrounding_points.size();
653 for (
unsigned int row = 0; row < weights.size(0); ++row)
654 Assert(std::abs(std::accumulate(&weights(row, 0),
655 &weights(row, 0) + n_points,
658 ExcMessage(
"The weights for each of the points should sum to "
661 constexpr std::size_t n_lanes =
664 const std::size_t n_regular_cols = (n_points / n_lanes) * n_lanes;
665 for (
unsigned int row = 0; row < weights.size(0); row += n_lanes)
667 std::array<unsigned int, n_lanes> offsets;
670 for (std::size_t i = 0; i < n_lanes; ++i)
672 std::min<unsigned int>((row + i) * n_points,
673 (weights.size(0) - 1) * n_points);
675 for (std::size_t col = 0; col < n_regular_cols; col += n_lanes)
677 std::array<VectorizedArrayType, n_lanes> vectorized_weights;
679 &weights(0, 0) + col,
681 vectorized_weights.data());
682 for (std::size_t i = 0; i < n_lanes; ++i)
683 point += vectorized_weights[i] * surrounding_points[col + i];
685 for (std::size_t col = n_regular_cols; col < n_points; ++col)
687 VectorizedArrayType vectorized_weights;
688 vectorized_weights.gather(&weights(0, 0) + col, offsets.data());
689 point += vectorized_weights * surrounding_points[col];
691 for (
unsigned int r = row;
692 r < std::min<unsigned int>(weights.size(0), row + n_lanes);
696 for (
unsigned int d = 0;
d < spacedim; ++
d)
697 new_points[r][
d] =
point[
d][r - row];
704 for (
unsigned int row = 0; row < weights.size(0); ++row)
712 template <
int dim,
int spacedim>
723 template <
int dim,
int spacedim>
732 template <
int dim,
int spacedim>
743 for (
unsigned int d = 0;
d < spacedim; ++
d)
796 const Tensor<1, 2> tangent = face->vertex(1) - face->vertex(0);
800 face_vertex_normals[vertex] =
Tensor<1, 2>({tangent[1], -tangent[0]});
824 static const unsigned int neighboring_vertices[4][2] = {{1, 2},
828 for (
unsigned int vertex = 0; vertex < vertices_per_face; ++vertex)
833 face->vertex(neighboring_vertices[vertex][0]) - face->vertex(vertex),
834 face->vertex(neighboring_vertices[vertex][1]) - face->vertex(vertex)};
838 face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
890 template <
int dim,
int spacedim>
926 const unsigned int facedim = dim - 1;
930 const auto face_reference_cell = face->reference_cell();
932 if (face_reference_cell == ReferenceCells::get_hypercube<facedim>())
934 for (
unsigned int i = 0; i < facedim; ++i)
939 for (
unsigned int i = 0; i < facedim; ++i)
943 const double eps = 1
e-12;
945 unsigned int iteration = 0;
949 for (
const unsigned int v : face->vertex_indices())
951 face->vertex(v) * face_reference_cell.d_linear_shape_function(xi, v);
953 for (
unsigned int i = 0; i < facedim; ++i)
956 for (
const unsigned int v : face->vertex_indices())
959 face_reference_cell.d_linear_shape_function_gradient(xi, v)[i];
963 for (
unsigned int i = 0; i < facedim; ++i)
964 for (
unsigned int j = 0; j < spacedim; ++j)
965 J[i] += grad_F[i][j] * (
F - p)[j];
968 for (
unsigned int i = 0; i < facedim; ++i)
969 for (
unsigned int j = 0; j < facedim; ++j)
970 for (
unsigned int k = 0; k < spacedim; ++k)
971 H[i][j] += grad_F[i][k] * grad_F[j][k];
979 "The Newton iteration to find the reference point "
980 "did not converge in 10 iterations. Do you have a "
981 "deformed cell? (See the glossary for a definition "
982 "of what a deformed cell is. You may want to output "
983 "the vertices of your cell."));
992 const double normalized_delta_world = (
F - p).
norm() / face->diameter();
994 if (delta_xi.
norm() <
eps || normalized_delta_world <
eps)
1002 return internal::normalized_alternating_product(grad_F);
1007 template <
int dim,
int spacedim,
int chartdim>
1010 : sub_manifold(periodicity)
1015 template <
int dim,
int spacedim,
int chartdim>
1020 const double w)
const
1022 const std::array<Point<spacedim>, 2> points{{p1, p2}};
1023 const std::array<double, 2> weights{{1. -
w,
w}};
1030 template <
int dim,
int spacedim,
int chartdim>
1036 const std::size_t n_points = surrounding_points.size();
1038 boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1040 for (
unsigned int i = 0; i < n_points; ++i)
1041 chart_points[i] =
pull_back(surrounding_points[i]);
1051 template <
int dim,
int spacedim,
int chartdim>
1061 const std::size_t n_points = surrounding_points.size();
1063 boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1064 for (std::size_t i = 0; i < n_points; ++i)
1065 chart_points[i] =
pull_back(surrounding_points[i]);
1067 boost::container::small_vector<Point<chartdim>, 200> new_points_on_chart(
1069 sub_manifold.get_new_points(
1072 make_array_view(new_points_on_chart.begin(), new_points_on_chart.end()));
1074 for (std::size_t row = 0; row < weights.size(0); ++row)
1075 new_points[row] =
push_forward(new_points_on_chart[row]);
1080 template <
int dim,
int spacedim,
int chartdim>
1093 template <
int dim,
int spacedim,
int chartdim>
1109 1
e-12 * F_prime.
norm(),
1111 "The derivative of a chart function must not be singular."));
1117 for (
unsigned int i = 0; i < spacedim; ++i)
1118 result[i] += F_prime[i] * delta;
1125 template <
int dim,
int spacedim,
int chartdim>
1129 return sub_manifold.get_periodicity();
1133 #include "manifold.inst"
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
ChartManifold(const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >())
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 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_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
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 void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Manifold< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const override
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
const Tensor< 1, spacedim > & get_periodicity() const
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &points, const Point< spacedim > &candidate) const override
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
const Tensor< 1, spacedim > periodicity
FlatManifold(const Tensor< 1, spacedim > &periodicity=Tensor< 1, spacedim >(), const double tolerance=1e-10)
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &surrounding_points, const Point< spacedim > &candidate) const
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
virtual Point< spacedim > get_new_point_on_hex(const typename Triangulation< dim, spacedim >::hex_iterator &hex) const
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
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 Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const
Point< spacedim > get_new_point_on_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcPureFunctionCalled()
static ::ExceptionBase & ExcPeriodicBox(int arg1, Point< spacedim > arg2, double arg3)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename IteratorSelector::hex_iterator hex_iterator
typename IteratorSelector::quad_iterator quad_iterator
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< dim > push_forward(const TopoDS_Shape &in_shape, const double u, const double v)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
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)
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)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
constexpr const ReferenceCell Line
static const unsigned int invalid_unsigned_int
constexpr DEAL_II_HOST 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)