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