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