uniform_refinement_estimator.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // C++ includes
19 #include <algorithm> // for std::fill
20 #include <sstream>
21 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for sqrt
23 
24 
25 // Local Includes
26 #include "libmesh/dof_map.h"
27 #include "libmesh/elem.h"
29 #include "libmesh/error_vector.h"
30 #include "libmesh/fe.h"
31 #include "libmesh/libmesh_common.h"
33 #include "libmesh/mesh_base.h"
35 #include "libmesh/numeric_vector.h"
36 #include "libmesh/quadrature.h"
37 #include "libmesh/system.h"
39 #include "libmesh/partitioner.h"
40 #include "libmesh/tensor_tools.h"
41 
42 #ifdef LIBMESH_ENABLE_AMR
43 
44 namespace libMesh
45 {
46 
47 //-----------------------------------------------------------------
48 // ErrorEstimator implementations
50  ErrorVector& error_per_cell,
51  const NumericVector<Number>* solution_vector,
52  bool estimate_parent_error)
53 {
54  START_LOG("estimate_error()", "UniformRefinementEstimator");
55  std::map<const System*, const NumericVector<Number>* > solution_vectors;
56  solution_vectors[&_system] = solution_vector;
57  this->_estimate_error (NULL, &_system, &error_per_cell, NULL, NULL,
58  &solution_vectors, estimate_parent_error);
59  STOP_LOG("estimate_error()", "UniformRefinementEstimator");
60 }
61 
63  ErrorVector& error_per_cell,
64  const std::map<const System*, SystemNorm>& error_norms,
65  const std::map<const System*, const NumericVector<Number>* >* solution_vectors,
66  bool estimate_parent_error)
67 {
68  START_LOG("estimate_errors()", "UniformRefinementEstimator");
69  this->_estimate_error (&_es, NULL, &error_per_cell, NULL,
70  &error_norms, solution_vectors,
71  estimate_parent_error);
72  STOP_LOG("estimate_errors()", "UniformRefinementEstimator");
73 }
74 
76  ErrorMap& errors_per_cell,
77  const std::map<const System*, const NumericVector<Number>* >* solution_vectors,
78  bool estimate_parent_error)
79 {
80  START_LOG("estimate_errors()", "UniformRefinementEstimator");
81  this->_estimate_error (&_es, NULL, NULL, &errors_per_cell, NULL,
82  solution_vectors, estimate_parent_error);
83  STOP_LOG("estimate_errors()", "UniformRefinementEstimator");
84 }
85 
87  const System* _system,
88  ErrorVector* error_per_cell,
89  ErrorMap* errors_per_cell,
90  const std::map<const System*, SystemNorm > *_error_norms,
91  const std::map<const System*, const NumericVector<Number>* >* solution_vectors,
92  bool)
93 {
94  // Get a vector of the Systems we're going to work on,
95  // and set up a error_norms map if necessary
96  std::vector<System *> system_list;
99  (new std::map<const System*, SystemNorm>);
100 
101  if (_es)
102  {
103  libmesh_assert(!_system);
104  libmesh_assert(_es->n_systems());
105  _system = &(_es->get_system(0));
106  libmesh_assert_equal_to (&(_system->get_equation_systems()), _es);
107 
108  libmesh_assert(_es->n_systems());
109  for (unsigned int i=0; i != _es->n_systems(); ++i)
110  // We have to break the rules here, because we can't refine a const System
111  system_list.push_back(const_cast<System *>(&(_es->get_system(i))));
112 
113  // If we're computing one vector, we need to know how to scale
114  // each variable's contributions to it.
115  if (_error_norms)
116  {
117  libmesh_assert(!errors_per_cell);
118  }
119  else
120  // If we're computing many vectors, we just need to know which
121  // variables to skip
122  {
123  libmesh_assert (errors_per_cell);
124 
125  _error_norms = error_norms.get();
126 
127  for (unsigned int i=0; i!= _es->n_systems(); ++i)
128  {
129  const System &sys = _es->get_system(i);
130  unsigned int n_vars = sys.n_vars();
131 
132  std::vector<Real> weights(n_vars, 0.0);
133  for (unsigned int v = 0; v != n_vars; ++v)
134  {
135  if (errors_per_cell->find(std::make_pair(&sys, v)) ==
136  errors_per_cell->end())
137  continue;
138 
139  weights[v] = 1.0;
140  }
141  (*error_norms)[&sys] =
142  SystemNorm(std::vector<FEMNormType>(n_vars, error_norm.type(0)),
143  weights);
144  }
145  }
146  }
147  else
148  {
149  libmesh_assert(_system);
150  // We have to break the rules here, because we can't refine a const System
151  system_list.push_back(const_cast<System *>(_system));
152 
153  libmesh_assert(!_error_norms);
154  (*error_norms)[_system] = error_norm;
155  _error_norms = error_norms.get();
156  }
157 
158  // An EquationSystems reference will be convenient.
159  // We have to break the rules here, because we can't refine a const System
160  EquationSystems& es =
161  const_cast<EquationSystems &>(_system->get_equation_systems());
162 
163  // The current mesh
164  MeshBase& mesh = es.get_mesh();
165 
166  // The dimensionality of the mesh
167  const unsigned int dim = mesh.mesh_dimension();
168 
169  // Resize the error_per_cell vectors to be
170  // the number of elements, initialize them to 0.
171  if (error_per_cell)
172  {
173  error_per_cell->clear();
174  error_per_cell->resize (mesh.max_elem_id(), 0.);
175  }
176  else
177  {
178  libmesh_assert(errors_per_cell);
179  for (ErrorMap::iterator i = errors_per_cell->begin();
180  i != errors_per_cell->end(); ++i)
181  {
182  ErrorVector *e = i->second;
183  e->clear();
184  e->resize(mesh.max_elem_id(), 0.);
185  }
186  }
187 
188  // We'll want to back up all coarse grid vectors
189  std::vector<std::map<std::string, NumericVector<Number> *> >
190  coarse_vectors(system_list.size());
191  std::vector<NumericVector<Number> *>
192  coarse_solutions(system_list.size());
193  std::vector<NumericVector<Number> *>
194  coarse_local_solutions(system_list.size());
195  // And make copies of projected solutions
196  std::vector<NumericVector<Number> *>
197  projected_solutions(system_list.size());
198 
199  // And we'll need to temporarily change solution projection settings
200  std::vector<bool> old_projection_settings(system_list.size());
201 
202  // And it'll be best to avoid any repartitioning
203  AutoPtr<Partitioner> old_partitioner = mesh.partitioner();
204  mesh.partitioner().reset(NULL);
205 
206  for (unsigned int i=0; i != system_list.size(); ++i)
207  {
208  System &system = *system_list[i];
209 
210  // Check for valid error_norms
211  libmesh_assert (_error_norms->find(&system) !=
212  _error_norms->end());
213 
214  // Back up the solution vector
215  coarse_solutions[i] = system.solution->clone().release();
216  coarse_local_solutions[i] =
217  system.current_local_solution->clone().release();
218 
219  // Back up all other coarse grid vectors
220  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
221  system.vectors_end(); ++vec)
222  {
223  // The (string) name of this vector
224  const std::string& var_name = vec->first;
225 
226  coarse_vectors[i][var_name] = vec->second->clone().release();
227  }
228 
229  // Use a non-standard solution vector if necessary
230  if (solution_vectors &&
231  solution_vectors->find(&system) != solution_vectors->end() &&
232  solution_vectors->find(&system)->second &&
233  solution_vectors->find(&system)->second != system.solution.get())
234  {
235  NumericVector<Number>* newsol =
236  const_cast<NumericVector<Number>*>
237  (solution_vectors->find(&system)->second);
238  newsol->swap(*system.solution);
239  system.update();
240  }
241 
242  // Make sure the solution is projected when we refine the mesh
243  old_projection_settings[i] = system.project_solution_on_reinit();
244  system.project_solution_on_reinit() = true;
245  }
246 
247  // Find the number of coarse mesh elements, to make it possible
248  // to find correct coarse elem ids later
249  const dof_id_type max_coarse_elem_id = mesh.max_elem_id();
250 #ifndef NDEBUG
251  // n_coarse_elem is only used in an assertion later so
252  // avoid declaring it unless asserts are active.
253  const dof_id_type n_coarse_elem = mesh.n_elem();
254 #endif
255 
256  // Uniformly refine the mesh
257  MeshRefinement mesh_refinement(mesh);
258 
260 
261  // FIXME: this may break if there is more than one System
262  // on this mesh but estimate_error was still called instead of
263  // estimate_errors
264  for (unsigned int i = 0; i != number_h_refinements; ++i)
265  {
266  mesh_refinement.uniformly_refine(1);
267  es.reinit();
268  }
269 
270  for (unsigned int i = 0; i != number_p_refinements; ++i)
271  {
272  mesh_refinement.uniformly_p_refine(1);
273  es.reinit();
274  }
275 
276  for (unsigned int i=0; i != system_list.size(); ++i)
277  {
278  System &system = *system_list[i];
279 
280  // Copy the projected coarse grid solutions, which will be
281  // overwritten by solve()
282 // projected_solutions[i] = system.solution->clone().release();
283  projected_solutions[i] = NumericVector<Number>::build(system.comm()).release();
284  projected_solutions[i]->init(system.solution->size(), true, SERIAL);
285  system.solution->localize(*projected_solutions[i],
286  system.get_dof_map().get_send_list());
287  }
288 
289  // Are we doing a forward or an adjoint solve?
290  bool solve_adjoint = false;
291  if (solution_vectors)
292  {
293  System *sys = system_list[0];
294  libmesh_assert (solution_vectors->find(sys) !=
295  solution_vectors->end());
296  const NumericVector<Number> *vec = solution_vectors->find(sys)->second;
297  for (unsigned int j=0; j != sys->qoi.size(); ++j)
298  {
299  std::ostringstream adjoint_name;
300  adjoint_name << "adjoint_solution" << j;
301 
302  if (vec == sys->request_vector(adjoint_name.str()))
303  {
304  solve_adjoint = true;
305  break;
306  }
307  }
308  }
309 
310  // Get the uniformly refined solution.
311 
312  if (_es)
313  {
314  // Even if we had a decent preconditioner, valid matrix etc. before
315  // refinement, we don't any more.
316  for (unsigned int i=0; i != es.n_systems(); ++i)
317  es.get_system(i).disable_cache();
318 
319  // No specified vectors == forward solve
320  if (!solution_vectors)
321  es.solve();
322  else
323  {
324  libmesh_assert_equal_to (solution_vectors->size(), es.n_systems());
325  libmesh_assert (solution_vectors->find(system_list[0]) !=
326  solution_vectors->end());
327  libmesh_assert(solve_adjoint ||
328  (solution_vectors->find(system_list[0])->second ==
329  system_list[0]->solution.get()) ||
330  !solution_vectors->find(system_list[0])->second);
331 
332 #ifdef DEBUG
333  for (unsigned int i=0; i != system_list.size(); ++i)
334  {
335  System *sys = system_list[i];
336  libmesh_assert (solution_vectors->find(sys) !=
337  solution_vectors->end());
338  const NumericVector<Number> *vec = solution_vectors->find(sys)->second;
339  if (solve_adjoint)
340  {
341  bool found_vec = false;
342  for (unsigned int j=0; j != sys->qoi.size(); ++j)
343  {
344  std::ostringstream adjoint_name;
345  adjoint_name << "adjoint_solution" << j;
346 
347  if (vec == sys->request_vector(adjoint_name.str()))
348  {
349  found_vec = true;
350  break;
351  }
352  }
353  libmesh_assert(found_vec);
354  }
355  else
356  libmesh_assert(vec == sys->solution.get() || !vec);
357  }
358 #endif
359 
360  if (solve_adjoint)
361  {
362  std::vector<unsigned int> adjs(system_list.size(),
364  // Set up proper initial guesses
365  for (unsigned int i=0; i != system_list.size(); ++i)
366  {
367  System *sys = system_list[i];
368  libmesh_assert (solution_vectors->find(sys) !=
369  solution_vectors->end());
370  const NumericVector<Number> *vec = solution_vectors->find(sys)->second;
371  for (unsigned int j=0; j != sys->qoi.size(); ++j)
372  {
373  std::ostringstream adjoint_name;
374  adjoint_name << "adjoint_solution" << j;
375 
376  if (vec == sys->request_vector(adjoint_name.str()))
377  {
378  adjs[i] = j;
379  break;
380  }
381  }
382  libmesh_assert_not_equal_to (adjs[i], libMesh::invalid_uint);
383  system_list[i]->get_adjoint_solution(adjs[i]) =
384  *system_list[i]->solution;
385  }
386 
387  es.adjoint_solve();
388 
389  // Put the adjoint_solution into solution for
390  // comparisons
391  for (unsigned int i=0; i != system_list.size(); ++i)
392  {
393  system_list[i]->get_adjoint_solution(adjs[i]).swap(*system_list[i]->solution);
394  system_list[i]->update();
395  }
396  }
397  else
398  es.solve();
399  }
400  }
401  else
402  {
403  System *sys = system_list[0];
404 
405  // Even if we had a decent preconditioner, valid matrix etc. before
406  // refinement, we don't any more.
407  sys->disable_cache();
408 
409  // No specified vectors == forward solve
410  if (!solution_vectors)
411  sys->solve();
412  else
413  {
414  libmesh_assert (solution_vectors->find(sys) !=
415  solution_vectors->end());
416 
417  const NumericVector<Number> *vec = solution_vectors->find(sys)->second;
418 
419  libmesh_assert(solve_adjoint ||
420  (solution_vectors->find(sys)->second ==
421  sys->solution.get()) ||
422  !solution_vectors->find(sys)->second);
423 
424  if (solve_adjoint)
425  {
426  unsigned int adj = libMesh::invalid_uint;
427  for (unsigned int j=0; j != sys->qoi.size(); ++j)
428  {
429  std::ostringstream adjoint_name;
430  adjoint_name << "adjoint_solution" << j;
431 
432  if (vec == sys->request_vector(adjoint_name.str()))
433  {
434  adj = j;
435  break;
436  }
437  }
438  libmesh_assert_not_equal_to (adj, libMesh::invalid_uint);
439 
440  // Set up proper initial guess
441  sys->get_adjoint_solution(adj) = *sys->solution;
442  sys->adjoint_solve();
443  // Put the adjoint_solution into solution for
444  // comparisons
445  sys->get_adjoint_solution(adj).swap(*sys->solution);
446  sys->update();
447  }
448  else
449  sys->solve();
450  }
451  }
452 
453  // Get the error in the uniformly refined solution(s).
454 
455  for (unsigned int sysnum=0; sysnum != system_list.size(); ++sysnum)
456  {
457  System &system = *system_list[sysnum];
458 
459  unsigned int n_vars = system.n_vars();
460 
461  DofMap &dof_map = system.get_dof_map();
462 
463  const SystemNorm &system_i_norm =
464  _error_norms->find(&system)->second;
465 
466  NumericVector<Number> *projected_solution = projected_solutions[sysnum];
467 
468  // Loop over all the variables in the system
469  for (unsigned int var=0; var<n_vars; var++)
470  {
471  // Get the error vector to fill for this system and variable
472  ErrorVector *err_vec = error_per_cell;
473  if (!err_vec)
474  {
475  libmesh_assert(errors_per_cell);
476  err_vec =
477  (*errors_per_cell)[std::make_pair(&system,var)];
478  }
479 
480  // The type of finite element to use for this variable
481  const FEType& fe_type = dof_map.variable_type (var);
482 
483  // Finite element object for each fine element
484  AutoPtr<FEBase> fe (FEBase::build (dim, fe_type));
485 
486  // Build and attach an appropriate quadrature rule
487  AutoPtr<QBase> qrule = fe_type.default_quadrature_rule(dim);
488  fe->attach_quadrature_rule (qrule.get());
489 
490  const std::vector<Real>& JxW = fe->get_JxW();
491  const std::vector<std::vector<Real> >& phi = fe->get_phi();
492  const std::vector<std::vector<RealGradient> >& dphi =
493  fe->get_dphi();
494 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
495  const std::vector<std::vector<RealTensor> >& d2phi =
496  fe->get_d2phi();
497 #endif
498 
499  // The global DOF indices for the fine element
500  std::vector<dof_id_type> dof_indices;
501 
502  // Iterate over all the active elements in the fine mesh
503  // that live on this processor.
506 
507  for (; elem_it != elem_end; ++elem_it)
508  {
509  // e is necessarily an active element on the local processor
510  const Elem* elem = *elem_it;
511 
512  // Find the element id for the corresponding coarse grid element
513  const Elem* coarse = elem;
514  dof_id_type e_id = coarse->id();
515  while (e_id >= max_coarse_elem_id)
516  {
517  libmesh_assert (coarse->parent());
518  coarse = coarse->parent();
519  e_id = coarse->id();
520  }
521 
522  double L2normsq = 0., H1seminormsq = 0., H2seminormsq = 0.;
523 
524  // reinitialize the element-specific data
525  // for the current element
526  fe->reinit (elem);
527 
528  // Get the local to global degree of freedom maps
529  dof_map.dof_indices (elem, dof_indices, var);
530 
531  // The number of quadrature points
532  const unsigned int n_qp = qrule->n_points();
533 
534  // The number of shape functions
535  const unsigned int n_sf =
536  libmesh_cast_int<unsigned int>(dof_indices.size());
537 
538  //
539  // Begin the loop over the Quadrature points.
540  //
541  for (unsigned int qp=0; qp<n_qp; qp++)
542  {
543  Number u_fine = 0., u_coarse = 0.;
544 
545  Gradient grad_u_fine, grad_u_coarse;
546 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
547  Tensor grad2_u_fine, grad2_u_coarse;
548 #endif
549 
550  // Compute solution values at the current
551  // quadrature point. This reqiures a sum
552  // over all the shape functions evaluated
553  // at the quadrature point.
554  for (unsigned int i=0; i<n_sf; i++)
555  {
556  u_fine += phi[i][qp]*system.current_solution (dof_indices[i]);
557  u_coarse += phi[i][qp]*(*projected_solution) (dof_indices[i]);
558  grad_u_fine += dphi[i][qp]*system.current_solution (dof_indices[i]);
559  grad_u_coarse += dphi[i][qp]*(*projected_solution) (dof_indices[i]);
560 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
561  grad2_u_fine += d2phi[i][qp]*system.current_solution (dof_indices[i]);
562  grad2_u_coarse += d2phi[i][qp]*(*projected_solution) (dof_indices[i]);
563 #endif
564  }
565 
566  // Compute the value of the error at this quadrature point
567  const Number val_error = u_fine - u_coarse;
568 
569  // Add the squares of the error to each contribution
570  if (system_i_norm.type(var) == L2 ||
571  system_i_norm.type(var) == H1 ||
572  system_i_norm.type(var) == H2)
573  {
574  L2normsq += JxW[qp] * system_i_norm.weight_sq(var) *
575  TensorTools::norm_sq(val_error);
576  libmesh_assert_greater_equal (L2normsq, 0.);
577  }
578 
579 
580  // Compute the value of the error in the gradient at this
581  // quadrature point
582  if (system_i_norm.type(var) == H1 ||
583  system_i_norm.type(var) == H2 ||
584  system_i_norm.type(var) == H1_SEMINORM)
585  {
586  Gradient grad_error = grad_u_fine - grad_u_coarse;
587 
588  H1seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
589  grad_error.size_sq();
590  libmesh_assert_greater_equal (H1seminormsq, 0.);
591  }
592 
593 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
594  // Compute the value of the error in the hessian at this
595  // quadrature point
596  if (system_i_norm.type(var) == H2 ||
597  system_i_norm.type(var) == H2_SEMINORM)
598  {
599  Tensor grad2_error = grad2_u_fine - grad2_u_coarse;
600 
601  H2seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
602  grad2_error.size_sq();
603  libmesh_assert_greater_equal (H2seminormsq, 0.);
604  }
605 #endif
606  } // end qp loop
607 
608  if (system_i_norm.type(var) == L2 ||
609  system_i_norm.type(var) == H1 ||
610  system_i_norm.type(var) == H2)
611  (*err_vec)[e_id] +=
612  static_cast<ErrorVectorReal>(L2normsq);
613  if (system_i_norm.type(var) == H1 ||
614  system_i_norm.type(var) == H2 ||
615  system_i_norm.type(var) == H1_SEMINORM)
616  (*err_vec)[e_id] +=
617  static_cast<ErrorVectorReal>(H1seminormsq);
618 
619  if (system_i_norm.type(var) == H2 ||
620  system_i_norm.type(var) == H2_SEMINORM)
621  (*err_vec)[e_id] +=
622  static_cast<ErrorVectorReal>(H2seminormsq);
623  } // End loop over active local elements
624  } // End loop over variables
625 
626  // Don't bother projecting the solution; we'll restore from backup
627  // after coarsening
628  system.project_solution_on_reinit() = false;
629  }
630 
631 
632  // Uniformly coarsen the mesh, without projecting the solution
634 
635  for (unsigned int i = 0; i != number_h_refinements; ++i)
636  {
637  mesh_refinement.uniformly_coarsen(1);
638  // FIXME - should the reinits here be necessary? - RHS
639  es.reinit();
640  }
641 
642  for (unsigned int i = 0; i != number_p_refinements; ++i)
643  {
644  mesh_refinement.uniformly_p_coarsen(1);
645  es.reinit();
646  }
647 
648  // We should be back where we started
649  libmesh_assert_equal_to (n_coarse_elem, mesh.n_elem());
650 
651  // Each processor has now computed the error contribuions
652  // for its local elements. We need to sum the vector
653  // and then take the square-root of each component. Note
654  // that we only need to sum if we are running on multiple
655  // processors, and we only need to take the square-root
656  // if the value is nonzero. There will in general be many
657  // zeros for the inactive elements.
658 
659  if (error_per_cell)
660  {
661  // First sum the vector of estimated error values
662  this->reduce_error(*error_per_cell, es.comm());
663 
664  // Compute the square-root of each component.
665  START_LOG("std::sqrt()", "UniformRefinementEstimator");
666  for (unsigned int i=0; i<error_per_cell->size(); i++)
667  if ((*error_per_cell)[i] != 0.)
668  (*error_per_cell)[i] = std::sqrt((*error_per_cell)[i]);
669  STOP_LOG("std::sqrt()", "UniformRefinementEstimator");
670  }
671  else
672  {
673  for (ErrorMap::iterator it = errors_per_cell->begin();
674  it != errors_per_cell->end(); ++it)
675  {
676  ErrorVector *e = it->second;
677  // First sum the vector of estimated error values
678  this->reduce_error(*e, es.comm());
679 
680  // Compute the square-root of each component.
681  START_LOG("std::sqrt()", "UniformRefinementEstimator");
682  for (unsigned int i=0; i<e->size(); i++)
683  if ((*e)[i] != 0.)
684  (*e)[i] = std::sqrt((*e)[i]);
685  STOP_LOG("std::sqrt()", "UniformRefinementEstimator");
686  }
687  }
688 
689  // Restore old solutions and clean up the heap
690  for (unsigned int i=0; i != system_list.size(); ++i)
691  {
692  System &system = *system_list[i];
693 
694  system.project_solution_on_reinit() = old_projection_settings[i];
695 
696  // Restore the coarse solution vectors and delete their copies
697  *system.solution = *coarse_solutions[i];
698  delete coarse_solutions[i];
699  *system.current_local_solution = *coarse_local_solutions[i];
700  delete coarse_local_solutions[i];
701  delete projected_solutions[i];
702 
703  for (System::vectors_iterator vec = system.vectors_begin(); vec !=
704  system.vectors_end(); ++vec)
705  {
706  // The (string) name of this vector
707  const std::string& var_name = vec->first;
708 
709  system.get_vector(var_name) = *coarse_vectors[i][var_name];
710 
711  coarse_vectors[i][var_name]->clear();
712  delete coarse_vectors[i][var_name];
713  }
714  }
715 
716  // Restore old partitioner settings
717  mesh.partitioner() = old_partitioner;
718 }
719 
720 } // namespace libMesh
721 
722 #endif // #ifdef LIBMESH_ENABLE_AMR

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo