darcy_elements.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifndef OOMPH_RAVIART_THOMAS_DARCY_HEADER
27 #define OOMPH_RAVIART_THOMAS_DARCY_HEADER
28 
29 // Config header generated by autoconfig
30 #ifdef HAVE_CONFIG_H
31 #include <oomph-lib-config.h>
32 #endif
33 
34 #include "../generic/elements.h"
35 #include "../generic/shape.h"
36 #include "../generic/error_estimator.h"
37 #include "../generic/projection.h"
38 
39 
40 namespace oomph
41 {
42  //===========================================================================
43  /// Class implementing the generic maths of the Darcy equations using
44  /// Raviart-Thomas elements with both edge and internal degrees of freedom
45  //===========================================================================
46  template<unsigned DIM>
47  class DarcyEquations : public virtual FiniteElement,
48  public virtual ElementWithZ2ErrorEstimator
49  {
50  public:
51  /// Source function pointer typedef
52  typedef void (*SourceFctPt)(const Vector<double>& x, Vector<double>& f);
53 
54  /// Mass source function pointer typedef
55  typedef void (*MassSourceFctPt)(const Vector<double>& x, double& f);
56 
57  /// Constructor
59 
60  /// Access function: Pointer to body force function
62  {
63  return Source_fct_pt;
64  }
65 
66  /// Access function: Pointer to body force function (const version)
68  {
69  return Source_fct_pt;
70  }
71 
72  /// Access function: Pointer to mass source function
74  {
75  return Mass_source_fct_pt;
76  }
77 
78  /// Access function: Pointer to mass source function (const version)
80  {
81  return Mass_source_fct_pt;
82  }
83 
84  /// Indirect access to the source function - returns 0 if no source
85  /// function has been set
86  void source(const Vector<double>& x, Vector<double>& b) const
87  {
88  // If no function has been set, return zero vector
89  if (Source_fct_pt == 0)
90  {
91  // Get spatial dimension of element
92  unsigned n = dim();
93  for (unsigned i = 0; i < n; i++)
94  {
95  b[i] = 0.0;
96  }
97  }
98  else
99  {
100  // Get body force
101  (*Source_fct_pt)(x, b);
102  }
103  }
104 
105  /// Indirect access to the mass source function - returns 0 if no
106  /// mass source function has been set
107  void mass_source(const Vector<double>& x, double& b) const
108  {
109  // If no function has been set, return zero vector
110  if (Mass_source_fct_pt == 0)
111  {
112  b = 0.0;
113  }
114  else
115  {
116  // Get body force
117  (*Mass_source_fct_pt)(x, b);
118  }
119  }
120 
121  /// Number of values required at node n
122  virtual unsigned required_nvalue(const unsigned& n) const = 0;
123 
124  /// Return the equation number of the n-th edge (flux) degree of freedom
125  virtual int q_edge_local_eqn(const unsigned& n) const = 0;
126 
127  /// Return the equation number of the n-th internal degree of freedom
128  virtual int q_internal_local_eqn(const unsigned& n) const = 0;
129 
130  /// Return vector of pointers to the Data objects that store the
131  /// edge flux values
132  virtual Vector<Data*> q_edge_data_pt() const = 0;
133 
134  /// Return pointer to the Data object that stores the internal flux values
135  virtual Data* q_internal_data_pt() const = 0;
136 
137  /// Return the nodal index at which the nth edge unknown is stored
138  virtual unsigned q_edge_index(const unsigned& n) const = 0;
139 
140  /// Return the index of the internal data where the q_internal
141  /// degrees of freedom are stored
142  virtual unsigned q_internal_index() const = 0;
143 
144  /// Return the number of the node where the nth edge unknown is stored
145  virtual unsigned q_edge_node_number(const unsigned& n) const = 0;
146 
147  /// Return the values of the n-th edge (flux) degree of freedom
148  virtual double q_edge(const unsigned& n) const = 0;
149 
150  /// Return the face index associated with edge flux degree of freedom
152  const unsigned& j) const = 0;
153 
154  /// Return the face index associated with specified edge
155  virtual unsigned face_index_of_edge(const unsigned& j) const = 0;
156 
157  /// Compute the face element coordinates of the nth flux
158  /// interpolation point along an edge
160  const unsigned& edge, const unsigned& n, Vector<double>& s) const = 0;
161 
162  /// Return the values of the internal degree of freedom
163  virtual double q_internal(const unsigned& n) const = 0;
164 
165  /// Set the values of the edge (flux) degree of freedom
166  virtual void set_q_edge(const unsigned& n, const double& value) = 0;
167 
168  /// Set the values of the internal degree of freedom
169  virtual void set_q_internal(const unsigned& n, const double& value) = 0;
170 
171  /// Return the total number of computational basis functions for q
172  unsigned nq_basis() const
173  {
174  return nq_basis_edge() + nq_basis_internal();
175  }
176 
177  /// Return the number of edge basis functions for q
178  virtual unsigned nq_basis_edge() const = 0;
179 
180  /// Return the number of internal basis functions for q
181  virtual unsigned nq_basis_internal() const = 0;
182 
183  /// Returns the local form of the q basis at local coordinate s
184  virtual void get_q_basis_local(const Vector<double>& s,
185  Shape& q_basis) const = 0;
186 
187  /// Returns the local form of the q basis and dbasis/ds at local coordinate
188  /// s
190  Shape& div_q_basis_ds) const = 0;
191 
192  /// Returns the transformed basis at local coordinate s
193  void get_q_basis(const Vector<double>& s, Shape& q_basis) const
194  {
195  const unsigned n_node = this->nnode();
196  Shape psi(n_node, DIM);
197  const unsigned n_q_basis = this->nq_basis();
198  Shape q_basis_local(n_q_basis, DIM);
199  this->get_q_basis_local(s, q_basis_local);
200  (void)this->transform_basis(s, q_basis_local, psi, q_basis);
201  }
202 
203  /// Returns the number of flux interpolation points along each
204  /// edge of the element
205  virtual unsigned nedge_flux_interpolation_point() const = 0;
206 
207  /// Compute the global coordinates of the flux_interpolation
208  /// point associated with the j-th edge-based q basis fct
210  const unsigned& j, Vector<double>& x) const = 0;
211 
212 
213  /// Returns the local coordinate of nth flux interpolation point along an
214  /// edge
216  const unsigned& edge, const unsigned& n) const = 0;
217 
218  /// Returns the global coordinates of the nth flux
219  /// interpolation point along an edge
221  const unsigned& edge, const unsigned& n, Vector<double>& x) const = 0;
222 
223  /// Pin the nth internal q value
224  virtual void pin_q_internal_value(const unsigned& n) = 0;
225 
226  /// Return the equation number of the n-th pressure degree of freedom
227  virtual int p_local_eqn(const unsigned& n) const = 0;
228 
229  /// Return the nth pressure value
230  virtual double p_value(const unsigned& n) const = 0;
231 
232  /// Return the total number of pressure basis functions
233  virtual unsigned np_basis() const = 0;
234 
235  /// Return the pressure basis
236  virtual void get_p_basis(const Vector<double>& s, Shape& p_basis) const = 0;
237 
238  /// Pin the nth pressure value
239  virtual void pin_p_value(const unsigned& n) = 0;
240 
241  /// Set the nth pressure value
242  virtual void set_p_value(const unsigned& n, const double& value) = 0;
243 
244  /// Return pointer to the Data object that stores the pressure values
245  virtual Data* p_data_pt() const = 0;
246 
247  /// Scale the edge basis to allow arbitrary edge mappings
248  virtual void scale_basis(Shape& basis) const = 0;
249 
250  /// Performs a div-conserving transformation of the vector basis
251  /// functions from the reference element to the actual element
252  double transform_basis(const Vector<double>& s,
253  const Shape& q_basis_local,
254  Shape& psi,
255  Shape& q_basis) const;
256 
257  /// Fill in contribution to residuals for the Darcy equations
259  {
261  residuals, GeneralisedElement::Dummy_matrix, 0);
262  }
263 
264  // mjr do finite differences for now
265  // void fill_in_contribution_to_jacobian(Vector<double> &residuals,
266  // DenseMatrix<double> &jacobian)
267  // {
268  // this->fill_in_generic_residual_contribution(residuals,jacobian,1);
269  // }
270 
271  /// Calculate the FE representation of q
273  {
274  unsigned n_q_basis = nq_basis();
275  unsigned n_q_basis_edge = nq_basis_edge();
276 
277  Shape q_basis(n_q_basis, DIM);
278 
279  get_q_basis(s, q_basis);
280  for (unsigned i = 0; i < DIM; i++)
281  {
282  q[i] = 0.0;
283  for (unsigned l = 0; l < n_q_basis_edge; l++)
284  {
285  q[i] += q_edge(l) * q_basis(l, i);
286  }
287  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
288  {
289  q[i] += q_internal(l - n_q_basis_edge) * q_basis(l, i);
290  }
291  }
292  }
293 
294  /// Calculate the FE representation of the i-th component of q
295  double interpolated_q(const Vector<double>& s, const unsigned i) const
296  {
297  unsigned n_q_basis = nq_basis();
298  unsigned n_q_basis_edge = nq_basis_edge();
299 
300  Shape q_basis(n_q_basis, DIM);
301 
302  get_q_basis(s, q_basis);
303  double q_i = 0.0;
304  for (unsigned l = 0; l < n_q_basis_edge; l++)
305  {
306  q_i += q_edge(l) * q_basis(l, i);
307  }
308  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
309  {
310  q_i += q_internal(l - n_q_basis_edge) * q_basis(l, i);
311  }
312 
313  return q_i;
314  }
315 
316  /// Calculate the FE representation of div q
317  void interpolated_div_q(const Vector<double>& s, double& div_q) const
318  {
319  // Zero the divergence
320  div_q = 0;
321 
322  // Get the number of nodes, q basis function, and q edge basis functions
323  unsigned n_node = nnode();
324  const unsigned n_q_basis = nq_basis();
325  const unsigned n_q_basis_edge = nq_basis_edge();
326 
327  // Storage for the divergence basis
328  Shape div_q_basis_ds(n_q_basis);
329 
330  // Storage for the geometric basis and it's derivatives
331  Shape psi(n_node);
332  DShape dpsi(n_node, DIM);
333 
334  // Call the geometric shape functions and their derivatives
335  this->dshape_local(s, psi, dpsi);
336 
337  // Storage for the inverse of the geometric jacobian (just so we can call
338  // the local to eulerian mapping)
339  DenseMatrix<double> inverse_jacobian(DIM);
340 
341  // Get the determinant of the geometric mapping
342  double det = local_to_eulerian_mapping(dpsi, inverse_jacobian);
343 
344  // Get the divergence basis (wrt local coords) at local coords s
345  get_div_q_basis_local(s, div_q_basis_ds);
346 
347  // Add the contribution to the divergence from the edge basis functions
348  for (unsigned l = 0; l < n_q_basis_edge; l++)
349  {
350  div_q += 1.0 / det * div_q_basis_ds(l) * q_edge(l);
351  }
352 
353  // Add the contribution to the divergence from the internal basis
354  // functions
355  for (unsigned l = n_q_basis_edge; l < n_q_basis; l++)
356  {
357  div_q += 1.0 / det * div_q_basis_ds(l) * q_internal(l - n_q_basis_edge);
358  }
359  }
360 
361  /// Calculate the FE representation of div q and return it
363  {
364  // Temporary storage for div q
365  double div_q = 0;
366 
367  // Get the intepolated divergence
368  interpolated_div_q(s, div_q);
369 
370  // Return it
371  return div_q;
372  }
373 
374  /// Calculate the FE representation of p
375  void interpolated_p(const Vector<double>& s, double& p) const
376  {
377  // Get the number of p basis functions
378  unsigned n_p_basis = np_basis();
379 
380  // Storage for the p basis
381  Shape p_basis(n_p_basis);
382 
383  // Call the p basis
384  get_p_basis(s, p_basis);
385 
386  // Zero the pressure
387  p = 0;
388 
389  // Add the contribution to the pressure from each basis function
390  for (unsigned l = 0; l < n_p_basis; l++)
391  {
392  p += p_value(l) * p_basis(l);
393  }
394  }
395 
396  /// Calculate the FE representation of p and return it
397  double interpolated_p(const Vector<double>& s) const
398  {
399  // Temporary storage for p
400  double p = 0;
401 
402  // Get the interpolated pressure
403  interpolated_p(s, p);
404 
405  // Return it
406  return p;
407  }
408 
409 
410  /// Helper function to pin superfluous dofs (empty; can be overloaded
411  /// in projectable elements where we introduce at least one
412  /// dof per node to allow projection during unstructured refinement)
413  virtual void pin_superfluous_darcy_dofs() {}
414 
415 
416  /// Self test -- empty for now.
417  unsigned self_test()
418  {
419  return 0;
420  }
421 
422  /// Output with default number of plot points
423  void output(std::ostream& outfile)
424  {
425  unsigned nplot = 5;
426  output(outfile, nplot);
427  }
428 
429  /// Output FE representation of soln: x,y,q1,q2,div_q,p at
430  /// Nplot^DIM plot points
431  void output(std::ostream& outfile, const unsigned& nplot);
432 
433 
434  /// Output incl. projection of fluxes into direction of
435  /// the specified unit vector
436  void output_with_projected_flux(std::ostream& outfile,
437  const unsigned& nplot,
438  const Vector<double> unit_normal);
439 
440 
441  /// Output FE representation of exact soln: x,y,q1,q2,div_q,p at
442  /// Nplot^DIM plot points
443  void output_fct(std::ostream& outfile,
444  const unsigned& nplot,
446 
447  /// Compute the error between the FE solution and the exact solution
448  /// using the H(div) norm for q and L^2 norm for p
449  void compute_error(std::ostream& outfile,
451  Vector<double>& error,
452  Vector<double>& norm);
453 
454 
455  // Z2 stuff:
456 
457 
458  /// Number off flux terms for Z2 error estimator (use actual flux)
459  unsigned num_Z2_flux_terms()
460  {
461  return DIM;
462  }
463 
464  /// Z2 flux (use actual flux)
466  {
467  interpolated_q(s, flux);
468  }
469 
470 
471  protected:
472  /// Returns the geometric basis, and the q, p and divergence basis
473  /// functions and test functions at local coordinate s
474  virtual double shape_basis_test_local(const Vector<double>& s,
475  Shape& psi,
476  Shape& q_basis,
477  Shape& q_test,
478  Shape& p_basis,
479  Shape& p_test,
480  Shape& div_q_basis_ds,
481  Shape& div_q_test_ds) const = 0;
482 
483  /// Returns the geometric basis, and the q, p and divergence basis
484  /// functions and test functions at integration point ipt
486  const unsigned& ipt,
487  Shape& psi,
488  Shape& q_basis,
489  Shape& q_test,
490  Shape& p_basis,
491  Shape& p_test,
492  Shape& div_q_basis_ds,
493  Shape& div_q_test_ds) const = 0;
494 
495  // fill in residuals and, if flag==true, jacobian
497  Vector<double>& residuals, DenseMatrix<double>& jacobian, bool flag);
498 
499  private:
500  /// Pointer to body force function
502 
503  /// Pointer to the mass source function
505  };
506 
507 
508  /// /////////////////////////////////////////////////////////////////////
509  /// /////////////////////////////////////////////////////////////////////
510  /// /////////////////////////////////////////////////////////////////////
511 
512 
513  //==========================================================
514  /// Darcy upgraded to become projectable
515  //==========================================================
516  template<class DARCY_ELEMENT>
518  : public virtual ProjectableElement<DARCY_ELEMENT>,
519  public virtual ProjectableElementBase
520  {
521  public:
522  /// Constructor [this was only required explicitly
523  /// from gcc 4.5.2 onwards...]
525 
526  /// Specify the values associated with field fld.
527  /// The information is returned in a vector of pairs which comprise
528  /// the Data object and the value within it, that correspond to field fld.
530  {
531 #ifdef PARANOID
532  if (fld > 1)
533  {
534  std::stringstream error_stream;
535  error_stream << "Darcy elements only store 2 fields so fld = " << fld
536  << " is illegal \n";
537  throw OomphLibError(
538  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
539  }
540 #endif
541 
542  // Create the vector
544 
545  // Pressure
546  if (fld == 0)
547  {
548  Data* dat_pt = this->p_data_pt();
549  unsigned nvalue = dat_pt->nvalue();
550  for (unsigned i = 0; i < nvalue; i++)
551  {
552  data_values.push_back(std::make_pair(dat_pt, i));
553  }
554  }
555  else
556  {
557  Vector<Data*> edge_dat_pt = this->q_edge_data_pt();
558  unsigned n = edge_dat_pt.size();
559  for (unsigned j = 0; j < n; j++)
560  {
561  unsigned nvalue = edge_dat_pt[j]->nvalue();
562  for (unsigned i = 0; i < nvalue; i++)
563  {
564  data_values.push_back(std::make_pair(edge_dat_pt[j], i));
565  }
566  }
567  if (this->nq_basis_internal() > 0)
568  {
569  Data* int_dat_pt = this->q_internal_data_pt();
570  unsigned nvalue = int_dat_pt->nvalue();
571  for (unsigned i = 0; i < nvalue; i++)
572  {
573  data_values.push_back(std::make_pair(int_dat_pt, i));
574  }
575  }
576  }
577 
578  // Return the vector
579  return data_values;
580  }
581 
582  /// Number of fields to be projected: 2 (pressure and flux)
584  {
585  return 2;
586  }
587 
588  /// Number of history values to be stored for fld-th field.
589  /// (Note: count includes current value!)
590  unsigned nhistory_values_for_projection(const unsigned& fld)
591  {
592 #ifdef PARANOID
593  if (fld > 1)
594  {
595  std::stringstream error_stream;
596  error_stream << "Darcy elements only store two fields so fld = " << fld
597  << " is illegal\n";
598  throw OomphLibError(
599  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
600  }
601 #endif
602  return this->node_pt(0)->ntstorage();
603  }
604 
605  /// Number of positional history values
606  /// (Note: count includes current value!)
608  {
609  return this->node_pt(0)->position_time_stepper_pt()->ntstorage();
610  }
611 
612  /// Return Jacobian of mapping and shape functions of field fld
613  /// at local coordinate s
614  double jacobian_and_shape_of_field(const unsigned& fld,
615  const Vector<double>& s,
616  Shape& psi)
617  {
618 #ifdef PARANOID
619  if (fld > 1)
620  {
621  std::stringstream error_stream;
622  error_stream << "Darcy elements only store two fields so fld = " << fld
623  << " is illegal.\n";
624  throw OomphLibError(
625  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
626  }
627 #endif
628 
629 
630  // Get the number of geometric nodes, total number of basis functions,
631  // and number of edges basis functions
632  const unsigned n_dim = this->dim();
633  const unsigned n_node = this->nnode();
634  const unsigned n_q_basis = this->nq_basis();
635  const unsigned n_p_basis = this->np_basis();
636 
637  // Storage for the geometric and computational bases and the test
638  // functions
639  Shape psi_geom(n_node), q_basis(n_q_basis, n_dim),
640  q_test(n_q_basis, n_dim), p_basis(n_p_basis), p_test(n_p_basis),
641  div_q_basis_ds(n_q_basis), div_q_test_ds(n_q_basis);
642  double J = this->shape_basis_test_local(s,
643  psi_geom,
644  q_basis,
645  q_test,
646  p_basis,
647  p_test,
648  div_q_basis_ds,
649  div_q_test_ds);
650  // Pressure basis functions
651  if (fld == 0)
652  {
653  unsigned n = p_basis.nindex1();
654  for (unsigned i = 0; i < n; i++)
655  {
656  psi[i] = p_basis[i];
657  }
658  }
659  // Flux basis functions
660  else
661  {
662  unsigned n = q_basis.nindex1();
663  unsigned m = q_basis.nindex2();
664  for (unsigned i = 0; i < n; i++)
665  {
666  for (unsigned j = 0; j < m; j++)
667  {
668  psi(i, j) = q_basis(i, j);
669  }
670  }
671  }
672 
673  return J;
674  }
675 
676 
677  /// Return interpolated field fld at local coordinate s, at time
678  /// level t (t=0: present; t>0: history values)
679  double get_field(const unsigned& t,
680  const unsigned& fld,
681  const Vector<double>& s)
682  {
683 #ifdef PARANOID
684  if (fld > 1)
685  {
686  std::stringstream error_stream;
687  error_stream << "Darcy elements only store two fields so fld = " << fld
688  << " is illegal\n";
689  throw OomphLibError(
690  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
691  }
692 #endif
693 
694  double return_value = 0.0;
695 
696  // Pressure
697  if (fld == 0)
698  {
699  this->interpolated_p(s, return_value);
700  }
701  else
702  {
703  throw OomphLibError("Don't call this function for Darcy!",
704  OOMPH_CURRENT_FUNCTION,
705  OOMPH_EXCEPTION_LOCATION);
706  }
707 
708  return return_value;
709  }
710 
711 
712  /// Return number of values in field fld
713  unsigned nvalue_of_field(const unsigned& fld)
714  {
715 #ifdef PARANOID
716  if (fld > 1)
717  {
718  std::stringstream error_stream;
719  error_stream << "Darcy elements only store two fields so fld = " << fld
720  << " is illegal\n";
721  throw OomphLibError(
722  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
723  }
724 #endif
725 
726  unsigned return_value = 0;
727  if (fld == 0)
728  {
729  return_value = this->np_basis();
730  }
731  else
732  {
733  return_value = this->nq_basis();
734  }
735 
736  return return_value;
737  }
738 
739 
740  /// Return local equation number of value j in field fld.
741  int local_equation(const unsigned& fld, const unsigned& j)
742  {
743 #ifdef PARANOID
744  if (fld > 1)
745  {
746  std::stringstream error_stream;
747  error_stream << "Darcy elements only store two fields so fld = " << fld
748  << " is illegal\n";
749  throw OomphLibError(
750  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
751  }
752 #endif
753 
754  int return_value = 0;
755  // Pressure
756  if (fld == 0)
757  {
758  return_value = this->p_local_eqn(j);
759  }
760  else
761  {
762  unsigned nedge = this->nq_basis_edge();
763  if (j < nedge)
764  {
765  return_value = this->q_edge_local_eqn(j);
766  }
767  else
768  {
769  return_value = this->q_internal_local_eqn(j - nedge);
770  }
771  }
772 
773  return return_value;
774  }
775 
776  /// Output FE representation of soln as in underlying element
777  void output(std::ostream& outfile, const unsigned& nplot)
778  {
779  DARCY_ELEMENT::output(outfile, nplot);
780  }
781 
782  /// Overload required_nvalue to return at least one value
783  unsigned required_nvalue(const unsigned& n) const
784  {
785  return std::max(this->Initial_Nvalue[n], unsigned(1));
786  }
787 
788 
789  /// Helper function to pin superfluous dofs; required because
790  /// we introduce at least one dof per node to allow projection
791  /// during unstructured refinement)
793  {
794  // Pin dofs at vertex nodes (because they're only used for projection)
795  for (unsigned j = 0; j < 3; j++)
796  {
797  this->node_pt(j)->pin(0);
798  }
799  }
800 
801  /// Residual for the projection step. Flag indicates if we
802  /// want the Jacobian (1) or not (0). Virtual so it can be
803  /// overloaded if necessary
805  DenseMatrix<double>& jacobian,
806  const unsigned& flag)
807  {
808  // Am I a solid element
809  SolidFiniteElement* solid_el_pt = dynamic_cast<SolidFiniteElement*>(this);
810 
811  unsigned n_dim = this->dim();
812 
813  // Allocate storage for local coordinates
814  Vector<double> s(n_dim);
815 
816  // Current field
817  unsigned fld = this->Projected_field;
818 
819  // Number of nodes
820  const unsigned n_node = this->nnode();
821  // Number of positional dofs
822  const unsigned n_position_type = this->nnodal_position_type();
823 
824  // Number of dof for current field
825  const unsigned n_value = nvalue_of_field(fld);
826 
827  // Set the value of n_intpt
828  const unsigned n_intpt = this->integral_pt()->nweight();
829 
830  // Loop over the integration points
831  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
832  {
833  // Get the local coordinates of Gauss point
834  for (unsigned i = 0; i < n_dim; i++)
835  s[i] = this->integral_pt()->knot(ipt, i);
836 
837  // Get the integral weight
838  double w = this->integral_pt()->weight(ipt);
839 
840  // Find same point in base mesh using external storage
841  FiniteElement* other_el_pt = 0;
842  other_el_pt = this->external_element_pt(0, ipt);
843  Vector<double> other_s(n_dim);
844  other_s = this->external_element_local_coord(0, ipt);
845 
847  dynamic_cast<ProjectableElement<DARCY_ELEMENT>*>(other_el_pt);
848 
849  // Storage for the local equation and local unknown
850  int local_eqn = 0, local_unknown = 0;
851 
852  // Now set up the three different projection problems
853  switch (Projection_type)
854  {
855  case Lagrangian:
856  {
857  // If we don't have a solid element there is a problem
858  if (solid_el_pt == 0)
859  {
860  throw OomphLibError("Trying to project Lagrangian coordinates in "
861  "non-SolidElement\n",
862  OOMPH_CURRENT_FUNCTION,
863  OOMPH_EXCEPTION_LOCATION);
864  }
865 
866  // Find the position shape functions
867  Shape psi(n_node, n_position_type);
868  // Get the position shape functions
869  this->shape(s, psi);
870  // Get the jacobian
871  double J = this->J_eulerian(s);
872 
873  // Premultiply the weights and the Jacobian
874  double W = w * J;
875 
876  // Get the value of the current position of the 0th coordinate
877  // in the current element
878  // at the current time level (which is the unkown)
879  double interpolated_xi_proj = this->interpolated_x(s, 0);
880 
881  // Get the Lagrangian position in the other element
882  double interpolated_xi_bar =
883  dynamic_cast<SolidFiniteElement*>(cast_el_pt)
884  ->interpolated_xi(other_s, Projected_lagrangian);
885 
886  // Now loop over the nodes and position dofs
887  for (unsigned l = 0; l < n_node; ++l)
888  {
889  // Loop over position unknowns
890  for (unsigned k = 0; k < n_position_type; ++k)
891  {
892  // The local equation is going to be the positional local
893  // equation
894  local_eqn = solid_el_pt->position_local_eqn(l, k, 0);
895 
896  // Now assemble residuals
897  if (local_eqn >= 0)
898  {
899  // calculate residuals
900  residuals[local_eqn] +=
901  (interpolated_xi_proj - interpolated_xi_bar) * psi(l, k) *
902  W;
903 
904  // Calculate the jacobian
905  if (flag == 1)
906  {
907  for (unsigned l2 = 0; l2 < n_node; l2++)
908  {
909  // Loop over position dofs
910  for (unsigned k2 = 0; k2 < n_position_type; k2++)
911  {
912  local_unknown =
913  solid_el_pt->position_local_eqn(l2, k2, 0);
914  if (local_unknown >= 0)
915  {
916  jacobian(local_eqn, local_unknown) +=
917  psi(l2, k2) * psi(l, k) * W;
918  }
919  }
920  }
921  } // end of jacobian
922  }
923  }
924  }
925  } // End of Lagrangian coordinate case
926 
927  break;
928 
929  // Now the coordinate history case
930  case Coordinate:
931  {
932  // Find the position shape functions
933  Shape psi(n_node, n_position_type);
934  // Get the position shape functions
935  this->shape(s, psi);
936  // Get the jacobian
937  double J = this->J_eulerian(s);
938 
939  // Premultiply the weights and the Jacobian
940  double W = w * J;
941 
942  // Get the value of the current position in the current element
943  // at the current time level (which is the unkown)
944  double interpolated_x_proj = 0.0;
945  // If we are a solid element read it out directly from the data
946  if (solid_el_pt != 0)
947  {
948  interpolated_x_proj =
950  }
951  // Otherwise it's the field value at the current time
952  else
953  {
954  interpolated_x_proj = this->get_field(0, fld, s);
955  }
956 
957  // Get the position in the other element
958  double interpolated_x_bar = cast_el_pt->interpolated_x(
960 
961  // Now loop over the nodes and position dofs
962  for (unsigned l = 0; l < n_node; ++l)
963  {
964  // Loop over position unknowns
965  for (unsigned k = 0; k < n_position_type; ++k)
966  {
967  // If I'm a solid element
968  if (solid_el_pt != 0)
969  {
970  // The local equation is going to be the positional local
971  // equation
972  local_eqn =
973  solid_el_pt->position_local_eqn(l, k, Projected_coordinate);
974  }
975  // Otherwise just pick the local unknown of a field to
976  // project into
977  else
978  {
979  // Complain if using generalised position types
980  // but this is not a solid element, because the storage
981  // is then not clear
982  if (n_position_type != 1)
983  {
984  throw OomphLibError("Trying to project generalised "
985  "positions not in SolidElement\n",
986  OOMPH_CURRENT_FUNCTION,
987  OOMPH_EXCEPTION_LOCATION);
988  }
989  local_eqn = local_equation(fld, l);
990  }
991 
992  // Now assemble residuals
993  if (local_eqn >= 0)
994  {
995  // calculate residuals
996  residuals[local_eqn] +=
997  (interpolated_x_proj - interpolated_x_bar) * psi(l, k) * W;
998 
999  // Calculate the jacobian
1000  if (flag == 1)
1001  {
1002  for (unsigned l2 = 0; l2 < n_node; l2++)
1003  {
1004  // Loop over position dofs
1005  for (unsigned k2 = 0; k2 < n_position_type; k2++)
1006  {
1007  // If I'm a solid element
1008  if (solid_el_pt != 0)
1009  {
1010  local_unknown = solid_el_pt->position_local_eqn(
1011  l2, k2, Projected_coordinate);
1012  }
1013  else
1014  {
1015  local_unknown = local_equation(fld, l2);
1016  }
1017 
1018  if (local_unknown >= 0)
1019  {
1020  jacobian(local_eqn, local_unknown) +=
1021  psi(l2, k2) * psi(l, k) * W;
1022  }
1023  }
1024  }
1025  } // end of jacobian
1026  }
1027  }
1028  }
1029  } // End of coordinate case
1030  break;
1031 
1032  // Now the value case
1033  case Value:
1034  {
1035  // Pressure -- "normal" procedure
1036  if (fld == 0)
1037  {
1038  // Field shape function
1039  Shape psi(n_value);
1040 
1041  // Calculate jacobian and shape functions for that field
1042  double J = jacobian_and_shape_of_field(fld, s, psi);
1043 
1044  // Premultiply the weights and the Jacobian
1045  double W = w * J;
1046 
1047  // Value of field in current element at current time level
1048  //(the unknown)
1049  unsigned now = 0;
1050  double interpolated_value_proj = this->get_field(now, fld, s);
1051 
1052  // Value of the interpolation of element located in base mesh
1053  double interpolated_value_bar =
1054  cast_el_pt->get_field(Time_level_for_projection, fld, other_s);
1055 
1056  // Loop over dofs of field fld
1057  for (unsigned l = 0; l < n_value; l++)
1058  {
1059  local_eqn = local_equation(fld, l);
1060  if (local_eqn >= 0)
1061  {
1062  // calculate residuals
1063  residuals[local_eqn] +=
1064  (interpolated_value_proj - interpolated_value_bar) *
1065  psi[l] * W;
1066 
1067  // Calculate the jacobian
1068  if (flag == 1)
1069  {
1070  for (unsigned l2 = 0; l2 < n_value; l2++)
1071  {
1072  local_unknown = local_equation(fld, l2);
1073  if (local_unknown >= 0)
1074  {
1075  jacobian(local_eqn, local_unknown) +=
1076  psi[l2] * psi[l] * W;
1077  }
1078  }
1079  } // end of jacobian
1080  }
1081  }
1082  }
1083  // Flux -- need inner product
1084  else if (fld == 1)
1085  {
1086  // Field shape function
1087  Shape psi(n_value, n_dim);
1088 
1089  // Calculate jacobian and shape functions for that field
1090  double J = jacobian_and_shape_of_field(fld, s, psi);
1091 
1092  // Premultiply the weights and the Jacobian
1093  double W = w * J;
1094 
1095  // Value of flux in current element at current time level
1096  //(the unknown)
1097  Vector<double> q_proj(n_dim);
1098  this->interpolated_q(s, q_proj);
1099 
1100  // Value of the interpolation of element located in base mesh
1101  Vector<double> q_bar(n_dim);
1102  cast_el_pt->interpolated_q(other_s, q_bar);
1103 
1104 #ifdef PARANOID
1105  if (Time_level_for_projection != 0)
1106  {
1107  std::stringstream error_stream;
1108  error_stream << "Darcy elements are steady!\n";
1109  throw OomphLibError(error_stream.str(),
1110  OOMPH_CURRENT_FUNCTION,
1111  OOMPH_EXCEPTION_LOCATION);
1112  }
1113 #endif
1114 
1115  // Loop over dofs of field fld
1116  for (unsigned l = 0; l < n_value; l++)
1117  {
1118  local_eqn = local_equation(fld, l);
1119  if (local_eqn >= 0)
1120  {
1121  // Loop over spatial dimension for dot product
1122  for (unsigned i = 0; i < n_dim; i++)
1123  {
1124  // Add to residuals
1125  residuals[local_eqn] +=
1126  (q_proj[i] - q_bar[i]) * psi(l, i) * W;
1127 
1128  // Calculate the jacobian
1129  if (flag == 1)
1130  {
1131  for (unsigned l2 = 0; l2 < n_value; l2++)
1132  {
1133  local_unknown = local_equation(fld, l2);
1134  if (local_unknown >= 0)
1135  {
1136  jacobian(local_eqn, local_unknown) +=
1137  psi(l2, i) * psi(l, i) * W;
1138  }
1139  }
1140  } // end of jacobian
1141  }
1142  }
1143  }
1144  }
1145  else
1146  {
1147  throw OomphLibError(
1148  "Wrong field specified. This should never happen\n",
1149  OOMPH_CURRENT_FUNCTION,
1150  OOMPH_EXCEPTION_LOCATION);
1151  }
1152 
1153 
1154  break;
1155 
1156  default:
1157  throw OomphLibError("Wrong type specificied in Projection_type. "
1158  "This should never happen\n",
1159  OOMPH_CURRENT_FUNCTION,
1160  OOMPH_EXCEPTION_LOCATION);
1161  }
1162  } // End of the switch statement
1163 
1164  } // End of loop over ipt
1165 
1166  } // End of residual_for_projection function
1167  };
1168 
1169 
1170  //=======================================================================
1171  /// Face geometry for element is the same as that for the underlying
1172  /// wrapped element
1173  //=======================================================================
1174  template<class ELEMENT>
1176  : public virtual FaceGeometry<ELEMENT>
1177  {
1178  public:
1179  FaceGeometry() : FaceGeometry<ELEMENT>() {}
1180  };
1181 
1182 
1183  /// /////////////////////////////////////////////////////////////////////
1184  /// /////////////////////////////////////////////////////////////////////
1185  /// /////////////////////////////////////////////////////////////////////
1186 
1187 
1188 } // namespace oomph
1189 
1190 #endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
Class implementing the generic maths of the Darcy equations using Raviart-Thomas elements with both e...
virtual unsigned required_nvalue(const unsigned &n) const =0
Number of values required at node n.
MassSourceFctPt mass_source_fct_pt() const
Access function: Pointer to mass source function (const version)
MassSourceFctPt Mass_source_fct_pt
Pointer to the mass source function.
virtual unsigned nq_basis_edge() const =0
Return the number of edge basis functions for q.
virtual int q_internal_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th internal degree of freedom.
double interpolated_div_q(const Vector< double > &s)
Calculate the FE representation of div q and return it.
virtual void edge_flux_interpolation_point_global(const unsigned &j, Vector< double > &x) const =0
Compute the global coordinates of the flux_interpolation point associated with the j-th edge-based q ...
virtual unsigned np_basis() const =0
Return the total number of pressure basis functions.
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output FE representation of exact soln: x,y,q1,q2,div_q,p at Nplot^DIM plot points.
virtual void scale_basis(Shape &basis) const =0
Scale the edge basis to allow arbitrary edge mappings.
virtual unsigned q_edge_index(const unsigned &n) const =0
Return the nodal index at which the nth edge unknown is stored.
void output_with_projected_flux(std::ostream &outfile, const unsigned &nplot, const Vector< double > unit_normal)
Output incl. projection of fluxes into direction of the specified unit vector.
double interpolated_q(const Vector< double > &s, const unsigned i) const
Calculate the FE representation of the i-th component of q.
virtual void fill_in_generic_residual_contribution(Vector< double > &residuals, DenseMatrix< double > &jacobian, bool flag)
Fill in residuals and, if flag==true, jacobian.
virtual int p_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th pressure degree of freedom.
virtual double shape_basis_test_local(const Vector< double > &s, Shape &psi, Shape &q_basis, Shape &q_test, Shape &p_basis, Shape &p_test, Shape &div_q_basis_ds, Shape &div_q_test_ds) const =0
Returns the geometric basis, and the q, p and divergence basis functions and test functions at local ...
SourceFctPt Source_fct_pt
Pointer to body force function.
virtual double shape_basis_test_local_at_knot(const unsigned &ipt, Shape &psi, Shape &q_basis, Shape &q_test, Shape &p_basis, Shape &p_test, Shape &div_q_basis_ds, Shape &div_q_test_ds) const =0
Returns the geometric basis, and the q, p and divergence basis functions and test functions at integr...
virtual void edge_flux_interpolation_point_global(const unsigned &edge, const unsigned &n, Vector< double > &x) const =0
Returns the global coordinates of the nth flux interpolation point along an edge.
void get_q_basis(const Vector< double > &s, Shape &q_basis) const
Returns the transformed basis at local coordinate s.
virtual unsigned nedge_flux_interpolation_point() const =0
Returns the number of flux interpolation points along each edge of the element.
virtual int q_edge_local_eqn(const unsigned &n) const =0
Return the equation number of the n-th edge (flux) degree of freedom.
double interpolated_p(const Vector< double > &s) const
Calculate the FE representation of p and return it.
DarcyEquations()
Constructor.
void get_Z2_flux(const Vector< double > &s, Vector< double > &flux)
Z2 flux (use actual flux)
void interpolated_q(const Vector< double > &s, Vector< double > &q) const
Calculate the FE representation of q.
SourceFctPt source_fct_pt() const
Access function: Pointer to body force function (const version)
virtual void set_q_edge(const unsigned &n, const double &value)=0
Set the values of the edge (flux) degree of freedom.
virtual void face_local_coordinate_of_flux_interpolation_point(const unsigned &edge, const unsigned &n, Vector< double > &s) const =0
Compute the face element coordinates of the nth flux interpolation point along an edge.
virtual void pin_q_internal_value(const unsigned &n)=0
Pin the nth internal q value.
virtual unsigned face_index_of_q_edge_basis_fct(const unsigned &j) const =0
Return the face index associated with edge flux degree of freedom.
virtual double p_value(const unsigned &n) const =0
Return the nth pressure value.
virtual void pin_superfluous_darcy_dofs()
Helper function to pin superfluous dofs (empty; can be overloaded in projectable elements where we in...
double transform_basis(const Vector< double > &s, const Shape &q_basis_local, Shape &psi, Shape &q_basis) const
Performs a div-conserving transformation of the vector basis functions from the reference element to ...
virtual void set_p_value(const unsigned &n, const double &value)=0
Set the nth pressure value.
void interpolated_p(const Vector< double > &s, double &p) const
Calculate the FE representation of p.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in contribution to residuals for the Darcy equations.
unsigned num_Z2_flux_terms()
Number off flux terms for Z2 error estimator (use actual flux)
virtual unsigned q_internal_index() const =0
Return the index of the internal data where the q_internal degrees of freedom are stored.
virtual unsigned q_edge_node_number(const unsigned &n) const =0
Return the number of the node where the nth edge unknown is stored.
virtual void pin_p_value(const unsigned &n)=0
Pin the nth pressure value.
virtual Data * q_internal_data_pt() const =0
Return pointer to the Data object that stores the internal flux values.
virtual void set_q_internal(const unsigned &n, const double &value)=0
Set the values of the internal degree of freedom.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, Vector< double > &error, Vector< double > &norm)
Compute the error between the FE solution and the exact solution using the H(div) norm for q and L^2 ...
virtual void get_p_basis(const Vector< double > &s, Shape &p_basis) const =0
Return the pressure basis.
virtual double q_edge(const unsigned &n) const =0
Return the values of the n-th edge (flux) degree of freedom.
void interpolated_div_q(const Vector< double > &s, double &div_q) const
Calculate the FE representation of div q.
void(* MassSourceFctPt)(const Vector< double > &x, double &f)
Mass source function pointer typedef.
unsigned nq_basis() const
Return the total number of computational basis functions for q.
SourceFctPt & source_fct_pt()
Access function: Pointer to body force function.
MassSourceFctPt & mass_source_fct_pt()
Access function: Pointer to mass source function.
void output(std::ostream &outfile)
Output with default number of plot points.
void source(const Vector< double > &x, Vector< double > &b) const
Indirect access to the source function - returns 0 if no source function has been set.
virtual Vector< Data * > q_edge_data_pt() const =0
Return vector of pointers to the Data objects that store the edge flux values.
virtual unsigned nq_basis_internal() const =0
Return the number of internal basis functions for q.
virtual void get_div_q_basis_local(const Vector< double > &s, Shape &div_q_basis_ds) const =0
Returns the local form of the q basis and dbasis/ds at local coordinate s.
void(* SourceFctPt)(const Vector< double > &x, Vector< double > &f)
Source function pointer typedef.
unsigned self_test()
Self test – empty for now.
virtual Data * p_data_pt() const =0
Return pointer to the Data object that stores the pressure values.
void mass_source(const Vector< double > &x, double &b) const
Indirect access to the mass source function - returns 0 if no mass source function has been set.
virtual unsigned face_index_of_edge(const unsigned &j) const =0
Return the face index associated with specified edge.
virtual void get_q_basis_local(const Vector< double > &s, Shape &q_basis) const =0
Returns the local form of the q basis at local coordinate s.
virtual Vector< double > edge_flux_interpolation_point(const unsigned &edge, const unsigned &n) const =0
Returns the local coordinate of nth flux interpolation point along an edge.
virtual double q_internal(const unsigned &n) const =0
Return the values of the internal degree of freedom.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:385
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
unsigned ntstorage() const
Return total number of doubles stored per value to record time history of each value (one for steady ...
Definition: nodes.cc:879
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
Base class for finite elements that can compute the quantities that are required for the Z2 error est...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
Definition: elements.cc:4103
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition: elements.h:2463
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual double local_to_eulerian_mapping(const DShape &dpsids, DenseMatrix< double > &jacobian, DenseMatrix< double > &inverse_jacobian) const
Calculate the mapping from local to Eulerian coordinates, given the derivatives of the shape function...
Definition: elements.h:1508
virtual void dshape_local(const Vector< double > &s, Shape &psi, DShape &dpsids) const
Function to compute the geometric shape functions and derivatives w.r.t. local coordinates at local c...
Definition: elements.h:1981
static DenseMatrix< double > Dummy_matrix
Empty dense matrix used as a dummy argument to combined residual and jacobian functions in the case w...
Definition: elements.h:227
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
An OomphLibError object which should be thrown when an run-time error is encountered....
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Vector< std::pair< Data *, unsigned > > data_values_of_field(const unsigned &fld)
Specify the values associated with field fld. The information is returned in a vector of pairs which ...
double jacobian_and_shape_of_field(const unsigned &fld, const Vector< double > &s, Shape &psi)
Return Jacobian of mapping and shape functions of field fld at local coordinate s.
unsigned nhistory_values_for_coordinate_projection()
Number of positional history values (Note: count includes current value!)
void output(std::ostream &outfile, const unsigned &nplot)
Output FE representation of soln as in underlying element.
double get_field(const unsigned &t, const unsigned &fld, const Vector< double > &s)
Return interpolated field fld at local coordinate s, at time level t (t=0: present; t>0: history valu...
unsigned nvalue_of_field(const unsigned &fld)
Return number of values in field fld.
int local_equation(const unsigned &fld, const unsigned &j)
Return local equation number of value j in field fld.
void pin_superfluous_darcy_dofs()
Helper function to pin superfluous dofs; required because we introduce at least one dof per node to a...
unsigned nhistory_values_for_projection(const unsigned &fld)
Number of history values to be stored for fld-th field. (Note: count includes current value!...
void residual_for_projection(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Residual for the projection step. Flag indicates if we want the Jacobian (1) or not (0)....
ProjectableDarcyElement()
Constructor [this was only required explicitly from gcc 4.5.2 onwards...].
unsigned required_nvalue(const unsigned &n) const
Overload required_nvalue to return at least one value.
unsigned nfields_for_projection()
Number of fields to be projected: 2 (pressure and flux)
Template-free Base class for projectable elements.
Definition: projection.h:55
unsigned Projected_lagrangian
When projecting the Lagrangain coordinates indicate which coordiante is to be projected.
Definition: projection.h:78
Projection_Type Projection_type
Variable to indicate if we're projecting the history values of the nodal coordinates (Coordinate) the...
Definition: projection.h:83
unsigned Time_level_for_projection
Time level we are projecting (0=current values; >0: history values)
Definition: projection.h:70
unsigned Projected_field
Field that is currently being projected.
Definition: projection.h:67
unsigned Projected_coordinate
When projecting the history values of the nodal coordinates, this is the coordinate we're projecting.
Definition: projection.h:74
virtual double get_field(const unsigned &time, const unsigned &fld, const Vector< double > &s)=0
Return the fld-th field at local coordinates s at time-level time (time=0: current value; time>0: his...
Wrapper class for projectable elements. Adds "projectability" to the underlying ELEMENT.
Definition: projection.h:183
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
unsigned nindex1() const
Return the range of index 1 of the shape function object.
Definition: shape.h:259
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: elements.h:3561
int position_local_eqn(const unsigned &n, const unsigned &k, const unsigned &j) const
Access function that returns the local equation number that corresponds to the j-th coordinate of the...
Definition: elements.h:4137
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
void output()
Doc the command line arguments.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...