[cig-commits] commit 2453 by bangerth to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Wed Apr 9 14:58:33 PDT 2014
Revision 2453
Minor edits.
U trunk/aspect/include/aspect/simulator.h
U trunk/aspect/source/simulator/nullspace.cc
U trunk/aspect/source/simulator/parameters.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2453&peg=2453
Diff:
Modified: trunk/aspect/include/aspect/simulator.h
===================================================================
--- trunk/aspect/include/aspect/simulator.h 2014-04-09 21:29:45 UTC (rev 2452)
+++ trunk/aspect/include/aspect/simulator.h 2014-04-09 21:58:30 UTC (rev 2453)
@@ -919,13 +919,19 @@
/**
- * interpolates the given function onto the velocity FE space and write it into the given vector.
+ * Interpolate the given function onto the velocity FE space and write it into the given vector.
+ *
+ * This function is implemented in
+ * <code>source/simulator/helper_functions.cc</code>.
*/
void interpolate_onto_velocity_system(const TensorFunction<1,dim> &func,
LinearAlgebra::Vector &vec);
/**
- * Sets up data structures for null space removal. Called after every mesh refinement.
+ * Set up data structures for null space removal. Called after every mesh refinement.
+ *
+ * This function is implemented in
+ * <code>source/simulator/nullspace.cc</code>.
*/
void setup_nullspace_removal();
@@ -935,6 +941,9 @@
*
* @param relevant_dst locally relevant vector for the whole FE, will be filled at the end.
* @param tmp_distributed_stokes only contains velocity and pressure.
+ *
+ * This function is implemented in
+ * <code>source/simulator/nullspace.cc</code>.
*/
void remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
LinearAlgebra::BlockVector &tmp_distributed_stokes);
Modified: trunk/aspect/source/simulator/nullspace.cc
===================================================================
--- trunk/aspect/source/simulator/nullspace.cc 2014-04-09 21:29:45 UTC (rev 2452)
+++ trunk/aspect/source/simulator/nullspace.cc 2014-04-09 21:58:30 UTC (rev 2453)
@@ -46,24 +46,54 @@
+ /**
+ * A class we use when setting up the data structures for nullspace removal
+ * of the rotations in spherical or annular shells.
+ */
template<int dim>
class Rotation : public TensorFunction<1,dim>
- {
- private:
- Tensor<1,dim> axis;
- public:
- Rotation(const unsigned int a) : axis(Tensor<1,dim>(Point<dim>::unit_vector(a))) {}
- virtual Tensor<1,dim> value (const Point<dim> &p) const { Tensor<1,dim> vel; if( dim == 2) cross_product(vel, p); else cross_product(vel, axis, p); return vel;}
- };
+ {
+ private:
+ Tensor<1,dim> axis;
+ public:
+ Rotation(const unsigned int a)
+ :
+ axis(Tensor<1,dim>(Point<dim>::unit_vector(a)))
+ {}
+
+ virtual Tensor<1,dim> value (const Point<dim> &p) const
+ {
+ Tensor<1,dim> vel;
+ if( dim == 2)
+ cross_product(vel, p);
+ else
+ cross_product(vel, axis, p);
+ return vel;
+ }
+ };
+
+
+ /**
+ * A class we use when setting up the data structures for nullspace removal
+ * of the translations in box-like geometries.
+ */
template <int dim>
class Translation : public TensorFunction<1,dim>
{
- private:
- const unsigned int direction;
- public:
- Translation(const unsigned int d) : direction(d) {}
- virtual Tensor<1,dim> value(const Point<dim> &) const { return Tensor<1,dim>(Point<dim>::unit_vector(direction)); }
+ private:
+ const unsigned int direction;
+
+ public:
+ Translation(const unsigned int d)
+ :
+ direction(d)
+ {}
+
+ virtual Tensor<1,dim> value(const Point<dim> &) const
+ {
+ return Point<dim>::unit_vector(direction);
+ }
};
}
@@ -71,37 +101,46 @@
template <int dim>
void Simulator<dim>::setup_nullspace_removal()
{
- if (parameters.nullspace_removal & NullspaceRemoval::translational_momentum)
- AssertThrow(false, ExcNotImplemented());
-
std::vector<std_cxx1x::shared_ptr<TensorFunction<1,dim> > > funcs;
- if (parameters.nullspace_removal & NullspaceRemoval::net_rotation)
+ switch (parameters.nullspace_removal)
+ {
+ case NullspaceRemoval::net_rotation:
{
if (dim==2)
funcs.push_back(std_cxx1x::shared_ptr<TensorFunction<1,dim> >(new internal::Rotation<dim>(0)));
if (dim==3)
for(unsigned int a=0; a<dim; ++a)
funcs.push_back(std_cxx1x::shared_ptr<TensorFunction<1,dim> >(new internal::Rotation<dim>(a)));
+
+ break;
}
- if (parameters.nullspace_removal & NullspaceRemoval::net_translation)
+ case NullspaceRemoval::net_translation:
{
for(unsigned int a=0; a<dim; ++a)
funcs.push_back(std_cxx1x::shared_ptr<TensorFunction<1,dim> >(new internal::Translation<dim>(a)));
+ break;
}
+ case NullspaceRemoval::translational_momentum:
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ break;
+ }
+
+
if (funcs.size()>0)
{
net_rotations_translations.resize(funcs.size());
for (unsigned int i=0;i<funcs.size();++i)
- net_rotations_translations[i].reinit(
- introspection.index_sets.system_partitioning[introspection.block_indices.velocities],
- mpi_communicator);
-
- unsigned int idx = 0;
- for (unsigned int i=0;i<funcs.size();++i)
{
+ // for each of the null space dimensions, set up
+ // a vector, fill it with an element of the null space,
+ // and normalize it
+ net_rotations_translations[i].reinit(
+ introspection.index_sets.system_partitioning[introspection.block_indices.velocities],
+ mpi_communicator);
interpolate_onto_velocity_system(*funcs[i],
net_rotations_translations[i]);
net_rotations_translations[i] /= net_rotations_translations[i].l2_norm();
@@ -110,40 +149,60 @@
}
}
+
+
template <int dim>
void Simulator<dim>::remove_nullspace(LinearAlgebra::BlockVector &relevant_dst, LinearAlgebra::BlockVector &tmp_distributed_stokes)
{
- if (parameters.nullspace_removal & NullspaceRemoval::net_rotation ||
- parameters.nullspace_removal & NullspaceRemoval::net_translation)
+ switch (parameters.nullspace_removal)
+ {
+ case NullspaceRemoval::net_rotation:
+ case NullspaceRemoval::net_translation:
{
for(unsigned int i=0; i<net_rotations_translations.size(); ++i)
{
- double power = net_rotations_translations[i]
+ // compute the magnitude of the solution vector in direction
+ // of this null space vector and subtract the corresponding multiple
+ const double power = net_rotations_translations[i]
* tmp_distributed_stokes.block(introspection.block_indices.velocities);
tmp_distributed_stokes.block(introspection.block_indices.velocities).sadd(1.0,
-1.0*power,
net_rotations_translations[i]);
}
relevant_dst.block(0) = tmp_distributed_stokes.block(0);
+ break;
}
- if (parameters.nullspace_removal & NullspaceRemoval::angular_momentum)
+
+ case NullspaceRemoval::angular_momentum:
{
remove_net_angular_momentum( relevant_dst, tmp_distributed_stokes);
+ break;
}
+
+ case NullspaceRemoval::none:
+ break;
+
+ default:
+ AssertThrow (false, ExcInternalError());
+ }
}
+
+
template <>
void
- Simulator<3>::remove_net_angular_momentum( LinearAlgebra::BlockVector &relevant_dst, LinearAlgebra::BlockVector &tmp_distributed_stokes )
+ Simulator<3>::remove_net_angular_momentum( LinearAlgebra::BlockVector &relevant_dst,
+ LinearAlgebra::BlockVector &tmp_distributed_stokes )
{
AssertThrow(false, ExcNotImplemented());
}
+
template <>
void
- Simulator<2>::remove_net_angular_momentum( LinearAlgebra::BlockVector &relevant_dst, LinearAlgebra::BlockVector &tmp_distributed_stokes )
+ Simulator<2>::remove_net_angular_momentum( LinearAlgebra::BlockVector &relevant_dst,
+ LinearAlgebra::BlockVector &tmp_distributed_stokes )
{
-
// compute and remove angular momentum from velocity field, by computing
// int rho V
More information about the CIG-COMMITS
mailing list