system_projection.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 
19 
20 // C++ includes
21 #include <vector>
22 
23 // Local includes
24 #include "libmesh/boundary_info.h"
25 #include "libmesh/dense_matrix.h"
26 #include "libmesh/dense_vector.h"
28 #include "libmesh/dof_map.h"
29 #include "libmesh/elem.h"
31 #include "libmesh/fe_base.h"
32 #include "libmesh/fe_interface.h"
34 #include "libmesh/mesh_base.h"
35 #include "libmesh/numeric_vector.h"
37 #include "libmesh/system.h"
38 #include "libmesh/threads.h"
40 
41 namespace libMesh
42 {
43 
44 // ------------------------------------------------------------
45 // Helper class definitions
46 
53  {
54  private:
55  const System &system;
58 
59  public:
60  ProjectVector (const System &system_in,
61  const NumericVector<Number> &old_v_in,
62  NumericVector<Number> &new_v_in) :
63  system(system_in),
64  old_vector(old_v_in),
65  new_vector(new_v_in)
66  {}
67 
68  void operator()(const ConstElemRange &range) const;
69  };
70 
71 
82  {
83  private:
84  const System &system;
85 
86  public:
87  BuildProjectionList (const System &system_in) :
88  system(system_in),
89  send_list()
90  {}
91 
93  system(other.system),
94  send_list()
95  {}
96 
97  void unique();
98  void operator()(const ConstElemRange &range);
99  void join (const BuildProjectionList &other);
100  std::vector<dof_id_type> send_list;
101  };
102 
103 
104 
111  {
112  private:
113  const System &system;
114 
119 
120  public:
121  ProjectSolution (const System &system_in,
122  FunctionBase<Number>* f_in,
124  const Parameters &parameters_in,
125  NumericVector<Number> &new_v_in) :
126  system(system_in),
127  f(f_in ? f_in->clone() : AutoPtr<FunctionBase<Number> >(NULL)),
128  g(g_in ? g_in->clone() : AutoPtr<FunctionBase<Gradient> >(NULL)),
129  parameters(parameters_in),
130  new_vector(new_v_in)
131  {
132  libmesh_assert(f.get());
133  f->init();
134  if (g.get())
135  g->init();
136  }
137 
139  system(in.system),
140  f(in.f.get() ? in.f->clone() : AutoPtr<FunctionBase<Number> >(NULL)),
141  g(in.g.get() ? in.g->clone() : AutoPtr<FunctionBase<Gradient> >(NULL)),
144  {
145  libmesh_assert(f.get());
146  f->init();
147  if (g.get())
148  g->init();
149  }
150 
151  void operator()(const ConstElemRange &range) const;
152  };
153 
154 
161  {
162  private:
163  const System &system;
164 
168 
169  public:
170  ProjectFEMSolution (const System &system_in,
173  NumericVector<Number> &new_v_in) :
174  system(system_in),
175  f(f_in ? f_in->clone() : AutoPtr<FEMFunctionBase<Number> >(NULL)),
176  g(g_in ? g_in->clone() : AutoPtr<FEMFunctionBase<Gradient> >(NULL)),
177  new_vector(new_v_in)
178  {
179  libmesh_assert(f.get());
180  }
181 
183  system(in.system),
184  f(in.f.get() ? in.f->clone() : AutoPtr<FEMFunctionBase<Number> >(NULL)),
185  g(in.g.get() ? in.g->clone() : AutoPtr<FEMFunctionBase<Gradient> >(NULL)),
187  {
188  libmesh_assert(f.get());
189  }
190 
191  void operator()(const ConstElemRange &range) const;
192  };
193 
194 
201  {
202  private:
203  const std::set<boundary_id_type> &b;
204  const std::vector<unsigned int> &variables;
205  const System &system;
210 
211  public:
212  BoundaryProjectSolution (const std::set<boundary_id_type> &b_in,
213  const std::vector<unsigned int> &variables_in,
214  const System &system_in,
215  FunctionBase<Number>* f_in,
217  const Parameters &parameters_in,
218  NumericVector<Number> &new_v_in) :
219  b(b_in),
220  variables(variables_in),
221  system(system_in),
222  f(f_in ? f_in->clone() : AutoPtr<FunctionBase<Number> >(NULL)),
223  g(g_in ? g_in->clone() : AutoPtr<FunctionBase<Gradient> >(NULL)),
224  parameters(parameters_in),
225  new_vector(new_v_in)
226  {
227  libmesh_assert(f.get());
228  f->init();
229  if (g.get())
230  g->init();
231  }
232 
234  b(in.b),
235  variables(in.variables),
236  system(in.system),
237  f(in.f.get() ? in.f->clone() : AutoPtr<FunctionBase<Number> >(NULL)),
238  g(in.g.get() ? in.g->clone() : AutoPtr<FunctionBase<Gradient> >(NULL)),
241  {
242  libmesh_assert(f.get());
243  f->init();
244  if (g.get())
245  g->init();
246  }
247 
248  void operator()(const ConstElemRange &range) const;
249  };
250 
251 
252 
253 // ------------------------------------------------------------
254 // System implementation
256 {
257  // Create a copy of the vector, which currently
258  // contains the old data.
260  old_vector (vector.clone());
261 
262  // Project the old vector to the new vector
263  this->project_vector (*old_vector, vector);
264 }
265 
266 
273  NumericVector<Number>& new_v) const
274 {
275  START_LOG ("project_vector()", "System");
276 
283  new_v.clear();
284 
285 #ifdef LIBMESH_ENABLE_AMR
286 
287  // Resize the new vector and get a serial version.
288  NumericVector<Number> *new_vector_ptr = NULL;
289  AutoPtr<NumericVector<Number> > new_vector_built;
290  NumericVector<Number> *local_old_vector;
291  AutoPtr<NumericVector<Number> > local_old_vector_built;
292  const NumericVector<Number> *old_vector_ptr = NULL;
293 
294  ConstElemRange active_local_elem_range
295  (this->get_mesh().active_local_elements_begin(),
296  this->get_mesh().active_local_elements_end());
297 
298  // If the old vector was uniprocessor, make the new
299  // vector uniprocessor
300  if (old_v.type() == SERIAL)
301  {
302  new_v.init (this->n_dofs(), false, SERIAL);
303  new_vector_ptr = &new_v;
304  old_vector_ptr = &old_v;
305  }
306 
307  // Otherwise it is a parallel, distributed vector, which
308  // we need to localize.
309  else if (old_v.type() == PARALLEL)
310  {
311  // Build a send list for efficient localization
312  BuildProjectionList projection_list(*this);
313  Threads::parallel_reduce (active_local_elem_range,
314  projection_list);
315 
316  // Create a sorted, unique send_list
317  projection_list.unique();
318 
319  new_v.init (this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
320  new_vector_built = NumericVector<Number>::build(this->comm());
321  local_old_vector_built = NumericVector<Number>::build(this->comm());
322  new_vector_ptr = new_vector_built.get();
323  local_old_vector = local_old_vector_built.get();
324  new_vector_ptr->init(this->n_dofs(), false, SERIAL);
325  local_old_vector->init(old_v.size(), false, SERIAL);
326  old_v.localize(*local_old_vector, projection_list.send_list);
327  local_old_vector->close();
328  old_vector_ptr = local_old_vector;
329  }
330  else if (old_v.type() == GHOSTED)
331  {
332  // Build a send list for efficient localization
333  BuildProjectionList projection_list(*this);
334  Threads::parallel_reduce (active_local_elem_range,
335  projection_list);
336 
337  // Create a sorted, unique send_list
338  projection_list.unique();
339 
340  new_v.init (this->n_dofs(), this->n_local_dofs(),
341  this->get_dof_map().get_send_list(), false, GHOSTED);
342 
343  local_old_vector_built = NumericVector<Number>::build(this->comm());
344  new_vector_ptr = &new_v;
345  local_old_vector = local_old_vector_built.get();
346  local_old_vector->init(old_v.size(), old_v.local_size(),
347  projection_list.send_list, false, GHOSTED);
348  old_v.localize(*local_old_vector, projection_list.send_list);
349  local_old_vector->close();
350  old_vector_ptr = local_old_vector;
351  }
352  else // unknown old_v.type()
353  {
354  libMesh::err << "ERROR: Unknown old_v.type() == " << old_v.type()
355  << std::endl;
356  libmesh_error();
357  }
358 
359  // Note that the above will have zeroed the new_vector.
360  // Just to be sure, assert that new_vector_ptr and old_vector_ptr
361  // were successfully set before trying to deref them.
362  libmesh_assert(new_vector_ptr);
363  libmesh_assert(old_vector_ptr);
364 
365  NumericVector<Number> &new_vector = *new_vector_ptr;
366  const NumericVector<Number> &old_vector = *old_vector_ptr;
367 
368  Threads::parallel_for (active_local_elem_range,
369  ProjectVector(*this,
370  old_vector,
371  new_vector)
372  );
373 
374  // Copy the SCALAR dofs from old_vector to new_vector
375  // Note: We assume that all SCALAR dofs are on the
376  // processor with highest ID
377  if(this->processor_id() == (this->n_processors()-1))
378  {
379  const DofMap& dof_map = this->get_dof_map();
380  for (unsigned int var=0; var<this->n_vars(); var++)
381  if(this->variable(var).type().family == SCALAR)
382  {
383  // We can just map SCALAR dofs directly across
384  std::vector<dof_id_type> new_SCALAR_indices, old_SCALAR_indices;
385  dof_map.SCALAR_dof_indices (new_SCALAR_indices, var, false);
386  dof_map.SCALAR_dof_indices (old_SCALAR_indices, var, true);
387  const unsigned int new_n_dofs =
388  libmesh_cast_int<unsigned int>(new_SCALAR_indices.size());
389 
390  for (unsigned int i=0; i<new_n_dofs; i++)
391  {
392  new_vector.set( new_SCALAR_indices[i], old_vector(old_SCALAR_indices[i]) );
393  }
394  }
395  }
396 
397  new_vector.close();
398 
399  // If the old vector was serial, we probably need to send our values
400  // to other processors
401  //
402  // FIXME: I'm not sure how to make a NumericVector do that without
403  // creating a temporary parallel vector to use localize! - RHS
404  if (old_v.type() == SERIAL)
405  {
407  dist_v->init(this->n_dofs(), this->n_local_dofs(), false, PARALLEL);
408  dist_v->close();
409 
410  for (dof_id_type i=0; i!=dist_v->size(); i++)
411  if (new_vector(i) != 0.0)
412  dist_v->set(i, new_vector(i));
413 
414  dist_v->close();
415 
416  dist_v->localize (new_v, this->get_dof_map().get_send_list());
417  new_v.close();
418  }
419  // If the old vector was parallel, we need to update it
420  // and free the localized copies
421  else if (old_v.type() == PARALLEL)
422  {
423  // We may have to set dof values that this processor doesn't
424  // own in certain special cases, like LAGRANGE FIRST or
425  // HERMITE THIRD elements on second-order meshes
426  for (dof_id_type i=0; i!=new_v.size(); i++)
427  if (new_vector(i) != 0.0)
428  new_v.set(i, new_vector(i));
429  new_v.close();
430  }
431 
432  this->get_dof_map().enforce_constraints_exactly(*this, &new_v);
433 
434 #else
435 
436  // AMR is disabled: simply copy the vector
437  new_v = old_v;
438 
439 #endif // #ifdef LIBMESH_ENABLE_AMR
440 
441  STOP_LOG("project_vector()", "System");
442 }
443 
444 
445 
451  const Parameters& parameters,
452  const std::string& sys_name,
453  const std::string& unknown_name),
454  Gradient gptr(const Point& p,
455  const Parameters& parameters,
456  const std::string& sys_name,
457  const std::string& unknown_name),
458  const Parameters& parameters) const
459 {
460  WrappedFunction<Number> f(*this, fptr, &parameters);
461  WrappedFunction<Gradient> g(*this, gptr, &parameters);
462  this->project_solution(&f, &g);
463 }
464 
465 
471  FunctionBase<Gradient> *g) const
472 {
473  this->project_vector(*solution, f, g);
474 
475  solution->localize(*current_local_solution, _dof_map->get_send_list());
476 }
477 
478 
484  FEMFunctionBase<Gradient> *g) const
485 {
486  this->project_vector(*solution, f, g);
487 
488  solution->localize(*current_local_solution, _dof_map->get_send_list());
489 }
490 
491 
496 void System::project_vector (Number fptr(const Point& p,
497  const Parameters& parameters,
498  const std::string& sys_name,
499  const std::string& unknown_name),
500  Gradient gptr(const Point& p,
501  const Parameters& parameters,
502  const std::string& sys_name,
503  const std::string& unknown_name),
504  const Parameters& parameters,
505  NumericVector<Number>& new_vector) const
506 {
507  WrappedFunction<Number> f(*this, fptr, &parameters);
508  WrappedFunction<Gradient> g(*this, gptr, &parameters);
509  this->project_vector(new_vector, &f, &g);
510 }
511 
518  FunctionBase<Gradient> *g) const
519 {
520  START_LOG ("project_vector()", "System");
521 
523  (ConstElemRange (this->get_mesh().active_local_elements_begin(),
524  this->get_mesh().active_local_elements_end() ),
525  ProjectSolution(*this, f, g,
526  this->get_equation_systems().parameters,
527  new_vector)
528  );
529 
530  // Also, load values into the SCALAR dofs
531  // Note: We assume that all SCALAR dofs are on the
532  // processor with highest ID
533  if(this->processor_id() == (this->n_processors()-1))
534  {
535  // We get different scalars as different
536  // components from a new-style f functor.
537  DenseVector<Number> fout(this->n_components());
538  bool filled_fout = false;
539 
540  const DofMap& dof_map = this->get_dof_map();
541  for (unsigned int var=0; var<this->n_vars(); var++)
542  if(this->variable(var).type().family == SCALAR)
543  {
544  if (!filled_fout)
545  {
546  (*f) (Point(), this->time, fout);
547  filled_fout = true;
548  }
549 
550  std::vector<dof_id_type> SCALAR_indices;
551  dof_map.SCALAR_dof_indices (SCALAR_indices, var);
552  const unsigned int n_SCALAR_dofs =
553  libmesh_cast_int<unsigned int>(SCALAR_indices.size());
554 
555  for (unsigned int i=0; i<n_SCALAR_dofs; i++)
556  {
557  const dof_id_type global_index = SCALAR_indices[i];
558  const unsigned int component_index =
559  this->variable_scalar_number(var,i);
560  new_vector.set(global_index, fout(component_index));
561  }
562  }
563  }
564 
565  new_vector.close();
566 
567 #ifdef LIBMESH_ENABLE_CONSTRAINTS
568  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
569 #endif
570 
571  STOP_LOG("project_vector()", "System");
572 }
573 
574 
581  FEMFunctionBase<Gradient> *g) const
582 {
583  START_LOG ("project_fem_vector()", "System");
584 
586  (ConstElemRange (this->get_mesh().active_local_elements_begin(),
587  this->get_mesh().active_local_elements_end() ),
588  ProjectFEMSolution(*this, f, g, new_vector)
589  );
590 
591  // Also, load values into the SCALAR dofs
592  // Note: We assume that all SCALAR dofs are on the
593  // processor with highest ID
594  if(this->processor_id() == (this->n_processors()-1))
595  {
596  // FIXME: Do we want to first check for SCALAR vars before building this? [PB]
597  FEMContext context( *this );
598 
599  const DofMap& dof_map = this->get_dof_map();
600  for (unsigned int var=0; var<this->n_vars(); var++)
601  if(this->variable(var).type().family == SCALAR)
602  {
603  // FIXME: We reinit with an arbitrary element in case the user
604  // doesn't override FEMFunctionBase::component. Is there
605  // any use case we're missing? [PB]
606  Elem *el = const_cast<Elem *>(*(this->get_mesh().active_local_elements_begin()));
607  context.pre_fe_reinit( *this, el );
608 
609  std::vector<dof_id_type> SCALAR_indices;
610  dof_map.SCALAR_dof_indices (SCALAR_indices, var);
611  const unsigned int n_SCALAR_dofs =
612  libmesh_cast_int<unsigned int>(SCALAR_indices.size());
613 
614  for (unsigned int i=0; i<n_SCALAR_dofs; i++)
615  {
616  const dof_id_type global_index = SCALAR_indices[i];
617  const unsigned int component_index =
618  this->variable_scalar_number(var,i);
619 
620  new_vector.set(global_index, f->component(context, component_index, Point(), this->time));
621  }
622  }
623  }
624 
625  new_vector.close();
626 
627 #ifdef LIBMESH_ENABLE_CONSTRAINTS
628  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
629 #endif
630 
631  STOP_LOG("project_fem_vector()", "System");
632 }
633 
634 
641  (const std::set<boundary_id_type> &b,
642  const std::vector<unsigned int> &variables,
643  Number fptr(const Point& p,
644  const Parameters& parameters,
645  const std::string& sys_name,
646  const std::string& unknown_name),
647  Gradient gptr(const Point& p,
648  const Parameters& parameters,
649  const std::string& sys_name,
650  const std::string& unknown_name),
651  const Parameters& parameters)
652 {
653  WrappedFunction<Number> f(*this, fptr, &parameters);
654  WrappedFunction<Gradient> g(*this, gptr, &parameters);
655  this->boundary_project_solution(b, variables, &f, &g);
656 }
657 
658 
665  (const std::set<boundary_id_type> &b,
666  const std::vector<unsigned int> &variables,
669 {
670  this->boundary_project_vector(b, variables, *solution, f, g);
671 
672  solution->localize(*current_local_solution);
673 }
674 
675 
676 
677 
678 
684  (const std::set<boundary_id_type> &b,
685  const std::vector<unsigned int> &variables,
686  Number fptr(const Point& p,
687  const Parameters& parameters,
688  const std::string& sys_name,
689  const std::string& unknown_name),
690  Gradient gptr(const Point& p,
691  const Parameters& parameters,
692  const std::string& sys_name,
693  const std::string& unknown_name),
694  const Parameters& parameters,
695  NumericVector<Number>& new_vector) const
696 {
697  WrappedFunction<Number> f(*this, fptr, &parameters);
698  WrappedFunction<Gradient> g(*this, gptr, &parameters);
699  this->boundary_project_vector(b, variables, new_vector, &f, &g);
700 }
701 
707  (const std::set<boundary_id_type> &b,
708  const std::vector<unsigned int> &variables,
709  NumericVector<Number>& new_vector,
711  FunctionBase<Gradient> *g) const
712 {
713  START_LOG ("boundary_project_vector()", "System");
714 
716  (ConstElemRange (this->get_mesh().active_local_elements_begin(),
717  this->get_mesh().active_local_elements_end() ),
718  BoundaryProjectSolution(b, variables, *this, f, g,
719  this->get_equation_systems().parameters,
720  new_vector)
721  );
722 
723  // We don't do SCALAR dofs when just projecting the boundary, so
724  // we're done here.
725 
726  new_vector.close();
727 
728 #ifdef LIBMESH_ENABLE_CONSTRAINTS
729  this->get_dof_map().enforce_constraints_exactly(*this, &new_vector);
730 #endif
731 
732  STOP_LOG("boundary_project_vector()", "System");
733 }
734 
735 
736 
737 #ifndef LIBMESH_ENABLE_AMR
739 {
740  libmesh_error();
741 }
742 #else
743 void ProjectVector::operator()(const ConstElemRange &range) const
744 {
745  START_LOG ("operator()","ProjectVector");
746 
747  // A vector for Lagrange element interpolation, indicating if we
748  // have visited a DOF yet. Note that this is thread-local storage,
749  // hence shared DOFS that live on thread boundaries may be doubly
750  // computed. It is expected that this will be more efficient
751  // than locking a thread-global version of already_done, though.
752  //
753  // FIXME: we should use this for non-Lagrange coarsening, too
754  std::vector<bool> already_done (system.n_dofs(), false);
755 
756  // The number of variables in this system
757  const unsigned int n_variables = system.n_vars();
758 
759  // The dimensionality of the current mesh
760  const unsigned int dim = system.get_mesh().mesh_dimension();
761 
762  // The DofMap for this system
763  const DofMap& dof_map = system.get_dof_map();
764 
765  // The element matrix and RHS for projections.
766  // Note that Ke is always real-valued, whereas
767  // Fe may be complex valued if complex number
768  // support is enabled
771  // The new element coefficients
773 
774 
775  // Loop over all the variables in the system
776  for (unsigned int var=0; var<n_variables; var++)
777  {
778  const Variable& variable = dof_map.variable(var);
779 
780  const FEType& base_fe_type = variable.type();
781 
782  if (base_fe_type.family == SCALAR)
783  continue;
784 
785  // Get FE objects of the appropriate type
786  AutoPtr<FEBase> fe (FEBase::build(dim, base_fe_type));
787  AutoPtr<FEBase> fe_coarse (FEBase::build(dim, base_fe_type));
788 
789  // Create FE objects with potentially different p_level
790  FEType fe_type, temp_fe_type;
791 
792  // Prepare variables for non-Lagrange projection
793  AutoPtr<QBase> qrule (base_fe_type.default_quadrature_rule(dim));
794  AutoPtr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
795  AutoPtr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
796  std::vector<Point> coarse_qpoints;
797 
798  // The values of the shape functions at the quadrature
799  // points
800  const std::vector<std::vector<Real> >& phi_values =
801  fe->get_phi();
802  const std::vector<std::vector<Real> >& phi_coarse =
803  fe_coarse->get_phi();
804 
805  // The Jacobian * quadrature weight at the quadrature points
806  const std::vector<Real>& JxW =
807  fe->get_JxW();
808 
809  // The XYZ locations of the quadrature points on the
810  // child element
811  const std::vector<Point>& xyz_values =
812  fe->get_xyz();
813 
814 
815  // The global DOF indices
816  std::vector<dof_id_type> new_dof_indices, old_dof_indices;
817  // Side/edge local DOF indices
818  std::vector<unsigned int> new_side_dofs, old_side_dofs;
819 
820  // Iterate over the elements in the range
821  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
822  {
823  const Elem* elem = *elem_it;
824  // If this element doesn't have an old_dof_object with dofs for the
825  // current system, then it must be newly added, so the user
826  // is responsible for setting the new dofs.
827 
828  // ... but we need a better way to test for that; the code
829  // below breaks on any FE type for which the elem stores no
830  // dofs.
831  // if (!elem->old_dof_object || !elem->old_dof_object->has_dofs(system.number()))
832  // continue;
833  const Elem* parent = elem->parent();
834 
835  // Per-subdomain variables don't need to be projected on
836  // elements where they're not active
837  if (!variable.active_on_subdomain(elem->subdomain_id()))
838  continue;
839 
840  // Adjust the FE type for p-refined elements
841  fe_type = base_fe_type;
842  fe_type.order = static_cast<Order>(fe_type.order +
843  elem->p_level());
844 
845  // We may need to remember the parent's p_level
846  unsigned int old_parent_level = 0;
847 
848  // Update the DOF indices for this element based on
849  // the new mesh
850  dof_map.dof_indices (elem, new_dof_indices, var);
851 
852  // The number of DOFs on the new element
853  const unsigned int new_n_dofs =
854  libmesh_cast_int<unsigned int>(new_dof_indices.size());
855 
856  // Fixed vs. free DoFs on edge/face projections
857  std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
858  std::vector<int> free_dof(new_n_dofs, 0);
859 
860  // The element type
861  const ElemType elem_type = elem->type();
862 
863  // The number of nodes on the new element
864  const unsigned int n_nodes = elem->n_nodes();
865 
866  // Zero the interpolated values
867  Ue.resize (new_n_dofs); Ue.zero();
868 
869  // Update the DOF indices based on the old mesh.
870  // This is done in one of three ways:
871  // 1.) If the child was just refined then it was not
872  // active in the previous mesh & hence has no solution
873  // values on it. In this case simply project or
874  // interpolate the solution from the parent, who was
875  // active in the previous mesh
876  // 2.) If the child was just coarsened, obtaining a
877  // well-defined solution may require doing independent
878  // projections on nodes, edges, faces, and interiors
879  // 3.) If the child was active in the previous
880  // mesh, we can just copy coefficients directly
881  if (elem->refinement_flag() == Elem::JUST_REFINED)
882  {
883  libmesh_assert(parent);
884  old_parent_level = parent->p_level();
885 
886  // We may have done p refinement or coarsening as well;
887  // if so then we need to reset the parent's p level
888  // so we can get the right DoFs from it
889  if (elem->p_refinement_flag() == Elem::JUST_REFINED)
890  {
891  libmesh_assert_greater (elem->p_level(), 0);
892  (const_cast<Elem *>(parent))->hack_p_level(elem->p_level() - 1);
893  }
894  else if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
895  {
896  (const_cast<Elem *>(parent))->hack_p_level(elem->p_level() + 1);
897  }
898 
899  dof_map.old_dof_indices (parent, old_dof_indices, var);
900  }
901  else
902  {
903  dof_map.old_dof_indices (elem, old_dof_indices, var);
904 
905  if (elem->p_refinement_flag() == Elem::DO_NOTHING)
906  libmesh_assert_equal_to (old_dof_indices.size(), new_n_dofs);
907 
908  // if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
909  // libmesh_assert (elem->has_children());
910  }
911 
912  unsigned int old_n_dofs =
913  libmesh_cast_int<unsigned int>(old_dof_indices.size());
914 
915  if (fe_type.family != LAGRANGE) {
916 
917  // For refined non-Lagrange elements, we do an L2
918  // projection
919  // FIXME: this will be a suboptimal and ill-defined
920  // result if we're using non-nested finite element
921  // spaces or if we're on a p-coarsened element!
922  if (elem->refinement_flag() == Elem::JUST_REFINED)
923  {
924  // Update the fe object based on the current child
925  fe->attach_quadrature_rule (qrule.get());
926  fe->reinit (elem);
927 
928  // The number of quadrature points on the child
929  const unsigned int n_qp = qrule->n_points();
930 
931  FEInterface::inverse_map (dim, fe_type, parent,
932  xyz_values, coarse_qpoints);
933 
934  fe_coarse->reinit(parent, &coarse_qpoints);
935 
936  // Reinitialize the element matrix and vector for
937  // the current element. Note that this will zero them
938  // before they are summed.
939  Ke.resize (new_n_dofs, new_n_dofs); Ke.zero();
940  Fe.resize (new_n_dofs); Fe.zero();
941 
942 
943  // Loop over the quadrature points
944  for (unsigned int qp=0; qp<n_qp; qp++)
945  {
946  // The solution value at the quadrature point
947  Number val = libMesh::zero;
948 
949  // Sum the function values * the DOF values
950  // at the point of interest to get the function value
951  // (Note that the # of DOFs on the parent need not be the
952  // same as on the child!)
953  for (unsigned int i=0; i<old_n_dofs; i++)
954  {
955  val += (old_vector(old_dof_indices[i])*
956  phi_coarse[i][qp]);
957  }
958 
959  // Now \p val contains the solution value of variable
960  // \p var at the qp'th quadrature point on the child
961  // element. It has been interpolated from the parent
962  // in case the child was just refined. Next we will
963  // construct the L2-projection matrix for the element.
964 
965  // Construct the Mass Matrix
966  for (unsigned int i=0; i<new_n_dofs; i++)
967  for (unsigned int j=0; j<new_n_dofs; j++)
968  Ke(i,j) += JxW[qp]*phi_values[i][qp]*phi_values[j][qp];
969 
970  // Construct the RHS
971  for (unsigned int i=0; i<new_n_dofs; i++)
972  Fe(i) += JxW[qp]*phi_values[i][qp]*val;
973 
974  } // end qp loop
975 
976  Ke.cholesky_solve(Fe, Ue);
977 
978  // Fix up the parent's p level in case we changed it
979  (const_cast<Elem *>(parent))->hack_p_level(old_parent_level);
980  }
981  else if (elem->refinement_flag() == Elem::JUST_COARSENED)
982  {
984  elem, Ue, var, true);
985  }
986  // For unrefined uncoarsened elements, we just copy DoFs
987  else
988  {
989  // FIXME - I'm sure this function would be about half
990  // the size if anyone ever figures out how to improve
991  // the DofMap interface... - RHS
992  if (elem->p_refinement_flag() == Elem::JUST_REFINED)
993  {
994  libmesh_assert_greater (elem->p_level(), 0);
995  // P refinement of non-hierarchic bases will
996  // require a whole separate code path
997  libmesh_assert (fe->is_hierarchic());
998  temp_fe_type = fe_type;
999  temp_fe_type.order =
1000  static_cast<Order>(temp_fe_type.order - 1);
1001  unsigned int old_index = 0, new_index = 0;
1002  for (unsigned int n=0; n != n_nodes; ++n)
1003  {
1004  const unsigned int nc =
1005  FEInterface::n_dofs_at_node (dim, temp_fe_type,
1006  elem_type, n);
1007  for (unsigned int i=0; i != nc; ++i)
1008  {
1009  Ue(new_index + i) =
1010  old_vector(old_dof_indices[old_index++]);
1011  }
1012  new_index +=
1013  FEInterface::n_dofs_at_node (dim, fe_type,
1014  elem_type, n);
1015  }
1016  const unsigned int nc =
1017  FEInterface::n_dofs_per_elem (dim, temp_fe_type,
1018  elem_type);
1019  for (unsigned int i=0; i != nc; ++i)
1020  {
1021  Ue(new_index++) =
1022  old_vector(old_dof_indices[old_index+i]);
1023  }
1024  }
1025  else if (elem->p_refinement_flag() ==
1027  {
1028  // P coarsening of non-hierarchic bases will
1029  // require a whole separate code path
1030  libmesh_assert (fe->is_hierarchic());
1031  temp_fe_type = fe_type;
1032  temp_fe_type.order =
1033  static_cast<Order>(temp_fe_type.order + 1);
1034  unsigned int old_index = 0, new_index = 0;
1035  for (unsigned int n=0; n != n_nodes; ++n)
1036  {
1037  const unsigned int nc =
1038  FEInterface::n_dofs_at_node (dim, fe_type,
1039  elem_type, n);
1040  for (unsigned int i=0; i != nc; ++i)
1041  {
1042  Ue(new_index++) =
1043  old_vector(old_dof_indices[old_index+i]);
1044  }
1045  old_index +=
1046  FEInterface::n_dofs_at_node (dim, temp_fe_type,
1047  elem_type, n);
1048  }
1049  const unsigned int nc =
1050  FEInterface::n_dofs_per_elem (dim, fe_type,
1051  elem_type);
1052  for (unsigned int i=0; i != nc; ++i)
1053  {
1054  Ue(new_index++) =
1055  old_vector(old_dof_indices[old_index+i]);
1056  }
1057  }
1058  else
1059  // If there's no p refinement, we can copy every DoF
1060  for (unsigned int i=0; i<new_n_dofs; i++)
1061  Ue(i) = old_vector(old_dof_indices[i]);
1062  }
1063  }
1064  else { // fe type is Lagrange
1065  // Loop over the DOFs on the element
1066  for (unsigned int new_local_dof=0;
1067  new_local_dof<new_n_dofs; new_local_dof++)
1068  {
1069  // The global DOF number for this local DOF
1070  const dof_id_type new_global_dof =
1071  new_dof_indices[new_local_dof];
1072 
1073  // The global DOF might lie outside of the bounds of a
1074  // distributed vector. Check for that and possibly continue
1075  // the loop
1076  if ((new_global_dof < new_vector.first_local_index()) ||
1077  (new_global_dof >= new_vector.last_local_index()))
1078  continue;
1079 
1080  // We might have already computed the solution for this DOF.
1081  // This is likely in the case of a shared node, particularly
1082  // at the corners of an element. Check to see if that is the
1083  // case
1084  if (already_done[new_global_dof] == true)
1085  continue;
1086 
1087  already_done[new_global_dof] = true;
1088 
1089  if (elem->refinement_flag() == Elem::JUST_REFINED ||
1091  {
1092  // The location of the child's node on the parent element
1093  const Point point =
1094  FEInterface::inverse_map (dim, fe_type, parent,
1095  elem->point(new_local_dof));
1096 
1097  // Sum the function values * the DOF values
1098  // at the point of interest to get the function value
1099  // (Note that the # of DOFs on the parent need not be the
1100  // same as on the child!)
1101  for (unsigned int old_local_dof=0;
1102  old_local_dof<old_n_dofs; old_local_dof++)
1103  {
1104  const dof_id_type old_global_dof =
1105  old_dof_indices[old_local_dof];
1106 
1107  Ue(new_local_dof) +=
1108  (old_vector(old_global_dof)*
1109  FEInterface::shape(dim, base_fe_type, parent,
1110  old_local_dof, point));
1111  }
1112  }
1113  else
1114  {
1115  // Get the old global DOF index
1116  const dof_id_type old_global_dof =
1117  old_dof_indices[new_local_dof];
1118 
1119  Ue(new_local_dof) = old_vector(old_global_dof);
1120  }
1121  } // end local DOF loop
1122 
1123  // We may have to clean up a parent's p_level
1124  if (elem->refinement_flag() == Elem::JUST_REFINED)
1125  (const_cast<Elem *>(parent))->hack_p_level(old_parent_level);
1126  } // end fe_type if()
1127 
1128  // Lock the new_vector since it is shared among threads.
1129  {
1130  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1131 
1132  for (unsigned int i = 0; i < new_n_dofs; i++)
1133  if (Ue(i) != 0.)
1134  new_vector.set(new_dof_indices[i], Ue(i));
1135  }
1136  } // end elem loop
1137  } // end variables loop
1138 
1139  STOP_LOG ("operator()","ProjectVector");
1140 }
1141 #endif // LIBMESH_ENABLE_AMR
1142 
1143 
1144 
1146 {
1147  // Sort the send list. After this duplicated
1148  // elements will be adjacent in the vector
1149  std::sort(this->send_list.begin(),
1150  this->send_list.end());
1151 
1152  // Now use std::unique to remove duplicate entries
1153  std::vector<dof_id_type>::iterator new_end =
1154  std::unique (this->send_list.begin(),
1155  this->send_list.end());
1156 
1157  // Remove the end of the send_list. Use the "swap trick"
1158  // from Effective STL
1159  std::vector<dof_id_type>
1160  (this->send_list.begin(), new_end).swap (this->send_list);
1161 }
1162 
1163 
1164 
1165 #ifndef LIBMESH_ENABLE_AMR
1167 {
1168  libmesh_error();
1169 }
1170 #else
1172 {
1173  // The DofMap for this system
1174  const DofMap& dof_map = system.get_dof_map();
1175 
1176  const dof_id_type first_old_dof = dof_map.first_old_dof();
1177  const dof_id_type end_old_dof = dof_map.end_old_dof();
1178 
1179  // We can handle all the variables at once.
1180  // The old global DOF indices
1181  std::vector<dof_id_type> di;
1182 
1183  // Iterate over the elements in the range
1184  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
1185  {
1186  const Elem* elem = *elem_it;
1187  // If this element doesn't have an old_dof_object with dofs for the
1188  // current system, then it must be newly added, so the user
1189  // is responsible for setting the new dofs.
1190 
1191  // ... but we need a better way to test for that; the code
1192  // below breaks on any FE type for which the elem stores no
1193  // dofs.
1194  // if (!elem->old_dof_object || !elem->old_dof_object->has_dofs(system.number()))
1195  // continue;
1196  const Elem* parent = elem->parent();
1197 
1198  if (elem->refinement_flag() == Elem::JUST_REFINED)
1199  {
1200  libmesh_assert(parent);
1201  unsigned int old_parent_level = parent->p_level();
1202 
1203  if (elem->p_refinement_flag() == Elem::JUST_REFINED)
1204  {
1205  // We may have done p refinement or coarsening as well;
1206  // if so then we need to reset the parent's p level
1207  // so we can get the right DoFs from it
1208  libmesh_assert_greater (elem->p_level(), 0);
1209  (const_cast<Elem *>(parent))->hack_p_level(elem->p_level() - 1);
1210  }
1211  else if (elem->p_refinement_flag() == Elem::JUST_COARSENED)
1212  {
1213  (const_cast<Elem *>(parent))->hack_p_level(elem->p_level() + 1);
1214  }
1215 
1216  dof_map.old_dof_indices (parent, di);
1217 
1218  // Fix up the parent's p level in case we changed it
1219  (const_cast<Elem *>(parent))->hack_p_level(old_parent_level);
1220  }
1221  else if (elem->refinement_flag() == Elem::JUST_COARSENED)
1222  {
1223  std::vector<dof_id_type> di_child;
1224  di.clear();
1225  for (unsigned int c=0; c != elem->n_children(); ++c)
1226  {
1227  dof_map.old_dof_indices (elem->child(c), di_child);
1228  di.insert(di.end(), di_child.begin(), di_child.end());
1229  }
1230  }
1231  else
1232  dof_map.old_dof_indices (elem, di);
1233 
1234  for (unsigned int i=0; i != di.size(); ++i)
1235  if (di[i] < first_old_dof || di[i] >= end_old_dof)
1236  this->send_list.push_back(di[i]);
1237  } // end elem loop
1238 }
1239 #endif // LIBMESH_ENABLE_AMR
1240 
1241 
1242 
1244 {
1245  // Joining simply requires I add the dof indices from the other object
1246  this->send_list.insert(this->send_list.end(),
1247  other.send_list.begin(),
1248  other.send_list.end());
1249 }
1250 
1251 
1253 {
1254  // We need data to project
1255  libmesh_assert(f.get());
1256 
1264  // The number of variables in this system
1265  const unsigned int n_variables = system.n_vars();
1266 
1267  // The dimensionality of the current mesh
1268  const unsigned int dim = system.get_mesh().mesh_dimension();
1269 
1270  // The DofMap for this system
1271  const DofMap& dof_map = system.get_dof_map();
1272 
1273  // The element matrix and RHS for projections.
1274  // Note that Ke is always real-valued, whereas
1275  // Fe may be complex valued if complex number
1276  // support is enabled
1277  DenseMatrix<Real> Ke;
1279  // The new element coefficients
1281 
1282 
1283  // Loop over all the variables in the system
1284  for (unsigned int var=0; var<n_variables; var++)
1285  {
1286  const Variable& variable = dof_map.variable(var);
1287 
1288  const FEType& fe_type = variable.type();
1289 
1290  if (fe_type.family == SCALAR)
1291  continue;
1292 
1293  const unsigned int var_component =
1295 
1296  // Get FE objects of the appropriate type
1297  AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
1298 
1299  // Prepare variables for projection
1300  AutoPtr<QBase> qrule (fe_type.default_quadrature_rule(dim));
1301  AutoPtr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
1302  AutoPtr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));
1303 
1304  // The values of the shape functions at the quadrature
1305  // points
1306  const std::vector<std::vector<Real> >& phi = fe->get_phi();
1307 
1308  // The gradients of the shape functions at the quadrature
1309  // points on the child element.
1310  const std::vector<std::vector<RealGradient> > *dphi = NULL;
1311 
1312  const FEContinuity cont = fe->get_continuity();
1313 
1314  if (cont == C_ONE)
1315  {
1316  // We'll need gradient data for a C1 projection
1317  libmesh_assert(g.get());
1318 
1319  const std::vector<std::vector<RealGradient> >&
1320  ref_dphi = fe->get_dphi();
1321  dphi = &ref_dphi;
1322  }
1323 
1324  // The Jacobian * quadrature weight at the quadrature points
1325  const std::vector<Real>& JxW =
1326  fe->get_JxW();
1327 
1328  // The XYZ locations of the quadrature points
1329  const std::vector<Point>& xyz_values =
1330  fe->get_xyz();
1331 
1332  // The global DOF indices
1333  std::vector<dof_id_type> dof_indices;
1334  // Side/edge DOF indices
1335  std::vector<unsigned int> side_dofs;
1336 
1337  // Iterate over all the elements in the range
1338  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
1339  {
1340  const Elem* elem = *elem_it;
1341 
1342  // Per-subdomain variables don't need to be projected on
1343  // elements where they're not active
1344  if (!variable.active_on_subdomain(elem->subdomain_id()))
1345  continue;
1346 
1347  // Update the DOF indices for this element based on
1348  // the current mesh
1349  dof_map.dof_indices (elem, dof_indices, var);
1350 
1351  // The number of DOFs on the element
1352  const unsigned int n_dofs =
1353  libmesh_cast_int<unsigned int>(dof_indices.size());
1354 
1355  // Fixed vs. free DoFs on edge/face projections
1356  std::vector<char> dof_is_fixed(n_dofs, false); // bools
1357  std::vector<int> free_dof(n_dofs, 0);
1358 
1359  // The element type
1360  const ElemType elem_type = elem->type();
1361 
1362  // The number of nodes on the new element
1363  const unsigned int n_nodes = elem->n_nodes();
1364 
1365  // Zero the interpolated values
1366  Ue.resize (n_dofs); Ue.zero();
1367 
1368  // In general, we need a series of
1369  // projections to ensure a unique and continuous
1370  // solution. We start by interpolating nodes, then
1371  // hold those fixed and project edges, then
1372  // hold those fixed and project faces, then
1373  // hold those fixed and project interiors
1374 
1375  // Interpolate node values first
1376  unsigned int current_dof = 0;
1377  for (unsigned int n=0; n!= n_nodes; ++n)
1378  {
1379  // FIXME: this should go through the DofMap,
1380  // not duplicate dof_indices code badly!
1381  const unsigned int nc =
1382  FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
1383  n);
1384  if (!elem->is_vertex(n))
1385  {
1386  current_dof += nc;
1387  continue;
1388  }
1389  if (cont == DISCONTINUOUS)
1390  {
1391  libmesh_assert_equal_to (nc, 0);
1392  }
1393  // Assume that C_ZERO elements have a single nodal
1394  // value shape function
1395  else if (cont == C_ZERO)
1396  {
1397  libmesh_assert_equal_to (nc, 1);
1398  Ue(current_dof) = f->component(var_component,
1399  elem->point(n),
1400  system.time);
1401  dof_is_fixed[current_dof] = true;
1402  current_dof++;
1403  }
1404  // The hermite element vertex shape functions are weird
1405  else if (fe_type.family == HERMITE)
1406  {
1407  Ue(current_dof) = f->component(var_component,
1408  elem->point(n),
1409  system.time);
1410  dof_is_fixed[current_dof] = true;
1411  current_dof++;
1412  Gradient grad = g->component(var_component,
1413  elem->point(n),
1414  system.time);
1415  // x derivative
1416  Ue(current_dof) = grad(0);
1417  dof_is_fixed[current_dof] = true;
1418  current_dof++;
1419  if (dim > 1)
1420  {
1421  // We'll finite difference mixed derivatives
1422  Point nxminus = elem->point(n),
1423  nxplus = elem->point(n);
1424  nxminus(0) -= TOLERANCE;
1425  nxplus(0) += TOLERANCE;
1426  Gradient gxminus = g->component(var_component,
1427  nxminus,
1428  system.time);
1429  Gradient gxplus = g->component(var_component,
1430  nxplus,
1431  system.time);
1432  // y derivative
1433  Ue(current_dof) = grad(1);
1434  dof_is_fixed[current_dof] = true;
1435  current_dof++;
1436  // xy derivative
1437  Ue(current_dof) = (gxplus(1) - gxminus(1))
1438  / 2. / TOLERANCE;
1439  dof_is_fixed[current_dof] = true;
1440  current_dof++;
1441 
1442  if (dim > 2)
1443  {
1444  // z derivative
1445  Ue(current_dof) = grad(2);
1446  dof_is_fixed[current_dof] = true;
1447  current_dof++;
1448  // xz derivative
1449  Ue(current_dof) = (gxplus(2) - gxminus(2))
1450  / 2. / TOLERANCE;
1451  dof_is_fixed[current_dof] = true;
1452  current_dof++;
1453  // We need new points for yz
1454  Point nyminus = elem->point(n),
1455  nyplus = elem->point(n);
1456  nyminus(1) -= TOLERANCE;
1457  nyplus(1) += TOLERANCE;
1458  Gradient gyminus = g->component(var_component,
1459  nyminus,
1460  system.time);
1461  Gradient gyplus = g->component(var_component,
1462  nyplus,
1463  system.time);
1464  // xz derivative
1465  Ue(current_dof) = (gyplus(2) - gyminus(2))
1466  / 2. / TOLERANCE;
1467  dof_is_fixed[current_dof] = true;
1468  current_dof++;
1469  // Getting a 2nd order xyz is more tedious
1470  Point nxmym = elem->point(n),
1471  nxmyp = elem->point(n),
1472  nxpym = elem->point(n),
1473  nxpyp = elem->point(n);
1474  nxmym(0) -= TOLERANCE;
1475  nxmym(1) -= TOLERANCE;
1476  nxmyp(0) -= TOLERANCE;
1477  nxmyp(1) += TOLERANCE;
1478  nxpym(0) += TOLERANCE;
1479  nxpym(1) -= TOLERANCE;
1480  nxpyp(0) += TOLERANCE;
1481  nxpyp(1) += TOLERANCE;
1482  Gradient gxmym = g->component(var_component,
1483  nxmym,
1484  system.time);
1485  Gradient gxmyp = g->component(var_component,
1486  nxmyp,
1487  system.time);
1488  Gradient gxpym = g->component(var_component,
1489  nxpym,
1490  system.time);
1491  Gradient gxpyp = g->component(var_component,
1492  nxpyp,
1493  system.time);
1494  Number gxzplus = (gxpyp(2) - gxmyp(2))
1495  / 2. / TOLERANCE;
1496  Number gxzminus = (gxpym(2) - gxmym(2))
1497  / 2. / TOLERANCE;
1498  // xyz derivative
1499  Ue(current_dof) = (gxzplus - gxzminus)
1500  / 2. / TOLERANCE;
1501  dof_is_fixed[current_dof] = true;
1502  current_dof++;
1503  }
1504  }
1505  }
1506  // Assume that other C_ONE elements have a single nodal
1507  // value shape function and nodal gradient component
1508  // shape functions
1509  else if (cont == C_ONE)
1510  {
1511  libmesh_assert_equal_to (nc, 1 + dim);
1512  Ue(current_dof) = f->component(var_component,
1513  elem->point(n),
1514  system.time);
1515  dof_is_fixed[current_dof] = true;
1516  current_dof++;
1517  Gradient grad = g->component(var_component,
1518  elem->point(n),
1519  system.time);
1520  for (unsigned int i=0; i!= dim; ++i)
1521  {
1522  Ue(current_dof) = grad(i);
1523  dof_is_fixed[current_dof] = true;
1524  current_dof++;
1525  }
1526  }
1527  else
1528  libmesh_error();
1529  }
1530 
1531  // In 3D, project any edge values next
1532  if (dim > 2 && cont != DISCONTINUOUS)
1533  for (unsigned int e=0; e != elem->n_edges(); ++e)
1534  {
1535  FEInterface::dofs_on_edge(elem, dim, fe_type, e,
1536  side_dofs);
1537 
1538  // Some edge dofs are on nodes and already
1539  // fixed, others are free to calculate
1540  unsigned int free_dofs = 0;
1541  for (unsigned int i=0; i != side_dofs.size(); ++i)
1542  if (!dof_is_fixed[side_dofs[i]])
1543  free_dof[free_dofs++] = i;
1544 
1545  // There may be nothing to project
1546  if (!free_dofs)
1547  continue;
1548 
1549  Ke.resize (free_dofs, free_dofs); Ke.zero();
1550  Fe.resize (free_dofs); Fe.zero();
1551  // The new edge coefficients
1552  DenseVector<Number> Uedge(free_dofs);
1553 
1554  // Initialize FE data on the edge
1555  fe->attach_quadrature_rule (qedgerule.get());
1556  fe->edge_reinit (elem, e);
1557  const unsigned int n_qp = qedgerule->n_points();
1558 
1559  // Loop over the quadrature points
1560  for (unsigned int qp=0; qp<n_qp; qp++)
1561  {
1562  // solution at the quadrature point
1563  Number fineval = f->component(var_component,
1564  xyz_values[qp],
1565  system.time);
1566  // solution grad at the quadrature point
1567  Gradient finegrad;
1568  if (cont == C_ONE)
1569  finegrad = g->component(var_component,
1570  xyz_values[qp],
1571  system.time);
1572 
1573  // Form edge projection matrix
1574  for (unsigned int sidei=0, freei=0;
1575  sidei != side_dofs.size(); ++sidei)
1576  {
1577  unsigned int i = side_dofs[sidei];
1578  // fixed DoFs aren't test functions
1579  if (dof_is_fixed[i])
1580  continue;
1581  for (unsigned int sidej=0, freej=0;
1582  sidej != side_dofs.size(); ++sidej)
1583  {
1584  unsigned int j = side_dofs[sidej];
1585  if (dof_is_fixed[j])
1586  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1587  JxW[qp] * Ue(j);
1588  else
1589  Ke(freei,freej) += phi[i][qp] *
1590  phi[j][qp] * JxW[qp];
1591  if (cont == C_ONE)
1592  {
1593  if (dof_is_fixed[j])
1594  Fe(freei) -= ((*dphi)[i][qp] *
1595  (*dphi)[j][qp]) *
1596  JxW[qp] * Ue(j);
1597  else
1598  Ke(freei,freej) += ((*dphi)[i][qp] *
1599  (*dphi)[j][qp])
1600  * JxW[qp];
1601  }
1602  if (!dof_is_fixed[j])
1603  freej++;
1604  }
1605  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1606  if (cont == C_ONE)
1607  Fe(freei) += (finegrad * (*dphi)[i][qp]) *
1608  JxW[qp];
1609  freei++;
1610  }
1611  }
1612 
1613  Ke.cholesky_solve(Fe, Uedge);
1614 
1615  // Transfer new edge solutions to element
1616  for (unsigned int i=0; i != free_dofs; ++i)
1617  {
1618  Number &ui = Ue(side_dofs[free_dof[i]]);
1620  std::abs(ui - Uedge(i)) < TOLERANCE);
1621  ui = Uedge(i);
1622  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1623  }
1624  }
1625 
1626  // Project any side values (edges in 2D, faces in 3D)
1627  if (dim > 1 && cont != DISCONTINUOUS)
1628  for (unsigned int s=0; s != elem->n_sides(); ++s)
1629  {
1630  FEInterface::dofs_on_side(elem, dim, fe_type, s,
1631  side_dofs);
1632 
1633  // Some side dofs are on nodes/edges and already
1634  // fixed, others are free to calculate
1635  unsigned int free_dofs = 0;
1636  for (unsigned int i=0; i != side_dofs.size(); ++i)
1637  if (!dof_is_fixed[side_dofs[i]])
1638  free_dof[free_dofs++] = i;
1639 
1640  // There may be nothing to project
1641  if (!free_dofs)
1642  continue;
1643 
1644  Ke.resize (free_dofs, free_dofs); Ke.zero();
1645  Fe.resize (free_dofs); Fe.zero();
1646  // The new side coefficients
1647  DenseVector<Number> Uside(free_dofs);
1648 
1649  // Initialize FE data on the side
1650  fe->attach_quadrature_rule (qsiderule.get());
1651  fe->reinit (elem, s);
1652  const unsigned int n_qp = qsiderule->n_points();
1653 
1654  // Loop over the quadrature points
1655  for (unsigned int qp=0; qp<n_qp; qp++)
1656  {
1657  // solution at the quadrature point
1658  Number fineval = f->component(var_component,
1659  xyz_values[qp],
1660  system.time);
1661  // solution grad at the quadrature point
1662  Gradient finegrad;
1663  if (cont == C_ONE)
1664  finegrad = g->component(var_component,
1665  xyz_values[qp],
1666  system.time);
1667 
1668  // Form side projection matrix
1669  for (unsigned int sidei=0, freei=0;
1670  sidei != side_dofs.size(); ++sidei)
1671  {
1672  unsigned int i = side_dofs[sidei];
1673  // fixed DoFs aren't test functions
1674  if (dof_is_fixed[i])
1675  continue;
1676  for (unsigned int sidej=0, freej=0;
1677  sidej != side_dofs.size(); ++sidej)
1678  {
1679  unsigned int j = side_dofs[sidej];
1680  if (dof_is_fixed[j])
1681  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1682  JxW[qp] * Ue(j);
1683  else
1684  Ke(freei,freej) += phi[i][qp] *
1685  phi[j][qp] * JxW[qp];
1686  if (cont == C_ONE)
1687  {
1688  if (dof_is_fixed[j])
1689  Fe(freei) -= ((*dphi)[i][qp] *
1690  (*dphi)[j][qp]) *
1691  JxW[qp] * Ue(j);
1692  else
1693  Ke(freei,freej) += ((*dphi)[i][qp] *
1694  (*dphi)[j][qp])
1695  * JxW[qp];
1696  }
1697  if (!dof_is_fixed[j])
1698  freej++;
1699  }
1700  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1701  if (cont == C_ONE)
1702  Fe(freei) += (finegrad * (*dphi)[i][qp]) *
1703  JxW[qp];
1704  freei++;
1705  }
1706  }
1707 
1708  Ke.cholesky_solve(Fe, Uside);
1709 
1710  // Transfer new side solutions to element
1711  for (unsigned int i=0; i != free_dofs; ++i)
1712  {
1713  Number &ui = Ue(side_dofs[free_dof[i]]);
1715  std::abs(ui - Uside(i)) < TOLERANCE);
1716  ui = Uside(i);
1717  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1718  }
1719  }
1720 
1721  // Project the interior values, finally
1722 
1723  // Some interior dofs are on nodes/edges/sides and
1724  // already fixed, others are free to calculate
1725  unsigned int free_dofs = 0;
1726  for (unsigned int i=0; i != n_dofs; ++i)
1727  if (!dof_is_fixed[i])
1728  free_dof[free_dofs++] = i;
1729 
1730  // There may be nothing to project
1731  if (free_dofs)
1732  {
1733 
1734  Ke.resize (free_dofs, free_dofs); Ke.zero();
1735  Fe.resize (free_dofs); Fe.zero();
1736  // The new interior coefficients
1737  DenseVector<Number> Uint(free_dofs);
1738 
1739  // Initialize FE data
1740  fe->attach_quadrature_rule (qrule.get());
1741  fe->reinit (elem);
1742  const unsigned int n_qp = qrule->n_points();
1743 
1744  // Loop over the quadrature points
1745  for (unsigned int qp=0; qp<n_qp; qp++)
1746  {
1747  // solution at the quadrature point
1748  Number fineval = f->component(var_component,
1749  xyz_values[qp],
1750  system.time);
1751  // solution grad at the quadrature point
1752  Gradient finegrad;
1753  if (cont == C_ONE)
1754  finegrad = g->component(var_component,
1755  xyz_values[qp],
1756  system.time);
1757 
1758  // Form interior projection matrix
1759  for (unsigned int i=0, freei=0; i != n_dofs; ++i)
1760  {
1761  // fixed DoFs aren't test functions
1762  if (dof_is_fixed[i])
1763  continue;
1764  for (unsigned int j=0, freej=0; j != n_dofs; ++j)
1765  {
1766  if (dof_is_fixed[j])
1767  Fe(freei) -= phi[i][qp] * phi[j][qp] * JxW[qp]
1768  * Ue(j);
1769  else
1770  Ke(freei,freej) += phi[i][qp] * phi[j][qp] *
1771  JxW[qp];
1772  if (cont == C_ONE)
1773  {
1774  if (dof_is_fixed[j])
1775  Fe(freei) -= ((*dphi)[i][qp] *
1776  (*dphi)[j][qp]) * JxW[qp] *
1777  Ue(j);
1778  else
1779  Ke(freei,freej) += ((*dphi)[i][qp] *
1780  (*dphi)[j][qp]) *
1781  JxW[qp];
1782  }
1783  if (!dof_is_fixed[j])
1784  freej++;
1785  }
1786  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1787  if (cont == C_ONE)
1788  Fe(freei) += (finegrad * (*dphi)[i][qp]) * JxW[qp];
1789  freei++;
1790  }
1791  }
1792  Ke.cholesky_solve(Fe, Uint);
1793 
1794  // Transfer new interior solutions to element
1795  for (unsigned int i=0; i != free_dofs; ++i)
1796  {
1797  Number &ui = Ue(free_dof[i]);
1799  std::abs(ui - Uint(i)) < TOLERANCE);
1800  ui = Uint(i);
1801  dof_is_fixed[free_dof[i]] = true;
1802  }
1803 
1804  } // if there are free interior dofs
1805 
1806  // Make sure every DoF got reached!
1807  for (unsigned int i=0; i != n_dofs; ++i)
1808  libmesh_assert(dof_is_fixed[i]);
1809 
1810  const dof_id_type
1811  first = new_vector.first_local_index(),
1812  last = new_vector.last_local_index();
1813 
1814  // Lock the new_vector since it is shared among threads.
1815  {
1816  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1817 
1818  for (unsigned int i = 0; i < n_dofs; i++)
1819  // We may be projecting a new zero value onto
1820  // an old nonzero approximation - RHS
1821  // if (Ue(i) != 0.)
1822  if ((dof_indices[i] >= first) &&
1823  (dof_indices[i] < last))
1824  {
1825  new_vector.set(dof_indices[i], Ue(i));
1826  }
1827  }
1828  } // end elem loop
1829  } // end variables loop
1830 }
1831 
1832 
1834 {
1835  // We need data to project
1836  libmesh_assert(f.get());
1837 
1845  FEMContext context( system );
1846 
1847  // The number of variables in this system
1848  const unsigned int n_variables = context.n_vars();
1849 
1850  // The dimensionality of the current mesh
1851  const unsigned int dim = context.get_dim();
1852 
1853  // The DofMap for this system
1854  const DofMap& dof_map = system.get_dof_map();
1855 
1856  // The element matrix and RHS for projections.
1857  // Note that Ke is always real-valued, whereas
1858  // Fe may be complex valued if complex number
1859  // support is enabled
1860  DenseMatrix<Real> Ke;
1862  // The new element coefficients
1864 
1865  // FIXME: Need to generalize this to vector-valued elements. [PB]
1866  FEBase* fe = NULL;
1867  FEBase* side_fe = NULL;
1868  FEBase* edge_fe = NULL;
1869 
1870  // First, loop over all variables and make sure the shape functions phi will
1871  // get built as well as the shape function gradients if the gradient functor
1872  // is supplied.
1873  for (unsigned int var=0; var<n_variables; var++)
1874  {
1875  context.get_element_fe( var, fe );
1876  if (fe->get_fe_type().family == SCALAR)
1877  continue;
1878  if( dim > 1 )
1879  context.get_side_fe( var, side_fe );
1880  if( dim > 2 )
1881  context.get_edge_fe( var, edge_fe );
1882 
1883  fe->get_phi();
1884  if( dim > 1 )
1885  side_fe->get_phi();
1886  if( dim > 2 )
1887  edge_fe->get_phi();
1888 
1889  const FEContinuity cont = fe->get_continuity();
1890  if(cont == C_ONE)
1891  {
1892  libmesh_assert(g.get());
1893  fe->get_dphi();
1894  side_fe->get_dphi();
1895  if( dim > 2 )
1896  edge_fe->get_dphi();
1897  }
1898  }
1899 
1900  // Now initialize any user requested shape functions
1901  f->init_context(context);
1902  if(g.get())
1903  g->init_context(context);
1904 
1905  std::vector<unsigned int> side_dofs;
1906 
1907  // Iterate over all the elements in the range
1908  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
1909  {
1910  const Elem* elem = *elem_it;
1911 
1912  context.pre_fe_reinit(system, elem);
1913 
1914  // Loop over all the variables in the system
1915  for (unsigned int var=0; var<n_variables; var++)
1916  {
1917  const Variable& variable = dof_map.variable(var);
1918 
1919  const FEType& fe_type = variable.type();
1920 
1921  if (fe_type.family == SCALAR)
1922  continue;
1923 
1924  // Per-subdomain variables don't need to be projected on
1925  // elements where they're not active
1926  if (!variable.active_on_subdomain(elem->subdomain_id()))
1927  continue;
1928 
1929  const FEContinuity cont = fe->get_continuity();
1930 
1931  const unsigned int var_component =
1933 
1934  const std::vector<dof_id_type>& dof_indices =
1935  context.get_dof_indices(var);
1936 
1937  // The number of DOFs on the element
1938  const unsigned int n_dofs =
1939  libmesh_cast_int<unsigned int>(dof_indices.size());
1940 
1941  // Fixed vs. free DoFs on edge/face projections
1942  std::vector<char> dof_is_fixed(n_dofs, false); // bools
1943  std::vector<int> free_dof(n_dofs, 0);
1944 
1945  // The element type
1946  const ElemType elem_type = elem->type();
1947 
1948  // The number of nodes on the new element
1949  const unsigned int n_nodes = elem->n_nodes();
1950 
1951  // Zero the interpolated values
1952  Ue.resize (n_dofs); Ue.zero();
1953 
1954  // In general, we need a series of
1955  // projections to ensure a unique and continuous
1956  // solution. We start by interpolating nodes, then
1957  // hold those fixed and project edges, then
1958  // hold those fixed and project faces, then
1959  // hold those fixed and project interiors
1960 
1961  // Interpolate node values first
1962  unsigned int current_dof = 0;
1963  for (unsigned int n=0; n!= n_nodes; ++n)
1964  {
1965  // FIXME: this should go through the DofMap,
1966  // not duplicate dof_indices code badly!
1967  const unsigned int nc =
1968  FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
1969  n);
1970  if (!elem->is_vertex(n))
1971  {
1972  current_dof += nc;
1973  continue;
1974  }
1975  if (cont == DISCONTINUOUS)
1976  {
1977  libmesh_assert_equal_to (nc, 0);
1978  }
1979  // Assume that C_ZERO elements have a single nodal
1980  // value shape function
1981  else if (cont == C_ZERO)
1982  {
1983  libmesh_assert_equal_to (nc, 1);
1984  Ue(current_dof) = f->component(context,
1985  var_component,
1986  elem->point(n),
1987  system.time);
1988  dof_is_fixed[current_dof] = true;
1989  current_dof++;
1990  }
1991  // The hermite element vertex shape functions are weird
1992  else if (fe_type.family == HERMITE)
1993  {
1994  Ue(current_dof) = f->component(context,
1995  var_component,
1996  elem->point(n),
1997  system.time);
1998  dof_is_fixed[current_dof] = true;
1999  current_dof++;
2000  Gradient grad = g->component(context,
2001  var_component,
2002  elem->point(n),
2003  system.time);
2004  // x derivative
2005  Ue(current_dof) = grad(0);
2006  dof_is_fixed[current_dof] = true;
2007  current_dof++;
2008  if (dim > 1)
2009  {
2010  // We'll finite difference mixed derivatives
2011  Point nxminus = elem->point(n),
2012  nxplus = elem->point(n);
2013  nxminus(0) -= TOLERANCE;
2014  nxplus(0) += TOLERANCE;
2015  Gradient gxminus = g->component(context,
2016  var_component,
2017  nxminus,
2018  system.time);
2019  Gradient gxplus = g->component(context,
2020  var_component,
2021  nxplus,
2022  system.time);
2023  // y derivative
2024  Ue(current_dof) = grad(1);
2025  dof_is_fixed[current_dof] = true;
2026  current_dof++;
2027  // xy derivative
2028  Ue(current_dof) = (gxplus(1) - gxminus(1))
2029  / 2. / TOLERANCE;
2030  dof_is_fixed[current_dof] = true;
2031  current_dof++;
2032 
2033  if (dim > 2)
2034  {
2035  // z derivative
2036  Ue(current_dof) = grad(2);
2037  dof_is_fixed[current_dof] = true;
2038  current_dof++;
2039  // xz derivative
2040  Ue(current_dof) = (gxplus(2) - gxminus(2))
2041  / 2. / TOLERANCE;
2042  dof_is_fixed[current_dof] = true;
2043  current_dof++;
2044  // We need new points for yz
2045  Point nyminus = elem->point(n),
2046  nyplus = elem->point(n);
2047  nyminus(1) -= TOLERANCE;
2048  nyplus(1) += TOLERANCE;
2049  Gradient gyminus = g->component(context,
2050  var_component,
2051  nyminus,
2052  system.time);
2053  Gradient gyplus = g->component(context,
2054  var_component,
2055  nyplus,
2056  system.time);
2057  // xz derivative
2058  Ue(current_dof) = (gyplus(2) - gyminus(2))
2059  / 2. / TOLERANCE;
2060  dof_is_fixed[current_dof] = true;
2061  current_dof++;
2062  // Getting a 2nd order xyz is more tedious
2063  Point nxmym = elem->point(n),
2064  nxmyp = elem->point(n),
2065  nxpym = elem->point(n),
2066  nxpyp = elem->point(n);
2067  nxmym(0) -= TOLERANCE;
2068  nxmym(1) -= TOLERANCE;
2069  nxmyp(0) -= TOLERANCE;
2070  nxmyp(1) += TOLERANCE;
2071  nxpym(0) += TOLERANCE;
2072  nxpym(1) -= TOLERANCE;
2073  nxpyp(0) += TOLERANCE;
2074  nxpyp(1) += TOLERANCE;
2075  Gradient gxmym = g->component(context,
2076  var_component,
2077  nxmym,
2078  system.time);
2079  Gradient gxmyp = g->component(context,
2080  var_component,
2081  nxmyp,
2082  system.time);
2083  Gradient gxpym = g->component(context,
2084  var_component,
2085  nxpym,
2086  system.time);
2087  Gradient gxpyp = g->component(context,
2088  var_component,
2089  nxpyp,
2090  system.time);
2091  Number gxzplus = (gxpyp(2) - gxmyp(2))
2092  / 2. / TOLERANCE;
2093  Number gxzminus = (gxpym(2) - gxmym(2))
2094  / 2. / TOLERANCE;
2095  // xyz derivative
2096  Ue(current_dof) = (gxzplus - gxzminus)
2097  / 2. / TOLERANCE;
2098  dof_is_fixed[current_dof] = true;
2099  current_dof++;
2100  }
2101  }
2102  }
2103  // Assume that other C_ONE elements have a single nodal
2104  // value shape function and nodal gradient component
2105  // shape functions
2106  else if (cont == C_ONE)
2107  {
2108  libmesh_assert_equal_to (nc, 1 + dim);
2109  Ue(current_dof) = f->component(context,
2110  var_component,
2111  elem->point(n),
2112  system.time);
2113  dof_is_fixed[current_dof] = true;
2114  current_dof++;
2115  Gradient grad = g->component(context,
2116  var_component,
2117  elem->point(n),
2118  system.time);
2119  for (unsigned int i=0; i!= dim; ++i)
2120  {
2121  Ue(current_dof) = grad(i);
2122  dof_is_fixed[current_dof] = true;
2123  current_dof++;
2124  }
2125  }
2126  else
2127  libmesh_error();
2128  }
2129 
2130  // In 3D, project any edge values next
2131  if (dim > 2 && cont != DISCONTINUOUS)
2132  {
2133  const std::vector<Point>& xyz_values = edge_fe->get_xyz();
2134  const std::vector<Real>& JxW = edge_fe->get_JxW();
2135 
2136  const std::vector<std::vector<Real> >& phi = edge_fe->get_phi();
2137  const std::vector<std::vector<RealGradient> >* dphi = NULL;
2138  if (cont == C_ONE)
2139  dphi = &(edge_fe->get_dphi());
2140 
2141  for (unsigned int e=0; e != elem->n_edges(); ++e)
2142  {
2143  context.edge = e;
2144  context.edge_fe_reinit();
2145 
2146  const QBase& qedgerule = context.get_edge_qrule();
2147  const unsigned int n_qp = qedgerule.n_points();
2148 
2149  FEInterface::dofs_on_edge(elem, dim, fe_type, e,
2150  side_dofs);
2151 
2152  // Some edge dofs are on nodes and already
2153  // fixed, others are free to calculate
2154  unsigned int free_dofs = 0;
2155  for (unsigned int i=0; i != side_dofs.size(); ++i)
2156  if (!dof_is_fixed[side_dofs[i]])
2157  free_dof[free_dofs++] = i;
2158 
2159  // There may be nothing to project
2160  if (!free_dofs)
2161  continue;
2162 
2163  Ke.resize (free_dofs, free_dofs); Ke.zero();
2164  Fe.resize (free_dofs); Fe.zero();
2165  // The new edge coefficients
2166  DenseVector<Number> Uedge(free_dofs);
2167 
2168  // Loop over the quadrature points
2169  for (unsigned int qp=0; qp<n_qp; qp++)
2170  {
2171  // solution at the quadrature point
2172  Number fineval = f->component(context,
2173  var_component,
2174  xyz_values[qp],
2175  system.time);
2176  // solution grad at the quadrature point
2177  Gradient finegrad;
2178  if (cont == C_ONE)
2179  finegrad = g->component(context,
2180  var_component,
2181  xyz_values[qp],
2182  system.time);
2183 
2184  // Form edge projection matrix
2185  for (unsigned int sidei=0, freei=0;
2186  sidei != side_dofs.size(); ++sidei)
2187  {
2188  unsigned int i = side_dofs[sidei];
2189  // fixed DoFs aren't test functions
2190  if (dof_is_fixed[i])
2191  continue;
2192  for (unsigned int sidej=0, freej=0;
2193  sidej != side_dofs.size(); ++sidej)
2194  {
2195  unsigned int j = side_dofs[sidej];
2196  if (dof_is_fixed[j])
2197  Fe(freei) -= phi[i][qp] * phi[j][qp] *
2198  JxW[qp] * Ue(j);
2199  else
2200  Ke(freei,freej) += phi[i][qp] *
2201  phi[j][qp] * JxW[qp];
2202  if (cont == C_ONE)
2203  {
2204  if (dof_is_fixed[j])
2205  Fe(freei) -= ( (*dphi)[i][qp] *
2206  (*dphi)[j][qp] ) *
2207  JxW[qp] * Ue(j);
2208  else
2209  Ke(freei,freej) += ( (*dphi)[i][qp] *
2210  (*dphi)[j][qp] )
2211  * JxW[qp];
2212  }
2213  if (!dof_is_fixed[j])
2214  freej++;
2215  }
2216  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
2217  if (cont == C_ONE)
2218  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
2219  JxW[qp];
2220  freei++;
2221  }
2222  }
2223 
2224  Ke.cholesky_solve(Fe, Uedge);
2225 
2226  // Transfer new edge solutions to element
2227  for (unsigned int i=0; i != free_dofs; ++i)
2228  {
2229  Number &ui = Ue(side_dofs[free_dof[i]]);
2231  std::abs(ui - Uedge(i)) < TOLERANCE);
2232  ui = Uedge(i);
2233  dof_is_fixed[side_dofs[free_dof[i]]] = true;
2234  }
2235  }
2236  } // end if (dim > 2 && cont != DISCONTINUOUS)
2237 
2238  // Project any side values (edges in 2D, faces in 3D)
2239  if (dim > 1 && cont != DISCONTINUOUS)
2240  {
2241  const std::vector<Point>& xyz_values = side_fe->get_xyz();
2242  const std::vector<Real>& JxW = side_fe->get_JxW();
2243 
2244  const std::vector<std::vector<Real> >& phi = side_fe->get_phi();
2245  const std::vector<std::vector<RealGradient> >* dphi = NULL;
2246  if (cont == C_ONE)
2247  dphi = &(side_fe->get_dphi());
2248 
2249  for (unsigned int s=0; s != elem->n_sides(); ++s)
2250  {
2251  FEInterface::dofs_on_side(elem, dim, fe_type, s,
2252  side_dofs);
2253 
2254  // Some side dofs are on nodes/edges and already
2255  // fixed, others are free to calculate
2256  unsigned int free_dofs = 0;
2257  for (unsigned int i=0; i != side_dofs.size(); ++i)
2258  if (!dof_is_fixed[side_dofs[i]])
2259  free_dof[free_dofs++] = i;
2260 
2261  // There may be nothing to project
2262  if (!free_dofs)
2263  continue;
2264 
2265  Ke.resize (free_dofs, free_dofs); Ke.zero();
2266  Fe.resize (free_dofs); Fe.zero();
2267  // The new side coefficients
2268  DenseVector<Number> Uside(free_dofs);
2269 
2270  context.side = s;
2271  context.side_fe_reinit();
2272 
2273  const QBase& qsiderule = context.get_side_qrule();
2274  const unsigned int n_qp = qsiderule.n_points();
2275 
2276  // Loop over the quadrature points
2277  for (unsigned int qp=0; qp<n_qp; qp++)
2278  {
2279  // solution at the quadrature point
2280  Number fineval = f->component(context,
2281  var_component,
2282  xyz_values[qp],
2283  system.time);
2284  // solution grad at the quadrature point
2285  Gradient finegrad;
2286  if (cont == C_ONE)
2287  finegrad = g->component(context,
2288  var_component,
2289  xyz_values[qp],
2290  system.time);
2291 
2292  // Form side projection matrix
2293  for (unsigned int sidei=0, freei=0;
2294  sidei != side_dofs.size(); ++sidei)
2295  {
2296  unsigned int i = side_dofs[sidei];
2297  // fixed DoFs aren't test functions
2298  if (dof_is_fixed[i])
2299  continue;
2300  for (unsigned int sidej=0, freej=0;
2301  sidej != side_dofs.size(); ++sidej)
2302  {
2303  unsigned int j = side_dofs[sidej];
2304  if (dof_is_fixed[j])
2305  Fe(freei) -= phi[i][qp] * phi[j][qp] *
2306  JxW[qp] * Ue(j);
2307  else
2308  Ke(freei,freej) += phi[i][qp] *
2309  phi[j][qp] * JxW[qp];
2310  if (cont == C_ONE)
2311  {
2312  if (dof_is_fixed[j])
2313  Fe(freei) -= ( (*dphi)[i][qp] *
2314  (*dphi)[j][qp] ) *
2315  JxW[qp] * Ue(j);
2316  else
2317  Ke(freei,freej) += ( (*dphi)[i][qp] *
2318  (*dphi)[j][qp] )
2319  * JxW[qp];
2320  }
2321  if (!dof_is_fixed[j])
2322  freej++;
2323  }
2324  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
2325  if (cont == C_ONE)
2326  Fe(freei) += (finegrad * (*dphi)[i][qp] ) *
2327  JxW[qp];
2328  freei++;
2329  }
2330  }
2331 
2332  Ke.cholesky_solve(Fe, Uside);
2333 
2334  // Transfer new side solutions to element
2335  for (unsigned int i=0; i != free_dofs; ++i)
2336  {
2337  Number &ui = Ue(side_dofs[free_dof[i]]);
2339  std::abs(ui - Uside(i)) < TOLERANCE);
2340  ui = Uside(i);
2341  dof_is_fixed[side_dofs[free_dof[i]]] = true;
2342  }
2343  }
2344  }// end if (dim > 1 && cont != DISCONTINUOUS)
2345 
2346  // Project the interior values, finally
2347 
2348  // Some interior dofs are on nodes/edges/sides and
2349  // already fixed, others are free to calculate
2350  unsigned int free_dofs = 0;
2351  for (unsigned int i=0; i != n_dofs; ++i)
2352  if (!dof_is_fixed[i])
2353  free_dof[free_dofs++] = i;
2354 
2355  // There may be nothing to project
2356  if (free_dofs)
2357  {
2358  context.elem_fe_reinit();
2359 
2360  const QBase& qrule = context.get_element_qrule();
2361  const unsigned int n_qp = qrule.n_points();
2362  const std::vector<Point>& xyz_values = fe->get_xyz();
2363  const std::vector<Real>& JxW = fe->get_JxW();
2364 
2365  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2366  const std::vector<std::vector<RealGradient> >* dphi = NULL;
2367  if (cont == C_ONE)
2368  dphi = &(side_fe->get_dphi());
2369 
2370  Ke.resize (free_dofs, free_dofs); Ke.zero();
2371  Fe.resize (free_dofs); Fe.zero();
2372  // The new interior coefficients
2373  DenseVector<Number> Uint(free_dofs);
2374 
2375  // Loop over the quadrature points
2376  for (unsigned int qp=0; qp<n_qp; qp++)
2377  {
2378  // solution at the quadrature point
2379  Number fineval = f->component(context,
2380  var_component,
2381  xyz_values[qp],
2382  system.time);
2383  // solution grad at the quadrature point
2384  Gradient finegrad;
2385  if (cont == C_ONE)
2386  finegrad = g->component(context,
2387  var_component,
2388  xyz_values[qp],
2389  system.time);
2390 
2391  // Form interior projection matrix
2392  for (unsigned int i=0, freei=0; i != n_dofs; ++i)
2393  {
2394  // fixed DoFs aren't test functions
2395  if (dof_is_fixed[i])
2396  continue;
2397  for (unsigned int j=0, freej=0; j != n_dofs; ++j)
2398  {
2399  if (dof_is_fixed[j])
2400  Fe(freei) -= phi[i][qp] * phi[j][qp] * JxW[qp]
2401  * Ue(j);
2402  else
2403  Ke(freei,freej) += phi[i][qp] * phi[j][qp] *
2404  JxW[qp];
2405  if (cont == C_ONE)
2406  {
2407  if (dof_is_fixed[j])
2408  Fe(freei) -= ( (*dphi)[i][qp] *
2409  (*dphi)[j][qp] ) * JxW[qp] *
2410  Ue(j);
2411  else
2412  Ke(freei,freej) += ( (*dphi)[i][qp] *
2413  (*dphi)[j][qp] ) *
2414  JxW[qp];
2415  }
2416  if (!dof_is_fixed[j])
2417  freej++;
2418  }
2419  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
2420  if (cont == C_ONE)
2421  Fe(freei) += (finegrad * (*dphi)[i][qp] ) * JxW[qp];
2422  freei++;
2423  }
2424  }
2425  Ke.cholesky_solve(Fe, Uint);
2426 
2427  // Transfer new interior solutions to element
2428  for (unsigned int i=0; i != free_dofs; ++i)
2429  {
2430  Number &ui = Ue(free_dof[i]);
2432  std::abs(ui - Uint(i)) < TOLERANCE);
2433  ui = Uint(i);
2434  dof_is_fixed[free_dof[i]] = true;
2435  }
2436 
2437  } // if there are free interior dofs
2438 
2439  // Make sure every DoF got reached!
2440  for (unsigned int i=0; i != n_dofs; ++i)
2441  libmesh_assert(dof_is_fixed[i]);
2442 
2443  const numeric_index_type
2444  first = new_vector.first_local_index(),
2445  last = new_vector.last_local_index();
2446 
2447  // Lock the new_vector since it is shared among threads.
2448  {
2449  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2450 
2451  for (unsigned int i = 0; i < n_dofs; i++)
2452  // We may be projecting a new zero value onto
2453  // an old nonzero approximation - RHS
2454  // if (Ue(i) != 0.)
2455  if ((dof_indices[i] >= first) &&
2456  (dof_indices[i] < last))
2457  {
2458  new_vector.set(dof_indices[i], Ue(i));
2459  }
2460  }
2461  } // end variables loop
2462  } // end elem loop
2463 }
2464 
2465 
2466 
2468 {
2469  // We need data to project
2470  libmesh_assert(f.get());
2471 
2479  // The dimensionality of the current mesh
2480  const unsigned int dim = system.get_mesh().mesh_dimension();
2481 
2482  // The DofMap for this system
2483  const DofMap& dof_map = system.get_dof_map();
2484 
2485  // Boundary info for the current mesh
2486  const BoundaryInfo& boundary_info = *system.get_mesh().boundary_info;
2487 
2488  // The element matrix and RHS for projections.
2489  // Note that Ke is always real-valued, whereas
2490  // Fe may be complex valued if complex number
2491  // support is enabled
2492  DenseMatrix<Real> Ke;
2494  // The new element coefficients
2496 
2497 
2498  // Loop over all the variables we've been requested to project
2499  for (unsigned int v=0; v!=variables.size(); v++)
2500  {
2501  const unsigned int var = variables[v];
2502 
2503  const Variable& variable = dof_map.variable(var);
2504 
2505  const FEType& fe_type = variable.type();
2506 
2507  if (fe_type.family == SCALAR)
2508  continue;
2509 
2510  const unsigned int var_component =
2512 
2513  // Get FE objects of the appropriate type
2514  AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
2515 
2516  // Prepare variables for projection
2517  AutoPtr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
2518  AutoPtr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));
2519 
2520  // The values of the shape functions at the quadrature
2521  // points
2522  const std::vector<std::vector<Real> >& phi = fe->get_phi();
2523 
2524  // The gradients of the shape functions at the quadrature
2525  // points on the child element.
2526  const std::vector<std::vector<RealGradient> > *dphi = NULL;
2527 
2528  const FEContinuity cont = fe->get_continuity();
2529 
2530  if (cont == C_ONE)
2531  {
2532  // We'll need gradient data for a C1 projection
2533  libmesh_assert(g.get());
2534 
2535  const std::vector<std::vector<RealGradient> >&
2536  ref_dphi = fe->get_dphi();
2537  dphi = &ref_dphi;
2538  }
2539 
2540  // The Jacobian * quadrature weight at the quadrature points
2541  const std::vector<Real>& JxW =
2542  fe->get_JxW();
2543 
2544  // The XYZ locations of the quadrature points
2545  const std::vector<Point>& xyz_values =
2546  fe->get_xyz();
2547 
2548  // The global DOF indices
2549  std::vector<dof_id_type> dof_indices;
2550  // Side/edge DOF indices
2551  std::vector<unsigned int> side_dofs;
2552 
2553  // Iterate over all the elements in the range
2554  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
2555  {
2556  const Elem* elem = *elem_it;
2557 
2558  // Per-subdomain variables don't need to be projected on
2559  // elements where they're not active
2560  if (!variable.active_on_subdomain(elem->subdomain_id()))
2561  continue;
2562 
2563  // Find out which nodes, edges and sides are on a requested
2564  // boundary:
2565  std::vector<bool> is_boundary_node(elem->n_nodes(), false),
2566  is_boundary_edge(elem->n_edges(), false),
2567  is_boundary_side(elem->n_sides(), false);
2568  for (unsigned char s=0; s != elem->n_sides(); ++s)
2569  {
2570  // First see if this side has been requested
2571  const std::vector<boundary_id_type>& bc_ids =
2572  boundary_info.boundary_ids (elem, s);
2573  bool do_this_side = false;
2574  for (unsigned int i=0; i != bc_ids.size(); ++i)
2575  if (b.count(bc_ids[i]))
2576  {
2577  do_this_side = true;
2578  break;
2579  }
2580  if (!do_this_side)
2581  continue;
2582 
2583  is_boundary_side[s] = true;
2584 
2585  // Then see what nodes and what edges are on it
2586  for (unsigned int n=0; n != elem->n_nodes(); ++n)
2587  if (elem->is_node_on_side(n,s))
2588  is_boundary_node[n] = true;
2589  for (unsigned int e=0; e != elem->n_edges(); ++e)
2590  if (elem->is_edge_on_side(e,s))
2591  is_boundary_edge[e] = true;
2592  }
2593 
2594  // Update the DOF indices for this element based on
2595  // the current mesh
2596  dof_map.dof_indices (elem, dof_indices, var);
2597 
2598  // The number of DOFs on the element
2599  const unsigned int n_dofs =
2600  libmesh_cast_int<unsigned int>(dof_indices.size());
2601 
2602  // Fixed vs. free DoFs on edge/face projections
2603  std::vector<char> dof_is_fixed(n_dofs, false); // bools
2604  std::vector<int> free_dof(n_dofs, 0);
2605 
2606  // The element type
2607  const ElemType elem_type = elem->type();
2608 
2609  // The number of nodes on the new element
2610  const unsigned int n_nodes = elem->n_nodes();
2611 
2612  // Zero the interpolated values
2613  Ue.resize (n_dofs); Ue.zero();
2614 
2615  // In general, we need a series of
2616  // projections to ensure a unique and continuous
2617  // solution. We start by interpolating boundary nodes, then
2618  // hold those fixed and project boundary edges, then hold
2619  // those fixed and project boundary faces,
2620 
2621  // Interpolate node values first
2622  unsigned int current_dof = 0;
2623  for (unsigned int n=0; n!= n_nodes; ++n)
2624  {
2625  // FIXME: this should go through the DofMap,
2626  // not duplicate dof_indices code badly!
2627  const unsigned int nc =
2628  FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
2629  n);
2630  if (!elem->is_vertex(n) || !is_boundary_node[n])
2631  {
2632  current_dof += nc;
2633  continue;
2634  }
2635  if (cont == DISCONTINUOUS)
2636  {
2637  libmesh_assert_equal_to (nc, 0);
2638  }
2639  // Assume that C_ZERO elements have a single nodal
2640  // value shape function
2641  else if (cont == C_ZERO)
2642  {
2643  libmesh_assert_equal_to (nc, 1);
2644  Ue(current_dof) = f->component(var_component,
2645  elem->point(n),
2646  system.time);
2647  dof_is_fixed[current_dof] = true;
2648  current_dof++;
2649  }
2650  // The hermite element vertex shape functions are weird
2651  else if (fe_type.family == HERMITE)
2652  {
2653  Ue(current_dof) = f->component(var_component,
2654  elem->point(n),
2655  system.time);
2656  dof_is_fixed[current_dof] = true;
2657  current_dof++;
2658  Gradient grad = g->component(var_component,
2659  elem->point(n),
2660  system.time);
2661  // x derivative
2662  Ue(current_dof) = grad(0);
2663  dof_is_fixed[current_dof] = true;
2664  current_dof++;
2665  if (dim > 1)
2666  {
2667  // We'll finite difference mixed derivatives
2668  Point nxminus = elem->point(n),
2669  nxplus = elem->point(n);
2670  nxminus(0) -= TOLERANCE;
2671  nxplus(0) += TOLERANCE;
2672  Gradient gxminus = g->component(var_component,
2673  nxminus,
2674  system.time);
2675  Gradient gxplus = g->component(var_component,
2676  nxplus,
2677  system.time);
2678  // y derivative
2679  Ue(current_dof) = grad(1);
2680  dof_is_fixed[current_dof] = true;
2681  current_dof++;
2682  // xy derivative
2683  Ue(current_dof) = (gxplus(1) - gxminus(1))
2684  / 2. / TOLERANCE;
2685  dof_is_fixed[current_dof] = true;
2686  current_dof++;
2687 
2688  if (dim > 2)
2689  {
2690  // z derivative
2691  Ue(current_dof) = grad(2);
2692  dof_is_fixed[current_dof] = true;
2693  current_dof++;
2694  // xz derivative
2695  Ue(current_dof) = (gxplus(2) - gxminus(2))
2696  / 2. / TOLERANCE;
2697  dof_is_fixed[current_dof] = true;
2698  current_dof++;
2699  // We need new points for yz
2700  Point nyminus = elem->point(n),
2701  nyplus = elem->point(n);
2702  nyminus(1) -= TOLERANCE;
2703  nyplus(1) += TOLERANCE;
2704  Gradient gyminus = g->component(var_component,
2705  nyminus,
2706  system.time);
2707  Gradient gyplus = g->component(var_component,
2708  nyplus,
2709  system.time);
2710  // xz derivative
2711  Ue(current_dof) = (gyplus(2) - gyminus(2))
2712  / 2. / TOLERANCE;
2713  dof_is_fixed[current_dof] = true;
2714  current_dof++;
2715  // Getting a 2nd order xyz is more tedious
2716  Point nxmym = elem->point(n),
2717  nxmyp = elem->point(n),
2718  nxpym = elem->point(n),
2719  nxpyp = elem->point(n);
2720  nxmym(0) -= TOLERANCE;
2721  nxmym(1) -= TOLERANCE;
2722  nxmyp(0) -= TOLERANCE;
2723  nxmyp(1) += TOLERANCE;
2724  nxpym(0) += TOLERANCE;
2725  nxpym(1) -= TOLERANCE;
2726  nxpyp(0) += TOLERANCE;
2727  nxpyp(1) += TOLERANCE;
2728  Gradient gxmym = g->component(var_component,
2729  nxmym,
2730  system.time);
2731  Gradient gxmyp = g->component(var_component,
2732  nxmyp,
2733  system.time);
2734  Gradient gxpym = g->component(var_component,
2735  nxpym,
2736  system.time);
2737  Gradient gxpyp = g->component(var_component,
2738  nxpyp,
2739  system.time);
2740  Number gxzplus = (gxpyp(2) - gxmyp(2))
2741  / 2. / TOLERANCE;
2742  Number gxzminus = (gxpym(2) - gxmym(2))
2743  / 2. / TOLERANCE;
2744  // xyz derivative
2745  Ue(current_dof) = (gxzplus - gxzminus)
2746  / 2. / TOLERANCE;
2747  dof_is_fixed[current_dof] = true;
2748  current_dof++;
2749  }
2750  }
2751  }
2752  // Assume that other C_ONE elements have a single nodal
2753  // value shape function and nodal gradient component
2754  // shape functions
2755  else if (cont == C_ONE)
2756  {
2757  libmesh_assert_equal_to (nc, 1 + dim);
2758  Ue(current_dof) = f->component(var_component,
2759  elem->point(n),
2760  system.time);
2761  dof_is_fixed[current_dof] = true;
2762  current_dof++;
2763  Gradient grad = g->component(var_component,
2764  elem->point(n),
2765  system.time);
2766  for (unsigned int i=0; i!= dim; ++i)
2767  {
2768  Ue(current_dof) = grad(i);
2769  dof_is_fixed[current_dof] = true;
2770  current_dof++;
2771  }
2772  }
2773  else
2774  libmesh_error();
2775  }
2776 
2777  // In 3D, project any edge values next
2778  if (dim > 2 && cont != DISCONTINUOUS)
2779  for (unsigned int e=0; e != elem->n_edges(); ++e)
2780  {
2781  if (!is_boundary_edge[e])
2782  continue;
2783 
2784  FEInterface::dofs_on_edge(elem, dim, fe_type, e,
2785  side_dofs);
2786 
2787  // Some edge dofs are on nodes and already
2788  // fixed, others are free to calculate
2789  unsigned int free_dofs = 0;
2790  for (unsigned int i=0; i != side_dofs.size(); ++i)
2791  if (!dof_is_fixed[side_dofs[i]])
2792  free_dof[free_dofs++] = i;
2793 
2794  // There may be nothing to project
2795  if (!free_dofs)
2796  continue;
2797 
2798  Ke.resize (free_dofs, free_dofs); Ke.zero();
2799  Fe.resize (free_dofs); Fe.zero();
2800  // The new edge coefficients
2801  DenseVector<Number> Uedge(free_dofs);
2802 
2803  // Initialize FE data on the edge
2804  fe->attach_quadrature_rule (qedgerule.get());
2805  fe->edge_reinit (elem, e);
2806  const unsigned int n_qp = qedgerule->n_points();
2807 
2808  // Loop over the quadrature points
2809  for (unsigned int qp=0; qp<n_qp; qp++)
2810  {
2811  // solution at the quadrature point
2812  Number fineval = f->component(var_component,
2813  xyz_values[qp],
2814  system.time);
2815  // solution grad at the quadrature point
2816  Gradient finegrad;
2817  if (cont == C_ONE)
2818  finegrad = g->component(var_component,
2819  xyz_values[qp],
2820  system.time);
2821 
2822  // Form edge projection matrix
2823  for (unsigned int sidei=0, freei=0;
2824  sidei != side_dofs.size(); ++sidei)
2825  {
2826  unsigned int i = side_dofs[sidei];
2827  // fixed DoFs aren't test functions
2828  if (dof_is_fixed[i])
2829  continue;
2830  for (unsigned int sidej=0, freej=0;
2831  sidej != side_dofs.size(); ++sidej)
2832  {
2833  unsigned int j = side_dofs[sidej];
2834  if (dof_is_fixed[j])
2835  Fe(freei) -= phi[i][qp] * phi[j][qp] *
2836  JxW[qp] * Ue(j);
2837  else
2838  Ke(freei,freej) += phi[i][qp] *
2839  phi[j][qp] * JxW[qp];
2840  if (cont == C_ONE)
2841  {
2842  if (dof_is_fixed[j])
2843  Fe(freei) -= ((*dphi)[i][qp] *
2844  (*dphi)[j][qp]) *
2845  JxW[qp] * Ue(j);
2846  else
2847  Ke(freei,freej) += ((*dphi)[i][qp] *
2848  (*dphi)[j][qp])
2849  * JxW[qp];
2850  }
2851  if (!dof_is_fixed[j])
2852  freej++;
2853  }
2854  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
2855  if (cont == C_ONE)
2856  Fe(freei) += (finegrad * (*dphi)[i][qp]) *
2857  JxW[qp];
2858  freei++;
2859  }
2860  }
2861 
2862  Ke.cholesky_solve(Fe, Uedge);
2863 
2864  // Transfer new edge solutions to element
2865  for (unsigned int i=0; i != free_dofs; ++i)
2866  {
2867  Number &ui = Ue(side_dofs[free_dof[i]]);
2869  std::abs(ui - Uedge(i)) < TOLERANCE);
2870  ui = Uedge(i);
2871  dof_is_fixed[side_dofs[free_dof[i]]] = true;
2872  }
2873  }
2874 
2875  // Project any side values (edges in 2D, faces in 3D)
2876  if (dim > 1 && cont != DISCONTINUOUS)
2877  for (unsigned int s=0; s != elem->n_sides(); ++s)
2878  {
2879  if (!is_boundary_side[s])
2880  continue;
2881 
2882  FEInterface::dofs_on_side(elem, dim, fe_type, s,
2883  side_dofs);
2884 
2885  // Some side dofs are on nodes/edges and already
2886  // fixed, others are free to calculate
2887  unsigned int free_dofs = 0;
2888  for (unsigned int i=0; i != side_dofs.size(); ++i)
2889  if (!dof_is_fixed[side_dofs[i]])
2890  free_dof[free_dofs++] = i;
2891 
2892  // There may be nothing to project
2893  if (!free_dofs)
2894  continue;
2895 
2896  Ke.resize (free_dofs, free_dofs); Ke.zero();
2897  Fe.resize (free_dofs); Fe.zero();
2898  // The new side coefficients
2899  DenseVector<Number> Uside(free_dofs);
2900 
2901  // Initialize FE data on the side
2902  fe->attach_quadrature_rule (qsiderule.get());
2903  fe->reinit (elem, s);
2904  const unsigned int n_qp = qsiderule->n_points();
2905 
2906  // Loop over the quadrature points
2907  for (unsigned int qp=0; qp<n_qp; qp++)
2908  {
2909  // solution at the quadrature point
2910  Number fineval = f->component(var_component,
2911  xyz_values[qp],
2912  system.time);
2913  // solution grad at the quadrature point
2914  Gradient finegrad;
2915  if (cont == C_ONE)
2916  finegrad = g->component(var_component,
2917  xyz_values[qp],
2918  system.time);
2919 
2920  // Form side projection matrix
2921  for (unsigned int sidei=0, freei=0;
2922  sidei != side_dofs.size(); ++sidei)
2923  {
2924  unsigned int i = side_dofs[sidei];
2925  // fixed DoFs aren't test functions
2926  if (dof_is_fixed[i])
2927  continue;
2928  for (unsigned int sidej=0, freej=0;
2929  sidej != side_dofs.size(); ++sidej)
2930  {
2931  unsigned int j = side_dofs[sidej];
2932  if (dof_is_fixed[j])
2933  Fe(freei) -= phi[i][qp] * phi[j][qp] *
2934  JxW[qp] * Ue(j);
2935  else
2936  Ke(freei,freej) += phi[i][qp] *
2937  phi[j][qp] * JxW[qp];
2938  if (cont == C_ONE)
2939  {
2940  if (dof_is_fixed[j])
2941  Fe(freei) -= ((*dphi)[i][qp] *
2942  (*dphi)[j][qp]) *
2943  JxW[qp] * Ue(j);
2944  else
2945  Ke(freei,freej) += ((*dphi)[i][qp] *
2946  (*dphi)[j][qp])
2947  * JxW[qp];
2948  }
2949  if (!dof_is_fixed[j])
2950  freej++;
2951  }
2952  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
2953  if (cont == C_ONE)
2954  Fe(freei) += (finegrad * (*dphi)[i][qp]) *
2955  JxW[qp];
2956  freei++;
2957  }
2958  }
2959 
2960  Ke.cholesky_solve(Fe, Uside);
2961 
2962  // Transfer new side solutions to element
2963  for (unsigned int i=0; i != free_dofs; ++i)
2964  {
2965  Number &ui = Ue(side_dofs[free_dof[i]]);
2967  std::abs(ui - Uside(i)) < TOLERANCE);
2968  ui = Uside(i);
2969  dof_is_fixed[side_dofs[free_dof[i]]] = true;
2970  }
2971  }
2972 
2973  const dof_id_type
2974  first = new_vector.first_local_index(),
2975  last = new_vector.last_local_index();
2976 
2977  // Lock the new_vector since it is shared among threads.
2978  {
2979  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2980 
2981  for (unsigned int i = 0; i < n_dofs; i++)
2982  if (dof_is_fixed[i] &&
2983  (dof_indices[i] >= first) &&
2984  (dof_indices[i] < last))
2985  {
2986  new_vector.set(dof_indices[i], Ue(i));
2987  }
2988  }
2989  } // end elem loop
2990  } // end variables loop
2991 }
2992 
2993 
2994 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo