* We collect here some helper functions used in the
* Manifold<dim,spacedim> classes.
*/
-namespace Manifolds
+namespace Manifolds
{
/**
* Given a hex iterator, construct a quadrature with the Laplace
* centers and face centers, which can be called when creating
* middle vertices in the manifold routines.
*/
- Quadrature<3>
- get_default_quadrature(const TriaIterator<CellAccessor<3, 3> >& hex);
-
+ Quadrature<3>
+ get_default_quadrature(const TriaIterator<CellAccessor<3, 3> > &hex);
+
/**
* Given a general mesh iterator, construct a quadrature with the
* Laplace weights or with uniform weights according the parameter
* when creating new vertices in the manifold routines.
*/
template <typename OBJECT, int spacedim>
- Quadrature<spacedim>
- get_default_quadrature(const OBJECT& obj, bool with_laplace = false);
+ Quadrature<spacedim>
+ get_default_quadrature(const OBJECT &obj, bool with_laplace = false);
}
* overload project_to_manifold() for simple situations.
*
* Should a finer control be necessary, then get_new_point() can be
- * overloaded.
+ * overloaded.
*
* FlatManifold is the specialization from which StraigthBoundary is
* derived, where the project_to_manifold() function is the identity.
/**
* Given a point which lies close to the given manifold, it modifies
* it and projects it to manifold itself.
- *
+ *
* This class is used by the default implementation of the function
* get_new_point(). It should be made pure virtual, but for
* historical reason, derived classes like Boundary<dim, spacedim>
* argument.
*/
FlatManifold (const Point<spacedim> periodicity=Point<spacedim>(),
- const double tolerance=1e-10);
+ const double tolerance=1e-10);
/**
* Let the new point be the average sum of surrounding vertices.
*/
virtual Point<spacedim>
get_new_point(const Quadrature<spacedim> &quad) const;
-
-
+
+
/**
* Project to FlatManifold. This is the identity function for flat,
* Euclidean spaces. Note however that this function can be
* similar (if not identical) to the one implemented in this class.
*/
virtual
- Point<spacedim> project_to_manifold (const std::vector<Point<spacedim> > &points,
- const Point<spacedim> &candidate) const;
+ Point<spacedim> project_to_manifold (const std::vector<Point<spacedim> > &points,
+ const Point<spacedim> &candidate) const;
protected:
/**
* Tolerance. This tolerance is used to compute distances in double
const Point<spacedim> periodicity;
DeclException4(ExcPeriodicBox, int, Point<spacedim>, Point<spacedim>, double,
- << "The component number " << arg1 << " of the point [ " << arg2
- << " ] is not in the interval [ " << -arg4
- << ", " << arg3[arg4] << "), bailing out.");
+ << "The component number " << arg1 << " of the point [ " << arg2
+ << " ] is not in the interval [ " << -arg4
+ << ", " << arg3[arg4] << "), bailing out.");
};
*
* Notice that the dimenisions #chartdim and #spacedim can be
* arbitrary, as long as the transformation from $\mathcal{B}$ to
- * $\mathcal{M}$ is invertible.
+ * $\mathcal{M}$ is invertible.
*
* @ingroup manifold
*
* these two points are at distance 2*eps and not (2*pi-eps)
*/
ManifoldChart(const Point<chartdim> periodicity=Point<chartdim>());
-
+
/**
* Destructor. Does nothing here, but needs to be declared to make
* it virtual.
* Refer to the general documentation of this class and the
* documentation of the base class for more information.
*/
- virtual Point<spacedim>
- get_new_point(const Quadrature<spacedim> &quad) const;
-
- /**
- * Pull back the given point in spacedim to the Euclidean chartdim
- * dimensional space.
- *
- * Refer to the general documentation of this class for more
- * information.
- */
- virtual Point<chartdim>
- pull_back(const Point<spacedim> &space_point) const = 0;
-
- /**
- * Given a point in the chartdim dimensional Euclidean space, this
- * method returns a point on the manifold embedded in the spacedim
- * Euclidean space.
- *
- * Refer to the general documentation of this class for more
- * information.
- */
- virtual Point<spacedim>
- push_forward(const Point<chartdim> &chart_point) const = 0;
-
- private:
+ virtual Point<spacedim>
+ get_new_point(const Quadrature<spacedim> &quad) const;
+
+ /**
+ * Pull back the given point in spacedim to the Euclidean chartdim
+ * dimensional space.
+ *
+ * Refer to the general documentation of this class for more
+ * information.
+ */
+ virtual Point<chartdim>
+ pull_back(const Point<spacedim> &space_point) const = 0;
+
+ /**
+ * Given a point in the chartdim dimensional Euclidean space, this
+ * method returns a point on the manifold embedded in the spacedim
+ * Euclidean space.
+ *
+ * Refer to the general documentation of this class for more
+ * information.
+ */
+ virtual Point<spacedim>
+ push_forward(const Point<chartdim> &chart_point) const = 0;
+
+private:
/**
* The sub_manifold object is used to compute the average of the
- * points in the chart coordinates system.
+ * points in the chart coordinates system.
*/
const FlatManifold<dim,chartdim> sub_manifold;
};
/*---Templated functions---*/
-namespace Manifolds {
-
+namespace Manifolds
+{
+
template <typename OBJECT, int spacedim>
- Quadrature<spacedim>
- get_default_quadrature(const OBJECT& obj, bool with_laplace)
+ Quadrature<spacedim>
+ get_default_quadrature(const OBJECT &obj, bool with_laplace)
{
std::vector<Point<spacedim> > sp;
std::vector<double> wp;
const int dim = OBJECT::AccessorType::structure_dimension;
-
+
// note that the exact weights are chosen such as to minimize the
// distortion of the four new quads from the optimal shape; their
// derivation and values is copied over from the
// @p{MappingQ::set_laplace_on_vector} function
AssertDimension(spacedim, OBJECT::AccessorType::space_dimension);
- switch(dim)
+ switch (dim)
{
- case 1:
- sp.resize(2);
- wp.resize(2);
- sp[0] = obj->vertex(0); wp[0] = .5;
- sp[1] = obj->vertex(1); wp[1] = .5;
- break;
- case 2:
- sp.resize(8);
- wp.resize(8);
- sp[0] = obj->vertex(0);
- sp[1] = obj->vertex(1);
- sp[2] = obj->vertex(2);
- sp[3] = obj->vertex(3);
-
- sp[4] = obj->line(0)->has_children() ?
- obj->line(0)->child(0)->vertex(1) :
- obj->line(0)->get_manifold().get_new_point_on_line(obj->line(0));
- sp[5] = obj->line(1)->has_children() ?
- obj->line(1)->child(0)->vertex(1) :
- obj->line(1)->get_manifold().get_new_point_on_line(obj->line(1));
- sp[6] = obj->line(2)->has_children() ?
- obj->line(2)->child(0)->vertex(1) :
- obj->line(2)->get_manifold().get_new_point_on_line(obj->line(2));
- sp[7] = obj->line(3)->has_children() ?
- obj->line(3)->child(0)->vertex(1) :
- obj->line(3)->get_manifold().get_new_point_on_line(obj->line(3));
- if(with_laplace)
- {
- std::fill(wp.begin(), wp.begin()+4, 1.0/16.0);
- std::fill(wp.begin()+4, wp.end(), 3.0/16.0);
- }
- else
- std::fill(wp.begin(), wp.end(), 1.0/8.0);
- break;
- default:
- Assert(false, ExcInternalError());
- break;
+ case 1:
+ sp.resize(2);
+ wp.resize(2);
+ sp[0] = obj->vertex(0);
+ wp[0] = .5;
+ sp[1] = obj->vertex(1);
+ wp[1] = .5;
+ break;
+ case 2:
+ sp.resize(8);
+ wp.resize(8);
+ sp[0] = obj->vertex(0);
+ sp[1] = obj->vertex(1);
+ sp[2] = obj->vertex(2);
+ sp[3] = obj->vertex(3);
+
+ sp[4] = obj->line(0)->has_children() ?
+ obj->line(0)->child(0)->vertex(1) :
+ obj->line(0)->get_manifold().get_new_point_on_line(obj->line(0));
+ sp[5] = obj->line(1)->has_children() ?
+ obj->line(1)->child(0)->vertex(1) :
+ obj->line(1)->get_manifold().get_new_point_on_line(obj->line(1));
+ sp[6] = obj->line(2)->has_children() ?
+ obj->line(2)->child(0)->vertex(1) :
+ obj->line(2)->get_manifold().get_new_point_on_line(obj->line(2));
+ sp[7] = obj->line(3)->has_children() ?
+ obj->line(3)->child(0)->vertex(1) :
+ obj->line(3)->get_manifold().get_new_point_on_line(obj->line(3));
+ if (with_laplace)
+ {
+ std::fill(wp.begin(), wp.begin()+4, 1.0/16.0);
+ std::fill(wp.begin()+4, wp.end(), 3.0/16.0);
+ }
+ else
+ std::fill(wp.begin(), wp.end(), 1.0/8.0);
+ break;
+ default:
+ Assert(false, ExcInternalError());
+ break;
}
return Quadrature<spacedim>(sp,wp);
}
DEAL_II_NAMESPACE_OPEN
/**
- * Manifold description for a spherical space coordinate system.
+ * Manifold description for a spherical space coordinate system.
*
* You can use this Manifold object to describe any sphere, circle,
* hypersphere or hyperdisc in two or three dimensions, both as a
* be attached to cells containing the center of the coordinate
* system. This point is a singular point of the coordinate
* transformation, and there taking averages does not make any sense.
- *
+ *
* @ingroup manifold
*
* @author Luca Heltai, 2014
*/
template <int dim, int spacedim>
-class SphericalManifold : public ManifoldChart<dim, spacedim, spacedim>
+class SphericalManifold : public ManifoldChart<dim, spacedim, spacedim>
{
public:
/**
*/
virtual Point<spacedim>
get_new_point(const Quadrature<spacedim> &quad) const;
-
+
/**
* The center of the spherical coordinate system.
*/
const Point<spacedim> center;
private:
-
+
/** Helper function which returns the periodicity associated with
this coordinate system, according to dim, chartdim, and
spacedim. */
* This class was developed to be used in conjunction with the @p
* cylinder or @p cylinder_shell functions of GridGenerator. This
* function will throw an exception whenever spacedim is not equal to
- * three.
+ * three.
*
* @ingroup manifold
*
* z-axis, respectively. The tolerance value is used to determine
* if a point is on the axis.
*/
- CylindricalManifold (const unsigned int axis = 0,
- const double tolerance = 1e-10);
+ CylindricalManifold (const unsigned int axis = 0,
+ const double tolerance = 1e-10);
/**
* Constructor. If constructed with this constructor, the manifold
* #direction and goes through the given #point_on_axis. The
* direction may be arbitrarily scaled, and the given point may be
* any point on the axis. The tolerance value is used to determine
- * if a point is on the axis.
+ * if a point is on the axis.
*/
CylindricalManifold (const Point<spacedim> &direction,
- const Point<spacedim> &point_on_axis,
- const double tolerance = 1e-10);
+ const Point<spacedim> &point_on_axis,
+ const double tolerance = 1e-10);
- /**
- * Compute new points on the CylindricalManifold. See the documentation
- * of the base class for a detailed description of what this
- * function does.
- */
+ /**
+ * Compute new points on the CylindricalManifold. See the documentation
+ * of the base class for a detailed description of what this
+ * function does.
+ */
virtual Point<spacedim>
get_new_point(const Quadrature<spacedim> &quad) const;
* Helper FlatManifold to compute temptative midpoints.
*/
FlatManifold<dim,spacedim> flat_manifold;
-
+
/**
* Relative tolerance to measure zero distances.
*/
* boundary description has been explicitly
* set using set_manifold().
*/
- static const StraightBoundary<dim,spacedim> straight_boundary;
+ static const StraightBoundary<dim,spacedim> straight_boundary;
/**
* Declare some symbolic names
virtual
void
create_notification (const Triangulation<dim, spacedim> &tria);
- };
+ };
/**
* A structure that is used as an
void set_manifold (const types::manifold_id number,
const Manifold<dim,spacedim> &manifold_object);
-
+
/**
* Reset those parts of the triangulation with the given manifold_id
* to use a FlatManifold object. This is the default state of a
* @ingroup manifold
*
* @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
- */
+ */
const Manifold<dim,spacedim> &get_manifold (const types::manifold_id number) const;
/**
* @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
*/
std::vector<types::boundary_id> get_boundary_indicators() const;
-
+
/**
* Returns a vector containing all manifold indicators assigned to
* the objects of this Triangulation. Note, that each manifold
* @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
*/
std::vector<types::manifold_id> get_manifold_ids() const;
-
+
/**
* Copy a triangulation. This operation is not cheap, so you should
* be careful with using this. We do not implement this function as
* a pointer.
*/
std::map<unsigned int, types::boundary_id> *vertex_to_boundary_id_map_1d;
-
+
/**
* A map that relates the number of a boundary vertex to the
ar &check_for_distorted_cells;
- if (dim == 1) {
- ar &vertex_to_boundary_id_map_1d;
- ar &vertex_to_manifold_id_map_1d;
- }
+ if (dim == 1)
+ {
+ ar &vertex_to_boundary_id_map_1d;
+ ar &vertex_to_manifold_id_map_1d;
+ }
}
"same setting with regard to reporting distorted "
"cell as the one previously stored."));
- if (dim == 1) {
- ar &vertex_to_boundary_id_map_1d;
- ar &vertex_to_manifold_id_map_1d;
- }
+ if (dim == 1)
+ {
+ ar &vertex_to_boundary_id_map_1d;
+ ar &vertex_to_manifold_id_map_1d;
+ }
// trigger the create signal to indicate
// that new content has been imported into
// -------------------------------------------------------------------
// Explicit invalid things...
-template <>
-const Manifold<2,1> & Triangulation<2, 1>::get_manifold(const types::manifold_id) const;
-template <>
-const Manifold<3,1> & Triangulation<3, 1>::get_manifold(const types::manifold_id) const;
-template <>
-const Manifold<3,2> & Triangulation<3, 2>::get_manifold(const types::manifold_id) const;
+template <>
+const Manifold<2,1> &Triangulation<2, 1>::get_manifold(const types::manifold_id) const;
+template <>
+const Manifold<3,1> &Triangulation<3, 1>::get_manifold(const types::manifold_id) const;
+template <>
+const Manifold<3,2> &Triangulation<3, 2>::get_manifold(const types::manifold_id) const;
#endif // DOXYGEN
* get_manifold() internally.
*/
const Boundary<dim,spacedim> &get_boundary () const;
-
+
/**
* Return a constant reference to the manifold object used for this object.
*
/**
* @}
- */
-
+ */
+
/**
* @name User data
/**
* Return the manifold indicator of this
- * object.
+ * object.
*
* @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
*/
*
* @see @ref GlossManifoldIndicator "Glossary entry on manifold indicators"
*/
-
+
void
set_all_manifold_ids (const types::manifold_id);
/**
template <int structdim, int dim, int spacedim>
const Boundary<dim,spacedim> &
TriaAccessor<structdim, dim, spacedim>::get_boundary () const
-{
+{
Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
-
- // Get the default (manifold_id)
+
+ // Get the default (manifold_id)
const types::manifold_id mi = this->objects().manifold_id[this->present_index];
-
- // In case this is not valid, check
- // the boundary id, after having
- // casted it to a manifold id
- if(mi == numbers::invalid_manifold_id)
+
+ // In case this is not valid, check
+ // the boundary id, after having
+ // casted it to a manifold id
+ if (mi == numbers::invalid_manifold_id)
return this->tria->get_boundary(structdim < dim ?
- this->objects().boundary_or_material_id[this->present_index].boundary_id:
- dim < spacedim ?
- this->objects().boundary_or_material_id[this->present_index].material_id:
- numbers::invalid_manifold_id);
+ this->objects().boundary_or_material_id[this->present_index].boundary_id:
+ dim < spacedim ?
+ this->objects().boundary_or_material_id[this->present_index].material_id:
+ numbers::invalid_manifold_id);
else
return this->tria->get_boundary(mi);
}
template <int structdim, int dim, int spacedim>
const Manifold<dim,spacedim> &
TriaAccessor<structdim, dim, spacedim>::get_manifold () const
-{
+{
Assert (this->used(), TriaAccessorExceptions::ExcCellNotUsed());
-
- // Get the default (manifold_id)
+
+ // Get the default (manifold_id)
const types::manifold_id mi = this->objects().manifold_id[this->present_index];
-
- // In case this is not valid, check
- // the boundary id, after having
- // casted it to a manifold id
- if(mi == numbers::invalid_manifold_id)
+
+ // In case this is not valid, check
+ // the boundary id, after having
+ // casted it to a manifold id
+ if (mi == numbers::invalid_manifold_id)
return this->tria->get_manifold(structdim < dim ?
- this->objects().boundary_or_material_id[this->present_index].boundary_id:
- dim < spacedim ?
- this->objects().boundary_or_material_id[this->present_index].material_id:
- numbers::invalid_manifold_id);
+ this->objects().boundary_or_material_id[this->present_index].boundary_id:
+ dim < spacedim ?
+ this->objects().boundary_or_material_id[this->present_index].material_id:
+ numbers::invalid_manifold_id);
else
return this->tria->get_manifold(mi);
}
switch (structdim)
{
case 1:
- if(dim == 1)
- {
- (*this->tria->vertex_to_manifold_id_map_1d)
- [vertex_index(0)] = manifold_ind;
- (*this->tria->vertex_to_manifold_id_map_1d)
- [vertex_index(1)] = manifold_ind;
- }
+ if (dim == 1)
+ {
+ (*this->tria->vertex_to_manifold_id_map_1d)
+ [vertex_index(0)] = manifold_ind;
+ (*this->tria->vertex_to_manifold_id_map_1d)
+ [vertex_index(1)] = manifold_ind;
+ }
break;
case 2:
types::manifold_id
TriaAccessor<0, 1, spacedim>::manifold_id () const
{
- if( tria->vertex_to_manifold_id_map_1d->find (this->vertex_index())
- != tria->vertex_to_manifold_id_map_1d->end())
+ if ( tria->vertex_to_manifold_id_map_1d->find (this->vertex_index())
+ != tria->vertex_to_manifold_id_map_1d->end())
return (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()];
else
return numbers::invalid_manifold_id;
* Store manifold ids. This field
* stores the manifold id of each object, which
* is a number between 0 and
- * numbers::invalid_manifold_id-1.
+ * numbers::invalid_manifold_id-1.
*/
std::vector<types::manifold_id> manifold_id;
DEAL_II_NAMESPACE_OPEN
-namespace Manifolds {
+namespace Manifolds
+{
- Quadrature<3>
- get_default_quadrature(const TriaIterator<CellAccessor<3, 3> >& obj)
+ Quadrature<3>
+ get_default_quadrature(const TriaIterator<CellAccessor<3, 3> > &obj)
{
- std::vector<Point<3> > sp;
- std::vector<double> wp;
-
- const int dim = 3;
-
- const unsigned int np =
- GeometryInfo<dim>::vertices_per_cell+
- GeometryInfo<dim>::lines_per_cell+
- GeometryInfo<dim>::faces_per_cell;
-
- sp.resize(np);
- wp.resize(np);
- unsigned int j=0;
-
- // note that the exact weights are chosen such as to minimize the
- // distortion of the eight new hexes from the optimal shape; their
- // derivation and values is copied over from the
- // @p{MappingQ::set_laplace_on_vector} function
- for(unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i, ++j)
- {
- sp[j] = obj->vertex(i);
- wp[j] = 1.0/128.0;
- }
- for(unsigned int i=0; i<GeometryInfo<dim>::lines_per_cell; ++i, ++j)
- {
- sp[j] = (obj->line(i)->has_children() ? obj->line(i)->child(0)->vertex(1) :
- obj->line(i)->get_manifold().get_new_point_on_line(obj->line(i)));
- wp[j] = 7.0/192.0;
- }
- for(unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i, ++j)
- {
- sp[j] = (obj->face(i)->has_children() ? obj->face(i)->isotropic_child(0)->vertex(3) :
- obj->face(i)->get_manifold().get_new_point_on_face(obj->face(i)));
- wp[j] = 1.0/12.0;
- }
- return Quadrature<3>(sp,wp);
+
std::vector<Point<3> > sp;
+ std::vector<double> wp;
+
+ const int dim = 3;
+
+ const unsigned int np =
+ GeometryInfo<dim>::vertices_per_cell+
+ GeometryInfo<dim>::lines_per_cell+
+ GeometryInfo<dim>::faces_per_cell;
+
+ sp.resize(np);
+ wp.resize(np);
+ unsigned int j=0;
+
+ // note that the exact weights are chosen such as to minimize the
+ // distortion of the eight new hexes from the optimal shape; their
+ // derivation and values is copied over from the
+ // @p{MappingQ::set_laplace_on_vector} function
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i, ++j)
+ {
+ sp[j] = obj->vertex(i);
+ wp[j] = 1.0/128.0;
+ }
+ for (unsigned int i=0; i<GeometryInfo<dim>::lines_per_cell; ++i, ++j)
+ {
+ sp[j] = (obj->line(i)->has_children() ? obj->line(i)->child(0)->vertex(1) :
+ obj->line(i)->get_manifold().get_new_point_on_line(obj->line(i)));
+ wp[j] = 7.0/192.0;
+ }
+ for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i, ++j)
+ {
+ sp[j] = (obj->face(i)->has_children() ? obj->face(i)->isotropic_child(0)->vertex(3) :
+ obj->face(i)->get_manifold().get_new_point_on_face(obj->face(i)));
+ wp[j] = 1.0/12.0;
+ }
+ return Quadrature<3>(sp,wp);
}
}
Point<spacedim>
Manifold<dim, spacedim>::
project_to_manifold (const std::vector<Point<spacedim> > &,
- const Point<spacedim> &) const
+ const Point<spacedim> &) const
{
Assert (false, ExcPureFunctionCalled());
return Point<spacedim>();
#ifdef DEBUG
double sum=0;
- for(unsigned int i=0; i<weights.size(); ++i)
+ for (unsigned int i=0; i<weights.size(); ++i)
sum+= weights[i];
Assert(std::abs(sum-1.0) < 1e-10, ExcMessage("Weights should sum to 1!"));
#endif
-
- for(unsigned int i=0; i<surrounding_points.size(); ++i)
+
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
p += surrounding_points[i]*weights[i];
return project_to_manifold(surrounding_points, p);
get_new_point_on_line (const typename Triangulation<dim, spacedim>::line_iterator &line) const
{
return get_new_point
- (get_default_quadrature<const typename Triangulation<dim, spacedim>::line_iterator,
- spacedim>(line, false));
+ (get_default_quadrature<const typename Triangulation<dim, spacedim>::line_iterator,
+ spacedim>(line, false));
}
get_new_point_on_quad (const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
{
return get_new_point
- (get_default_quadrature<const typename Triangulation<dim, spacedim>::quad_iterator,
- spacedim>(quad, false));
+ (get_default_quadrature<const typename Triangulation<dim, spacedim>::quad_iterator,
+ spacedim>(quad, false));
}
template <>
Point<3>
Manifold<3,3>::
-get_new_point_on_hex (const typename Triangulation<3, 3>::hex_iterator &hex) const{
+get_new_point_on_hex (const typename Triangulation<3, 3>::hex_iterator &hex) const
+{
return get_new_point(get_default_quadrature(hex));
}
template <int dim, int spacedim>
-FlatManifold<dim,spacedim>::FlatManifold (const Point<spacedim> periodicity,
- const double tolerance) :
+FlatManifold<dim,spacedim>::FlatManifold (const Point<spacedim> periodicity,
+ const double tolerance) :
tolerance(tolerance),
periodicity(periodicity)
{}
#ifdef DEBUG
double sum=0;
- for(unsigned int i=0; i<weights.size(); ++i)
+ for (unsigned int i=0; i<weights.size(); ++i)
sum+= weights[i];
Assert(std::abs(sum-1.0) < tolerance, ExcMessage("Weights should sum to 1!"));
#endif
-
-
+
+
Point<spacedim> p;
Point<spacedim> dp;
Point<spacedim> minP = periodicity;
const bool check_period = (periodicity.norm() > tolerance);
- if(check_period)
- for(unsigned int i=0; i<surrounding_points.size(); ++i)
- for(unsigned int d=0; d<spacedim; ++d) {
- minP[d] = std::min(minP[d], surrounding_points[i][d]);
- if(periodicity[d] > 0)
- Assert( (surrounding_points[i][d] < periodicity[d]+tolerance) ||
- (surrounding_points[i][d] >= -tolerance),
- ExcPeriodicBox(d, surrounding_points[i], periodicity, tolerance));
- }
-
- for(unsigned int i=0; i<surrounding_points.size(); ++i) {
- dp = Point<spacedim>();
- if(check_period) {
- for(unsigned int d=0; d<spacedim; ++d)
- if(periodicity[d] > 0)
- dp[d] = ( (surrounding_points[i][d]-minP[d]) > periodicity[d]/2.0 ?
- -periodicity[d] : 0.0 );
+ if (check_period)
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
+ for (unsigned int d=0; d<spacedim; ++d)
+ {
+ minP[d] = std::min(minP[d], surrounding_points[i][d]);
+ if (periodicity[d] > 0)
+ Assert( (surrounding_points[i][d] < periodicity[d]+tolerance) ||
+ (surrounding_points[i][d] >= -tolerance),
+ ExcPeriodicBox(d, surrounding_points[i], periodicity, tolerance));
+ }
+
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
+ {
+ dp = Point<spacedim>();
+ if (check_period)
+ {
+ for (unsigned int d=0; d<spacedim; ++d)
+ if (periodicity[d] > 0)
+ dp[d] = ( (surrounding_points[i][d]-minP[d]) > periodicity[d]/2.0 ?
+ -periodicity[d] : 0.0 );
+ }
+ p += (surrounding_points[i]+dp)*weights[i];
}
- p += (surrounding_points[i]+dp)*weights[i];
- }
- if(check_period)
- for(unsigned int d=0; d<spacedim; ++d)
- if(periodicity[d] > 0)
- p[d] = (p[d] < 0 ? p[d] + periodicity[d] : p[d]);
+ if (check_period)
+ for (unsigned int d=0; d<spacedim; ++d)
+ if (periodicity[d] > 0)
+ p[d] = (p[d] < 0 ? p[d] + periodicity[d] : p[d]);
return project_to_manifold(surrounding_points, p);
}
template <int dim, int spacedim>
-Point<spacedim>
-FlatManifold<dim, spacedim>::project_to_manifold (const std::vector<Point<spacedim> > & vertices,
- const Point<spacedim> &candidate) const
+Point<spacedim>
+FlatManifold<dim, spacedim>::project_to_manifold (const std::vector<Point<spacedim> > &vertices,
+ const Point<spacedim> &candidate) const
{
return candidate;
}
const std::vector<Point<spacedim> > &surrounding_points = quad.get_points();
const std::vector<double> &weights = quad.get_weights();
std::vector<Point<chartdim> > chart_points(surrounding_points.size());
-
- for(unsigned int i=0; i<surrounding_points.size(); ++i)
+
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
chart_points[i] = pull_back(surrounding_points[i]);
const Quadrature<chartdim> chart_quad(chart_points, weights);
const Point<chartdim> p_chart = sub_manifold.get_new_point(chart_quad);
-
+
return push_forward(p_chart);
}
template <int dim, int spacedim>
Point<spacedim>
-SphericalManifold<dim,spacedim>::get_periodicity() {
+SphericalManifold<dim,spacedim>::get_periodicity()
+{
Point<spacedim> periodicity;
periodicity[spacedim-1] = 2*numbers::PI; // theta and phi period.
return periodicity;
template <int dim, int spacedim>
Point<spacedim>
-SphericalManifold<dim,spacedim>::get_new_point(const Quadrature<spacedim> &quad) const {
- if(spacedim == 2)
+SphericalManifold<dim,spacedim>::get_new_point(const Quadrature<spacedim> &quad) const
+{
+ if (spacedim == 2)
return ManifoldChart<dim,spacedim,spacedim>::get_new_point(quad);
- else {
- double rho_average = 0;
- Point<spacedim> mid_point;
- for(unsigned int i=0; i<quad.size(); ++i) {
- rho_average += quad.weight(i)*(quad.point(i)-center).norm();
- mid_point += quad.weight(i)*quad.point(i);
+ else
+ {
+ double rho_average = 0;
+ Point<spacedim> mid_point;
+ for (unsigned int i=0; i<quad.size(); ++i)
+ {
+ rho_average += quad.weight(i)*(quad.point(i)-center).norm();
+ mid_point += quad.weight(i)*quad.point(i);
+ }
+ // Project the mid_pont back to the right location
+ Point<spacedim> R = mid_point-center;
+ // Scale it to have radius rho_average
+ R *= rho_average/R.norm();
+ // And return it.
+ return center+R;
}
- // Project the mid_pont back to the right location
- Point<spacedim> R = mid_point-center;
- // Scale it to have radius rho_average
- R *= rho_average/R.norm();
- // And return it.
- return center+R;
- }
}
template <int dim, int spacedim>
Point<spacedim>
-SphericalManifold<dim,spacedim>::push_forward(const Point<spacedim> &spherical_point) const {
- Assert(spherical_point[0] >=0.0,
- ExcMessage("Negative radius for given point."));
+SphericalManifold<dim,spacedim>::push_forward(const Point<spacedim> &spherical_point) const
+{
+ Assert(spherical_point[0] >=0.0,
+ ExcMessage("Negative radius for given point."));
const double rho = spherical_point[0];
const double theta = spherical_point[1];
-
+
Point<spacedim> p;
- if(rho > 1e-10)
- switch(spacedim) {
- case 2:
- p[0] = rho*cos(theta);
- p[1] = rho*sin(theta);
- break;
- case 3:
+ if (rho > 1e-10)
+ switch (spacedim)
+ {
+ case 2:
+ p[0] = rho*cos(theta);
+ p[1] = rho*sin(theta);
+ break;
+ case 3:
{
- const double &phi= spherical_point[2];
- p[0] = rho*sin(theta)*cos(phi);
- p[1] = rho*sin(theta)*sin(phi);
- p[2] = rho*cos(theta);
+ const double &phi= spherical_point[2];
+ p[0] = rho*sin(theta)*cos(phi);
+ p[1] = rho*sin(theta)*sin(phi);
+ p[2] = rho*cos(theta);
}
break;
- default:
- Assert(false, ExcInternalError());
- }
+ default:
+ Assert(false, ExcInternalError());
+ }
return p+center;
}
template <int dim, int spacedim>
Point<spacedim>
-SphericalManifold<dim,spacedim>::pull_back(const Point<spacedim> &space_point) const {
+SphericalManifold<dim,spacedim>::pull_back(const Point<spacedim> &space_point) const
+{
const Point<spacedim> R = space_point-center;
const double rho = R.norm();
const double x = R[0];
const double y = R[1];
-
+
Point<spacedim> p;
p[0] = rho;
-
- switch(spacedim) {
- case 2:
- p[1] = atan2(y,x);
- if(p[1] < 0)
- p[1] += 2*numbers::PI;
- break;
- case 3:
+
+ switch (spacedim)
+ {
+ case 2:
+ p[1] = atan2(y,x);
+ if (p[1] < 0)
+ p[1] += 2*numbers::PI;
+ break;
+ case 3:
{
const double z = R[2];
p[2] = atan2(y,x); // phi
- if(p[2] < 0)
- p[2] += 2*numbers::PI; // phi is periodic
+ if (p[2] < 0)
+ p[2] += 2*numbers::PI; // phi is periodic
p[1] = atan2(sqrt(x*x+y*y),z); // theta
}
break;
- default:
- Assert(false, ExcInternalError());
- }
+ default:
+ Assert(false, ExcInternalError());
+ }
return p;
}
// ============================================================
template <int dim, int spacedim>
-CylindricalManifold<dim,spacedim>::CylindricalManifold(const unsigned int axis,
- const double tolerance) :
+CylindricalManifold<dim,spacedim>::CylindricalManifold(const unsigned int axis,
+ const double tolerance) :
direction (Point<spacedim>::unit_vector(axis)),
point_on_axis (Point<spacedim>()),
tolerance(tolerance)
template <int dim, int spacedim>
CylindricalManifold<dim,spacedim>::CylindricalManifold(const Point<spacedim> &direction,
- const Point<spacedim> &point_on_axis,
- const double tolerance) :
+ const Point<spacedim> &point_on_axis,
+ const double tolerance) :
direction (direction),
point_on_axis (point_on_axis),
tolerance(tolerance)
const std::vector<Point<spacedim> > &surrounding_points = quad.get_points();
const std::vector<double> &weights = quad.get_weights();
- // compute a proposed new point
+ // compute a proposed new point
Point<spacedim> middle = flat_manifold.get_new_point(quad);
double radius = 0;
Point<spacedim> on_plane;
-
- for(unsigned int i=0; i<surrounding_points.size(); ++i)
+
+ for (unsigned int i=0; i<surrounding_points.size(); ++i)
{
on_plane = surrounding_points[i]-point_on_axis;
on_plane = on_plane - (on_plane*direction) * direction;
radius += weights[i]*on_plane.norm();
}
-
+
// we then have to project this point out to the given radius from
// the axis. to this end, we have to take into account the offset
// point_on_axis and the direction of the axis
const Point<spacedim> vector_from_axis = (middle-point_on_axis) -
- ((middle-point_on_axis) * direction) * direction;
+ ((middle-point_on_axis) * direction) * direction;
// scale it to the desired length and put everything back together,
// unless we have a point on the axis
else
return (vector_from_axis / vector_from_axis.norm() * radius +
- ((middle-point_on_axis) * direction) * direction +
- point_on_axis);
+ ((middle-point_on_axis) * direction) * direction +
+ point_on_axis);
}
-
+
// explicit instantiations
#include "manifold_lib.inst"
cells[cell].vertices[1]));
next_free_line->set_used_flag ();
next_free_line->set_material_id (cells[cell].material_id);
- next_free_line->set_manifold_id (cells[cell].manifold_id);
+ next_free_line->set_manifold_id (cells[cell].manifold_id);
next_free_line->clear_user_data ();
next_free_line->set_subdomain_id (0);
// finally set the
// vertex_to_boundary_id_map_1d
- // and vertex_to_manifold_id_map_1d
+ // and vertex_to_manifold_id_map_1d
// maps
triangulation.vertex_to_boundary_id_map_1d->clear();
- triangulation.vertex_to_manifold_id_map_1d->clear();
+ triangulation.vertex_to_manifold_id_map_1d->clear();
for (typename Triangulation<dim,spacedim>::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end(); ++cell)
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
- {
- (*triangulation
- .vertex_to_manifold_id_map_1d)[cell->face(f)->vertex_index()]
- = numbers::invalid_manifold_id;
-
- if (cell->at_boundary(f))
- (*triangulation
- .vertex_to_boundary_id_map_1d)[cell->face(f)->vertex_index()]
- = f;
- }
+ {
+ (*triangulation
+ .vertex_to_manifold_id_map_1d)[cell->face(f)->vertex_index()]
+ = numbers::invalid_manifold_id;
+
+ if (cell->at_boundary(f))
+ (*triangulation
+ .vertex_to_boundary_id_map_1d)[cell->face(f)->vertex_index()]
+ = f;
+ }
}
cell->set_used_flag ();
cell->set_material_id (cells[c].material_id);
- cell->set_manifold_id (cells[c].manifold_id);
+ cell->set_manifold_id (cells[c].manifold_id);
cell->clear_user_data ();
cell->set_subdomain_id (0);
else
// interior line -> numbers::internal_face_boundary_id
line->set_boundary_indicator (numbers::internal_face_boundary_id);
- line->set_manifold_id(numbers::flat_manifold_id);
+ line->set_manifold_id(numbers::flat_manifold_id);
}
// set boundary indicators where
ExcInteriorLineCantBeBoundary());
line->set_boundary_indicator (boundary_line->boundary_id);
- line->set_manifold_id (boundary_line->manifold_id);
+ line->set_manifold_id (boundary_line->manifold_id);
}
cell->set_used_flag ();
cell->set_material_id (cells[c].material_id);
- cell->set_manifold_id (cells[c].manifold_id);
+ cell->set_manifold_id (cells[c].manifold_id);
cell->clear_user_flag ();
cell->clear_user_data ();
cell->set_subdomain_id (0);
else
// interior quad -> numbers::internal_face_boundary_id
quad->set_boundary_indicator (numbers::internal_face_boundary_id);
- // Manifold ids are set
- // independently of where
- // they are
- quad->set_manifold_id(numbers::flat_manifold_id);
+ // Manifold ids are set
+ // independently of where
+ // they are
+ quad->set_manifold_id(numbers::flat_manifold_id);
}
/////////////////////////////////////////
// for this: first mark all lines as interior. use this loop
// to also set all manifold ids of all lines
for (typename Triangulation<dim,spacedim>::line_iterator
- line=triangulation.begin_line(); line!=triangulation.end_line(); ++line)
- {
- line->set_boundary_indicator (numbers::internal_face_boundary_id);
- line->set_manifold_id(numbers::flat_manifold_id);
- }
-
+ line=triangulation.begin_line(); line!=triangulation.end_line(); ++line)
+ {
+ line->set_boundary_indicator (numbers::internal_face_boundary_id);
+ line->set_manifold_id(numbers::flat_manifold_id);
+ }
+
// next reset all lines bounding
// boundary quads as on the
// boundary also. note that since
"if they carry the same boundary indicator."));
line->set_boundary_indicator (boundary_line->boundary_id);
- // Set manifold id if given
- line->set_manifold_id(boundary_line->manifold_id);
+ // Set manifold id if given
+ line->set_manifold_id(boundary_line->manifold_id);
}
"if they carry the same boundary indicator."));
quad->set_boundary_indicator (boundary_quad->boundary_id);
- quad->set_manifold_id (boundary_quad->manifold_id);
+ quad->set_manifold_id (boundary_quad->manifold_id);
}
switch_1->set_line_orientation(2, switch_2->line_orientation(2));
switch_1->set_line_orientation(3, switch_2->line_orientation(3));
switch_1->set_boundary_indicator(switch_2->boundary_indicator());
- switch_1->set_manifold_id(switch_2->manifold_id());
+ switch_1->set_manifold_id(switch_2->manifold_id());
switch_1->set_user_index(switch_2->user_index());
if (switch_2->user_flag_set())
switch_1->set_user_flag();
switch_2->set_line_orientation(2, switch_1_line_orientations[2]);
switch_2->set_line_orientation(3, switch_1_line_orientations[3]);
switch_2->set_boundary_indicator(switch_1_boundary_indicator);
- switch_2->set_manifold_id(switch_1->manifold_id());
+ switch_2->set_manifold_id(switch_1->manifold_id());
switch_2->set_user_index(switch_1_user_index);
if (switch_1_user_flag)
switch_2->set_user_flag();
if (dim == spacedim)
{
// triangulation.vertices[next_unused_vertex] = new_point;
- triangulation.vertices[next_unused_vertex] =
- cell->get_manifold().get_new_point
- (Manifolds::get_default_quadrature
- <typename Triangulation<dim,spacedim>::active_cell_iterator,
- spacedim>(cell));
+ triangulation.vertices[next_unused_vertex] =
+ cell->get_manifold().get_new_point
+ (Manifolds::get_default_quadrature
+ <typename Triangulation<dim,spacedim>::active_cell_iterator,
+ spacedim>(cell));
// if the user_flag is set, i.e. if the
// cell is at the boundary, use a
// points on this face and on the opposite face,
// as returned by the underlying manifold
// object.
- {
- std::vector<Point<spacedim> > ps(2);
- std::vector<double> ws(2, 0.5);
- ps[0] = cell->face(boundary_face)
- ->child(0)->vertex(1);
- ps[1] = cell->face(GeometryInfo<dim>
- ::opposite_face[boundary_face])
- ->child(0)->vertex(1);
- Quadrature<spacedim> qs(ps,ws);
- triangulation.vertices[next_unused_vertex]
- = cell->get_manifold().get_new_point(qs);
- }
- }
+ {
+
std::vector<Point<spacedim> > ps(2);
+ std::vector<double> ws(2, 0.5);
+ ps[0] = cell->face(boundary_face)
+ ->child(0)->vertex(1);
+ ps[1] = cell->face(GeometryInfo<dim>
+ ::opposite_face[boundary_face])
+ ->child(0)->vertex(1);
+ Quadrature<spacedim> qs(ps,ws);
+ triangulation.vertices[next_unused_vertex]
+ = cell->get_manifold().get_new_point(qs);
+ }
+ }
}
else
{
// new vertex is placed on the surface according to
// the information stored in the boundary class
- const Manifold<dim,spacedim> &manifold = cell->get_manifold();
-
+ const Manifold<dim,spacedim> &manifold = cell->get_manifold();
+
triangulation.vertices[next_unused_vertex] =
manifold.get_new_point_on_quad (cell);
}
new_lines[l]->clear_children();
// interior line
new_lines[l]->set_boundary_indicator(numbers::internal_face_boundary_id);
- new_lines[l]->set_manifold_id(cell->manifold_id());
+ new_lines[l]->set_manifold_id(cell->manifold_id());
}
// Now add the four (two)
// inherit material
// properties
subcells[i]->set_material_id (cell->material_id());
- subcells[i]->set_manifold_id (cell->manifold_id());
+ subcells[i]->set_manifold_id (cell->manifold_id());
subcells[i]->set_subdomain_id (subdomainid);
if (i%2==0)
Assert (next_unused_vertex < triangulation.vertices.size(),
ExcTooFewVerticesAllocated());
- // Now we always ask the manifold where to put the
- // new point. The get_manifold function will return
- // a flat boundary if invalid_manifold_id is set,
- // behaving as before in the flat case. Backward
- // compatibility requires us to use the material id
- // of the cell in the codimension one case, however
- // we only do this if the manifold_id is the invalid
- // one, otherwise use the material id. This is done
- // internally in the get_manifold() function.
- triangulation.vertices[next_unused_vertex] =
- cell->get_manifold().get_new_point_on_line(cell);
- triangulation.vertices_used[next_unused_vertex] = true;
-
+ // Now we always ask the manifold where to put the
+ // new point. The get_manifold function will return
+ // a flat boundary if invalid_manifold_id is set,
+ // behaving as before in the flat case. Backward
+ // compatibility requires us to use the material id
+ // of the cell in the codimension one case, however
+ // we only do this if the manifold_id is the invalid
+ // one, otherwise use the material id. This is done
+ // internally in the get_manifold() function.
+ triangulation.vertices[next_unused_vertex] =
+ cell->get_manifold().get_new_point_on_line(cell);
+ triangulation.vertices_used[next_unused_vertex] = true;
+
// search for next two unused cell (++ takes care of
// the end of the vector)
typename Triangulation<dim,spacedim>::raw_cell_iterator
::TriaObject<1> (cell->vertex_index(0),
next_unused_vertex));
first_child->set_material_id (cell->material_id());
- first_child->set_manifold_id (cell->manifold_id());
+ first_child->set_manifold_id (cell->manifold_id());
first_child->set_subdomain_id (subdomainid);
first_child->set_direction_flag (cell->direction_flag());
-
+
first_child->set_parent (cell->index ());
- // Set manifold id of the right face. Only do this
- // on the first child.
- first_child->face(1)->set_manifold_id(cell->manifold_id());
-
+ // Set manifold id of the right face. Only do this
+ // on the first child.
+ first_child->face(1)->set_manifold_id(cell->manifold_id());
+
// reset neighborship info (refer to
// internal::Triangulation::TriaLevel<0> for
// details)
cell->vertex_index(1)));
second_child->set_neighbor (0, first_child);
second_child->set_material_id (cell->material_id());
- second_child->set_manifold_id (cell->manifold_id());
+ second_child->set_manifold_id (cell->manifold_id());
second_child->set_subdomain_id (subdomainid);
second_child->set_direction_flag (cell->direction_flag());
// boundary lines differently; for interior
// lines we can compute the midpoint as the mean
// of the two vertices: if (line->at_boundary())
- triangulation.vertices[next_unused_vertex]
- = line->get_manifold().get_new_point_on_line (line);
+ triangulation.vertices[next_unused_vertex]
+ = line->get_manifold().get_new_point_on_line (line);
}
else
// however, if spacedim>dim, we always have to ask
// same object of the cell (which was stored in
// line->user_index() before) unless a manifold_id
// has been set on this very line.
- if(line->manifold_id() == numbers::invalid_manifold_id)
- triangulation.vertices[next_unused_vertex]
- = triangulation.get_manifold(line->user_index()).get_new_point_on_line (line);
- else
- triangulation.vertices[next_unused_vertex]
- = line->get_manifold().get_new_point_on_line (line);
-
+ if (line->manifold_id() == numbers::invalid_manifold_id)
+ triangulation.vertices[next_unused_vertex]
+ = triangulation.get_manifold(line->user_index()).get_new_point_on_line (line);
+ else
+ triangulation.vertices[next_unused_vertex]
+ = line->get_manifold().get_new_point_on_line (line);
+
// now that we created the right point, make up the
// two child lines. To this end, find a pair of
// unused lines
ExcTooFewVerticesAllocated());
triangulation.vertices_used[next_unused_vertex] = true;
- triangulation.vertices[next_unused_vertex]
- = line->get_manifold().get_new_point_on_line (line);
+ triangulation.vertices[next_unused_vertex]
+ = line->get_manifold().get_new_point_on_line (line);
// now that we created the right point, make up the
// two child lines (++ takes care of the end of the
children[0]->set_boundary_indicator (line->boundary_indicator());
children[1]->set_boundary_indicator (line->boundary_indicator());
-
- children[0]->set_manifold_id (line->manifold_id());
+
+ children[0]->set_manifold_id (line->manifold_id());
children[1]->set_manifold_id (line->manifold_id());
// finally clear flag
new_line->clear_user_data();
new_line->clear_children();
new_line->set_boundary_indicator(quad->boundary_indicator());
- new_line->set_manifold_id(quad->manifold_id());
+ new_line->set_manifold_id(quad->manifold_id());
// child 0 and 1 of a line are switched if the
// line orientation is false. set up a miniature
new_quads[i]->clear_user_data();
new_quads[i]->clear_children();
new_quads[i]->set_boundary_indicator (quad->boundary_indicator());
- new_quads[i]->set_manifold_id (quad->manifold_id());
+ new_quads[i]->set_manifold_id (quad->manifold_id());
// set all line orientations to true, change
// this after the loop, as we have to consider
// different lines for each child
new_child[i]->set(internal::Triangulation::TriaObject<1>(old_child[i]->vertex_index(0),
old_child[i]->vertex_index(1)));
new_child[i]->set_boundary_indicator(old_child[i]->boundary_indicator());
- new_child[i]->set_manifold_id(old_child[i]->manifold_id());
+ new_child[i]->set_manifold_id(old_child[i]->manifold_id());
new_child[i]->set_user_index(old_child[i]->user_index());
if (old_child[i]->user_flag_set())
new_child[i]->set_user_flag();
switch_1->set_line_orientation(2, switch_2->line_orientation(2));
switch_1->set_line_orientation(3, switch_2->line_orientation(3));
switch_1->set_boundary_indicator(switch_2->boundary_indicator());
- switch_1->set_manifold_id(switch_2->manifold_id());
+ switch_1->set_manifold_id(switch_2->manifold_id());
switch_1->set_user_index(switch_2->user_index());
if (switch_2->user_flag_set())
switch_1->set_user_flag();
switch_2->set_line_orientation(2, switch_1_line_orientations[2]);
switch_2->set_line_orientation(3, switch_1_line_orientations[3]);
switch_2->set_boundary_indicator(switch_1_boundary_indicator);
- switch_2->set_manifold_id(switch_1->manifold_id());
+ switch_2->set_manifold_id(switch_1->manifold_id());
switch_2->set_user_index(switch_1_user_index);
if (switch_1_user_flag)
switch_2->set_user_flag();
// set the middle vertex
// appropriately
- if (quad->at_boundary() ||
- (quad->manifold_id() != numbers::invalid_manifold_id) )
+ if (quad->at_boundary() ||
+ (quad->manifold_id() != numbers::invalid_manifold_id) )
triangulation.vertices[next_unused_vertex]
- = quad->get_manifold().get_new_point_on_quad (quad);
+ = quad->get_manifold().get_new_point_on_quad (quad);
else
- {
- // it might be that the quad itself is not at
- // the boundary, but that one of its lines
- // actually is. in this case, the newly
- // created vertices at the centers of the
- // lines are not necessarily the mean values
- // of the adjacent vertices, so do not compute
- // the new vertex as the mean value of the 4
- // vertices of the face, but rather as a
- // weighted mean value of the 8 vertices which
- // we already have (the four old ones, and the
- // four ones inserted as middle points for the
- // four lines). summing up some more points is
- // generally cheaper than first asking whether
- // one of the lines is at the boundary
- //
- // note that the exact weights are chosen such
- // as to minimize the distortion of the four
- // new quads from the optimal shape; their
- // derivation and values is copied over from
- // the @p{MappingQ::set_laplace_on_vector}
- // function
- triangulation.vertices[next_unused_vertex] =
- quad->get_manifold().get_new_point
- (Manifolds::get_default_quadrature
- <typename Triangulation<dim,spacedim>::quad_iterator,
- spacedim>(quad, true));
- }
- triangulation.vertices_used[next_unused_vertex] = true;
+ {
+ // it might be that the quad itself is not at
+ // the boundary, but that one of its lines
+ // actually is. in this case, the newly
+ // created vertices at the centers of the
+ // lines are not necessarily the mean values
+ // of the adjacent vertices, so do not compute
+ // the new vertex as the mean value of the 4
+ // vertices of the face, but rather as a
+ // weighted mean value of the 8 vertices which
+ // we already have (the four old ones, and the
+ // four ones inserted as middle points for the
+ // four lines). summing up some more points is
+ // generally cheaper than first asking whether
+ // one of the lines is at the boundary
+ //
+ // note that the exact weights are chosen such
+ // as to minimize the distortion of the four
+ // new quads from the optimal shape; their
+ // derivation and values is copied over from
+ // the @p{MappingQ::set_laplace_on_vector}
+ // function
+ triangulation.vertices[next_unused_vertex] =
+ quad->get_manifold().get_new_point
+ (Manifolds::get_default_quadrature
+ <typename Triangulation<dim,spacedim>::quad_iterator,
+ spacedim>(quad, true));
+ }
+ triangulation.vertices_used[next_unused_vertex] = true;
// now that we created the right point, make up
// the four lines interior to the quad (++ takes
// care of the end of the vector)
// *--2--*
//
// the lines are numbered as follows:
- //
+ //
// *--*--*
// | 1 |
// *2-*-3*
new_lines[i]->clear_user_data();
new_lines[i]->clear_children();
new_lines[i]->set_boundary_indicator(quad->boundary_indicator());
- new_lines[i]->set_manifold_id(quad->manifold_id());
+ new_lines[i]->set_manifold_id(quad->manifold_id());
}
// now for the quads. again, first collect some
// data about the indices of the lines, with the
// following numbering:
- //
+ //
// .-6-.-7-.
// 1 9 3
// .-10.11-.
new_quads[i]->clear_user_data();
new_quads[i]->clear_children();
new_quads[i]->set_boundary_indicator (quad->boundary_indicator());
- new_quads[i]->set_manifold_id (quad->manifold_id());
+ new_quads[i]->set_manifold_id (quad->manifold_id());
// set all line orientations to true, change
// this after the loop, as we have to consider
// different lines for each child
// interior quad
new_quads[i]->set_boundary_indicator (numbers::internal_face_boundary_id);
// they inherit geometry description of the hex they belong to
- new_quads[i]->set_manifold_id (hex->manifold_id());
+ new_quads[i]->set_manifold_id (hex->manifold_id());
// set all line orientation flags to true by
// default, change this afterwards, if necessary
for (unsigned int j=0; j<GeometryInfo<dim>::lines_per_face; ++j)
// inherit material
// properties
new_hexes[i]->set_material_id (hex->material_id());
- new_hexes[i]->set_manifold_id (hex->manifold_id());
+ new_hexes[i]->set_manifold_id (hex->manifold_id());
new_hexes[i]->set_subdomain_id (subdomainid);
if (i%2)
// the quads are numbered as follows:
//
// planes in the interior of the old hex:
- //
+ //
// *
// /|
// / | x
// *
//
// children of the faces of the old hex
- //
+ //
// *---*---* *---*---*
// /| | | / / /|
// / | | | / 9 / 10/ |
// the quads are numbered as follows:
//
// planes in the interior of the old hex:
- //
+ //
// *
// /|
// / | x
// *
//
// children of the faces of the old hex
- //
+ //
// *-------* *-------*
// /| | / 10 /|
// * | | *-------* |
// the quads are numbered as follows:
//
// planes in the interior of the old hex:
- //
+ //
// *
// /|
// / | x
// *
//
// children of the faces of the old hex
- //
+ //
// *---*---* *-------*
// /| 8 | / /|
// / | | / 10 / |
// the quads are numbered as follows:
//
// planes in the interior of the old hex:
- //
+ //
// *
// /|
// * | x
// *
//
// children of the faces of the old hex
- //
+ //
// *---*---* *---*---*
// /| | | /18 / 19/|
// * |10 | 11| /---/---* |
// the quads are numbered as follows:
//
// planes in the interior of the old hex:
- //
+ //
// *
// /|
// / | x
// the quads are numbered as follows:
//
// planes in the interior of the old hex:
- //
+ //
// *
// /|
// / | x
// function, and like the new vertex at the
// center of the quad is weighted (see above)
- // triangulation.vertices[next_unused_vertex] = Point<dim>();
- // // first add corners of hex
- // for (unsigned int vertex=0;
- // vertex<GeometryInfo<dim>::vertices_per_cell; ++vertex)
- // triangulation.vertices[next_unused_vertex] += hex->vertex(vertex) / 128;
- // // now add center of lines
- // for (unsigned int line=0;
- // line<GeometryInfo<dim>::lines_per_cell; ++line)
- // triangulation.vertices[next_unused_vertex] += hex->line(line)->child(0)->vertex(1) *
- // 7./192.;
- // // finally add centers of faces. note that
- // // vertex 3 of child 0 is an invariant with
- // // respect to the face orientation, flip and
- // // rotation
- // for (unsigned int face=0;
- // face<GeometryInfo<dim>::faces_per_cell; ++face)
- // triangulation.vertices[next_unused_vertex] += hex->face(face)->isotropic_child(0)->vertex(3) *
- // 1./12.;
-
- triangulation.vertices[next_unused_vertex] =
- hex->get_manifold().get_new_point_on_hex(hex);
+ triangulation.vertices[next_unused_vertex] =
+ hex->get_manifold().get_new_point_on_hex(hex);
// set the data of the six lines. first collect
// the indices of the seven vertices (consider
// the two planes to be crossed to form the
// planes cutting the hex in two vertically and
// horizontally)
- //
+ //
// *--3--* *--5--*
// / / / | | |
// 0--6--1 0--6--1
// the quads are numbered as follows:
//
// planes in the interior of the old hex:
- //
+ //
// *
// /|
// * |
if (dim == 1)
{
vertex_to_boundary_id_map_1d
- = new std::map<unsigned int, types::boundary_id>();
+ = new std::map<unsigned int, types::boundary_id>();
vertex_to_manifold_id_map_1d
- = new std::map<unsigned int, types::manifold_id>();
+ = new std::map<unsigned int, types::manifold_id>();
}
// connect the any_change signal to the other signals
const Manifold<dim, spacedim> &manifold_object)
{
Assert(m_number < numbers::invalid_manifold_id,
- ExcIndexRange(m_number,0,numbers::invalid_manifold_id));
+ ExcIndexRange(m_number,0,numbers::invalid_manifold_id));
manifold[m_number] = &manifold_object;
}
Triangulation<dim, spacedim>::set_manifold (const types::manifold_id m_number)
{
Assert(m_number < numbers::invalid_manifold_id,
- ExcIndexRange(m_number,0,numbers::invalid_manifold_id));
+ ExcIndexRange(m_number,0,numbers::invalid_manifold_id));
//delete the entry located at number.
manifold.erase(m_number);
const Boundary<dim,spacedim> &
Triangulation<dim, spacedim>::get_boundary (const types::manifold_id m_number) const
{
- const Boundary<dim, spacedim> * man =
+ const Boundary<dim, spacedim> *man =
dynamic_cast<const Boundary<dim, spacedim> *>(&get_manifold(m_number));
Assert(man != NULL,
- ExcMessage("You tried to get a Boundary, but I only have a Manifold."));
-
+ ExcMessage("You tried to get a Boundary, but I only have a Manifold."));
+
return *man;
}
const Manifold<dim,spacedim> &
Triangulation<dim, spacedim>::get_manifold (const types::manifold_id m_number) const
{
- //look, if there is a manifold stored at
- //manifold_id number.
+ //look, if there is a manifold stored at
+ //manifold_id number.
typename std::map<types::manifold_id, SmartPointer<const Manifold<dim,spacedim>, Triangulation<dim, spacedim> > >::const_iterator it
= manifold.find(m_number);
{
std::set<types::manifold_id> m_ids;
active_cell_iterator cell=begin_active();
- for (; cell!=end(); ++cell) {
- m_ids.insert(cell->manifold_id());
- if(dim>1)
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->at_boundary(face))
- m_ids.insert(cell->face(face)->manifold_id());
- }
+ for (; cell!=end(); ++cell)
+ {
+ m_ids.insert(cell->manifold_id());
+ if (dim>1)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->at_boundary(face))
+ m_ids.insert(cell->face(face)->manifold_id());
+ }
std::vector<types::manifold_id> manifold_indicators(m_ids.begin(), m_ids.end());
return manifold_indicators;
}
faces = new internal::Triangulation::TriaFaces<dim>(*old_tria.faces);
typename std::map<types::manifold_id,
- SmartPointer<const Manifold<dim,spacedim> , Triangulation<dim, spacedim> > >::const_iterator
- bdry_iterator = old_tria.manifold.begin();
+ SmartPointer<const Manifold<dim,spacedim> , Triangulation<dim, spacedim> > >::const_iterator
+ bdry_iterator = old_tria.manifold.begin();
for (; bdry_iterator != old_tria.manifold.end() ; bdry_iterator++)
manifold[bdry_iterator->first] = bdry_iterator->second;
template <int dim, int spacedim>
bool Triangulation<dim, spacedim>::has_hanging_nodes () const
{
- for (unsigned int lvl = 0; lvl<n_global_levels()-1;lvl++)
+ for (unsigned int lvl = 0; lvl<n_global_levels()-1; lvl++)
if (n_active_cells(lvl) != 0)
return true;
// coarsen flags on level 1 to avoid level 0 cells being created
// by coarsening. As coarsen flags will never be added, this can
// be done once and for all before the actual loop starts.
- //
+ //
// 1/ do not coarsen a cell if 'most of the neighbors' will be
// refined after the step. This is to prevent occurrence of
// unrefined islands.
- //
+ //
// 2/ eliminate refined islands in the interior and at the
// boundary. since they don't do much harm besides increasing the
// number of degrees of freedom, doing this has a rather low
// priority.
- //
+ //
// 3/ limit the level difference of neighboring cells at each
// vertex.
- //
+ //
// 4/ eliminate unrefined islands. this has higher priority since
// this diminishes the approximation properties not only of the
// unrefined island, but also of the surrounding patch.
- //
+ //
// 5/ ensure patch level 1. Then the triangulation consists of
// patches, i.e. of cells that are refined once. It follows that
// if at least one of the children of a cell is or will be
// eliminate_unrefined_islands, eliminate_refined_inner_islands
// and eliminate_refined_boundary_islands will be fulfilled
// automatically and do not need to be enforced separately.
- //
+ //
// 6/ take care of the requirement that no double refinement is done
// at each face
- //
+ //
// 7/ take care that no double refinement is done at each line in 3d
// or higher dimensions.
- //
+ //
// 8/ make sure that all children of each cell are either flagged
// for coarsening or none of the children is
//
}
//////////////////////////////////////
- // STEP 3:
+ // STEP 3:
// limit the level difference of neighboring cells at each
// vertex.
//
}
template <>
-const Manifold<2,1> & Triangulation<2, 1>::get_manifold(const types::manifold_id) const {
+const Manifold<2,1> &Triangulation<2, 1>::get_manifold(const types::manifold_id) const
+{
Assert(false, ExcImpossibleInDim(1));
static FlatManifold<2,1> flat;
return flat;
}
template <>
-const Manifold<3,1> & Triangulation<3, 1>::get_manifold(const types::manifold_id) const {
+const Manifold<3,1> &Triangulation<3, 1>::get_manifold(const types::manifold_id) const
+{
Assert(false, ExcImpossibleInDim(1));
static FlatManifold<3,1> flat;
return flat;
}
template <>
-const Manifold<3,2> & Triangulation<3, 2>::get_manifold(const types::manifold_id) const {
+const Manifold<3,2> &Triangulation<3, 2>::get_manifold(const types::manifold_id) const
+{
Assert(false, ExcImpossibleInDim(2));
static FlatManifold<3,2> flat;
return flat;
for (unsigned int c=0; c<this->n_children(); ++c)
this->child(c)->set_all_manifold_ids (manifold_ind);
- // for hexes also set manifold_id
- // of bounding quads and lines
-
- // Six bonding quads
- for(unsigned int i=0; i<6; ++i)
+ // for hexes also set manifold_id
+ // of bounding quads and lines
+
+ // Six bonding quads
+ for (unsigned int i=0; i<6; ++i)
this->quad(i)->set_manifold_id(manifold_ind);
- // Twelve bounding lines
+ // Twelve bounding lines
for (unsigned int i=0; i<12; ++i)
this->line(i)->set_manifold_id (manifold_ind);
}
manifold_id.reserve (new_size);
manifold_id.insert (manifold_id.end(),
- new_size-manifold_id.size(),
- numbers::flat_manifold_id);
+ new_size-manifold_id.size(),
+ numbers::flat_manifold_id);
}
manifold_id.reserve (new_size);
manifold_id.insert (manifold_id.end(),
- new_size-manifold_id.size(),
- numbers::flat_manifold_id);
+ new_size-manifold_id.size(),
+ numbers::flat_manifold_id);
user_data.reserve (new_size);
user_data.resize (new_size);
MemoryConsumption::memory_consumption (used) +
MemoryConsumption::memory_consumption (user_flags) +
MemoryConsumption::memory_consumption (boundary_or_material_id) +
- MemoryConsumption::memory_consumption (manifold_id) +
+ MemoryConsumption::memory_consumption (manifold_id) +
MemoryConsumption::memory_consumption (refinement_cases) +
user_data.capacity() * sizeof(UserData) + sizeof(user_data));
}
projects / dealii.git / commitdiff