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-2025 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
30#ifdef HAVE_CONFIG_H
31#include <oomph-lib-config.h>
32#endif
33
34#include "generic/elements.h"
35#include "generic/shape.h"
37#include "generic/projection.h"
38
39
40namespace 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
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
77
78 /// Access function: Pointer to mass source function (const version)
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
194 {
195 const unsigned n_node = this->nnode();
197 const unsigned n_q_basis = this->nq_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
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
278
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
301
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 {
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
329
330 // Storage for the geometric basis and it's derivatives
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)
340
341 // Get the determinant of the geometric mapping
343
344 // Get the divergence basis (wrt local coords) at local coords s
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 {
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
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
382
383 // Call the 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
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)
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;
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,
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,
452 Vector<double>& norm);
453
454
455 // Z2 stuff:
456
457
458 /// Number off flux terms for Z2 error estimator (use actual flux)
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
475 Shape& psi,
476 Shape& q_basis,
477 Shape& q_test,
478 Shape& p_basis,
479 Shape& p_test,
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,
493 Shape& div_q_test_ds) const = 0;
494
495 // fill in residuals and, if flag==true, jacobian
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(
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(
600 }
601#endif
602 return this->node_pt(0)->ntstorage();
603 }
604
605 /// Number of positional history values
606 /// (Note: count includes current value!)
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(
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
642 double J = this->shape_basis_test_local(s,
643 psi_geom,
644 q_basis,
645 q_test,
646 p_basis,
647 p_test,
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(
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!",
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(
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(
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
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
845
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",
864 }
865
866 // Find the position shape functions
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
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 {
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
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 {
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
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",
988 }
990 }
991
992 // Now assemble residuals
993 if (local_eqn >= 0)
994 {
995 // calculate residuals
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(
1012 }
1013 else
1014 {
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
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 =
1055
1056 // Loop over dofs of field fld
1057 for (unsigned l = 0; l < n_value; l++)
1058 {
1060 if (local_eqn >= 0)
1061 {
1062 // calculate residuals
1065 psi[l] * W;
1066
1067 // Calculate the jacobian
1068 if (flag == 1)
1069 {
1070 for (unsigned l2 = 0; l2 < n_value; l2++)
1071 {
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
1088
1089 // Calculate jacobian and shape functions for that field
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)
1098 this->interpolated_q(s, q_proj);
1099
1100 // Value of the interpolation of element located in base mesh
1102 cast_el_pt->interpolated_q(other_s, q_bar);
1103
1104#ifdef PARANOID
1106 {
1107 std::stringstream error_stream;
1108 error_stream << "Darcy elements are steady!\n";
1109 throw OomphLibError(error_stream.str(),
1112 }
1113#endif
1114
1115 // Loop over dofs of field fld
1116 for (unsigned l = 0; l < n_value; l++)
1117 {
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
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 {
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",
1151 }
1152
1153
1154 break;
1155
1156 default:
1157 throw OomphLibError("Wrong type specificied in Projection_type. "
1158 "This should never happen\n",
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 Vector< Data * > q_edge_data_pt() const =0
Return vector of pointers to the Data objects that store the edge flux values.
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.
virtual Data * q_internal_data_pt() const =0
Return pointer to the Data object that stores the internal flux values.
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.
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.
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 Data * p_data_pt() const =0
Return pointer to the Data object that stores the pressure values.
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 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.
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 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.
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.
SourceFctPt & source_fct_pt()
Access function: Pointer to body force function.
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 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...
FaceGeometry class definition: This policy class is used to allow construction of face elements that ...
Definition elements.h:5002
A general Finite Element class.
Definition elements.h:1317
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:4133
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition elements.h:1967
unsigned nnodal_position_type() const
Return the number of coordinate types that the element requires to interpolate the geometry between t...
Definition elements.h:2467
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition elements.cc:4320
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:3992
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition elements.h:2615
unsigned nnode() const
Return the number of nodes.
Definition elements.h:2214
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition elements.h:1763
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:1512
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
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:1985
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....
Darcy upgraded to become projectable.
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...
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 ...
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
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
SolidFiniteElement class.
Definition elements.h:3565
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).