[cig-commits] r1348 - in trunk/aspect: include/aspect/particle source/postprocess

[heien at dealii.org]

Author: heien
Date: 2012-11-06 17:44:54 -0700 (Tue, 06 Nov 2012)
New Revision: 1348

Modified:
   trunk/aspect/include/aspect/particle/integrator.h
   trunk/aspect/source/postprocess/tracer.cc
Log:
More work on hybrid integration scheme (currently only for research purposes)


Modified: trunk/aspect/include/aspect/particle/integrator.h
===================================================================
--- trunk/aspect/include/aspect/particle/integrator.h	2012-11-06 20:29:13 UTC (rev 1347)
+++ trunk/aspect/include/aspect/particle/integrator.h	2012-11-07 00:44:54 UTC (rev 1348)
@@ -392,14 +392,16 @@
         std::map<double, IntegrationScheme>        _scheme;
         const parallel::distributed::Triangulation<dim>   *_tria;
         const DoFHandler<dim>           *_dh;
+        const Mapping<dim>              *_mapping;
+        const TrilinosWrappers::MPI::BlockVector *_solution;
 
-        virtual IntegrationScheme select_scheme(const std::vector<Point<dim> > &cell_vertices, const std::vector<Point<dim> > &cell_velocities)
+        virtual IntegrationScheme select_scheme(const std::vector<Point<dim> > &cell_vertices, const std::vector<Point<dim> > &cell_velocities, const double timestep)
         {
           return cell_vertices[0][0] > 0.5 ? SCHEME_RK4 : SCHEME_EULER;
         };
 
       public:
-        HybridIntegrator(const parallel::distributed::Triangulation<dim> *new_tria, const DoFHandler<dim> *new_dh)
+        HybridIntegrator(const parallel::distributed::Triangulation<dim> *new_tria, const DoFHandler<dim> *new_dh, const Mapping<dim> *new_mapping, const TrilinosWrappers::MPI::BlockVector *new_solution)
         {
           _step = 0;
           _loc0.clear();
@@ -409,6 +411,8 @@
           _scheme.clear();
           _tria = new_tria;
           _dh = new_dh;
+          _mapping = new_mapping;
+          _solution = new_solution;
         };
 
         virtual bool integrate_step(std::multimap<LevelInd, T> &particles, double dt)
@@ -418,17 +422,23 @@
           double                              id_num;
           LevelInd                            cur_level_ind;
           IntegrationScheme                   cur_scheme;
-          typename parallel::distributed::Triangulation<dim>::cell_iterator found_cell;
+          //typename parallel::distributed::Triangulation<dim>::cell_iterator found_cell;
+          typename DoFHandler<dim>::active_cell_iterator found_cell;
+          Functions::FEFieldFunction<dim, DoFHandler<dim>, TrilinosWrappers::MPI::BlockVector> fe_value(*_dh, *_solution, *_mapping);
 
           // If this is the first step, go through all the cells and determine
           // which integration scheme the particles in each cell should use
           if (_step == 0)
             {
+              Vector<double>                single_res(dim+2);
               std::vector<Point<dim> >    cell_vertices, cell_velocities;
+              std::vector<Vector<double> > temp_vals;
+
               cell_vertices.resize(GeometryInfo<dim>::vertices_per_cell);
               cell_velocities.resize(GeometryInfo<dim>::vertices_per_cell);
               cur_level_ind.first = cur_level_ind.second = -1;
               cur_scheme = SCHEME_UNDEFINED;
+              temp_vals.resize(GeometryInfo<dim>::vertices_per_cell, single_res);
 
               for (it=particles.begin(); it!=particles.end(); ++it)
                 {
@@ -436,9 +446,21 @@
                     {
                       cur_level_ind = it->first;
                       found_cell = typename DoFHandler<dim>::active_cell_iterator(_tria, cur_level_ind.first, cur_level_ind.second, _dh);
+
+                      // Ideally we should use all quadrature point velocities, but for now
+                      // we just evaluate them at the vertices
+                      fe_value.set_active_cell(found_cell);
                       for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i) cell_vertices[i] = found_cell->vertex(i);
+                      fe_value.vector_value_list(cell_vertices, temp_vals);
+                      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+                        {
+                          for (unsigned int d=0; d<dim; ++d)
+                            {
+                              cell_velocities[i][d] = temp_vals[i][d];
+                            }
+                        }
 
-                      cur_scheme = select_scheme(cell_vertices, cell_velocities);
+                      cur_scheme = select_scheme(cell_vertices, cell_velocities, dt);
                     }
                   _scheme[it->second.id_num()] = cur_scheme;
                 }

Modified: trunk/aspect/source/postprocess/tracer.cc
===================================================================
--- trunk/aspect/source/postprocess/tracer.cc	2012-11-06 20:29:13 UTC (rev 1347)
+++ trunk/aspect/source/postprocess/tracer.cc	2012-11-07 00:44:54 UTC (rev 1348)
@@ -70,7 +70,10 @@
             }
           else if (_integration_scheme == "hybrid")
             {
-              _integrator = new Particle::HybridIntegrator<dim, Particle::BaseParticle<dim> >(&(this->get_triangulation()), &(this->get_dof_handler()));
+              _integrator = new Particle::HybridIntegrator<dim, Particle::BaseParticle<dim> >(&(this->get_triangulation()),
+                                                                                              &(this->get_dof_handler()),
+                                                                                              &(this->get_mapping()),
+                                                                                              &(this->get_solution()));
             }
 
           // Set up the particle world with the appropriate simulation objects