Merge pull request #35 from nicola-giuliani/spheroid

added the possibility of using a spheroid

Author: nicola-giuliani

include/bem_problem.h

@@ -239,6 +239,7 @@ class BEMProblem : public ParameterAcceptor
   void adaptive_refinement(const TrilinosWrappers::MPI::Vector &error_vector);
 
 
+
   ConditionalOStream pcout;
   ComputationalDomain<dim> &comp_dom;
 

include/computational_domain.h

@@ -249,6 +249,10 @@ class ComputationalDomain : public ParameterAcceptor
   std::map<unsigned int, std::vector<typename Triangulation<dim-1,dim>::active_cell_iterator> > vert_to_elems;
   std::set<typename Triangulation<dim-1,dim>::active_cell_iterator> edge_cells;
   Manifold<dim-1, dim> *manifold;
+
+  bool spheroid_bool, used_spherical_manifold;
+  double spheroid_x_axis, spheroid_y_axis, spheroid_z_axis;
+
   ConditionalOStream pcout;
 
 

source/bem_problem.cc

@@ -148,6 +148,16 @@ void BEMProblem<dim>::reinit()
   // mapping_degree = fe->get_degree();
   if (!mapping)
     {
+      if (comp_dom.spheroid_bool && comp_dom.used_spherical_manifold)
+        {
+          for (types::global_dof_index ii=0; ii<gradient_dh.n_dofs()/dim; ++ii)
+            {
+              map_vector[vec_original_to_sub_wise[ii]] *= comp_dom.spheroid_x_axis;
+              map_vector[vec_original_to_sub_wise[ii+gradient_dh.n_dofs()/dim]] *= comp_dom.spheroid_y_axis;
+              if (dim == 3)
+                map_vector[vec_original_to_sub_wise[ii+gradient_dh.n_dofs()/dim]] *= comp_dom.spheroid_z_axis;
+            }
+        }
       if (mapping_type == "FE")
         mapping = SP(new MappingFEField<dim-1, dim> (gradient_dh, map_vector));
       else

source/computational_domain.cc

@@ -39,6 +39,7 @@ ComputationalDomain<dim>::ComputationalDomain(MPI_Comm comm)
   mpi_communicator (comm),
   n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)),
   this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)),
+  used_spherical_manifold(false),
   pcout(std::cout)
 
 {
@@ -115,6 +116,18 @@ void ComputationalDomain<dim>::declare_parameters (ParameterHandler &prm)
   }
   prm.leave_subsection();
 
+  prm.declare_entry("Use a spheroid", "false",
+                    Patterns::Bool());
+
+  prm.declare_entry("Axis x dimension", "2.",
+                    Patterns::Double());
+
+  prm.declare_entry("Axis y dimension", "3.",
+                    Patterns::Double());
+
+  prm.declare_entry("Axis z dimension", "4.",
+                    Patterns::Double());
+
 
 }
 
@@ -134,6 +147,12 @@ void ComputationalDomain<dim>::parse_parameters (ParameterHandler &prm)
   max_curvature_ref_cycles = prm.get_integer("Maximum number of curvature adaptive refinement cycles");
   cad_to_projectors_tolerance_ratio = prm.get_double("Cad tolerance to projectors tolerance ratio");
 
+
+  spheroid_bool = prm.get_bool("Use a spheroid");
+  spheroid_x_axis = prm.get_double("Axis x dimension");
+  spheroid_y_axis = prm.get_double("Axis y dimension");
+  spheroid_z_axis = prm.get_double("Axis z dimension");
+
   prm.enter_subsection("Boundary Conditions ID Numbers");
   {
     std::vector<std::string> dirichlet_string_list = Utilities::split_string_list(prm.get("Dirichlet boundary ids"));
@@ -266,6 +285,8 @@ void ComputationalDomain<dim>::read_domain()
       tria.set_all_manifold_ids(0);
       tria.set_manifold(0, *manifold);
 
+      used_spherical_manifold = true;
+
     }
 
 }