solid_traction_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// Header file for elements that are used to apply surface loads to
27// the equations of elasticity
28
29#ifndef OOMPH_SOLID_TRACTION_ELEMENTS_HEADER
30#define OOMPH_SOLID_TRACTION_ELEMENTS_HEADER
31
32// Config header generated by autoconfig
33#ifdef HAVE_CONFIG_H
34#include <oomph-lib-config.h>
35#endif
36
37
38// OOMPH-LIB headers
39#include "../generic/Qelements.h"
40#include "../generic/hermite_elements.h"
41
42namespace oomph
43{
44 //=======================================================================
45 /// Namespace containing the zero traction function for solid traction
46 /// elements
47 //=======================================================================
48 namespace SolidTractionElementHelper
49 {
50 //=======================================================================
51 /// Default load function (zero traction)
52 //=======================================================================
54 const Vector<double>& x,
55 const Vector<double>& N,
56 Vector<double>& load)
57 {
58 unsigned n_dim = load.size();
59 for (unsigned i = 0; i < n_dim; i++)
60 {
61 load[i] = 0.0;
62 }
63 }
64
65 } // namespace SolidTractionElementHelper
66
67
68 //======================================================================
69 /// A class for elements that allow the imposition of an applied traction
70 /// in the principle of virtual displacements.
71 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
72 /// class and and thus, we can be generic enough without the need to have
73 /// a separate equations class.
74 //======================================================================
75 template<class ELEMENT>
76 class SolidTractionElement : public virtual FaceGeometry<ELEMENT>,
77 public virtual SolidFaceElement
78 {
79 protected:
80 /// Pointer to an imposed traction function. Arguments:
81 /// Lagrangian coordinate; Eulerian coordinate; outer unit normal;
82 /// applied traction. (Not all of the input arguments will be
83 /// required for all specific load functions but the list should
84 /// cover all cases)
86 const Vector<double>& x,
87 const Vector<double>& n,
88 Vector<double>& result);
89
90
91 /// Get the traction vector: Pass number of integration point
92 /// (dummy), Lagr. coordinate and normal vector and return the load vector
93 /// (not all of the input arguments will be
94 /// required for all specific load functions but the list should
95 /// cover all cases). This function is virtual so it can be
96 /// overloaded for FSI.
97 virtual void get_traction(const unsigned& intpt,
98 const Vector<double>& xi,
99 const Vector<double>& x,
100 const Vector<double>& n,
102 {
103 Traction_fct_pt(xi, x, n, traction);
104 }
105
106
107 /// Helper function that actually calculates the residuals
108 // This small level of indirection is required to avoid calling
109 // fill_in_contribution_to_residuals in fill_in_contribution_to_jacobian
110 // which causes all kinds of pain if overloading later on
112 Vector<double>& residuals);
113
114
115 public:
116 /// Constructor, which takes a "bulk" element and the
117 /// value of the index and its limit
119 const int& face_index,
120 const bool& called_from_refineable_constructor = false)
121 : FaceGeometry<ELEMENT>(), FaceElement()
122 {
123 // Attach the geometrical information to the element. N.B. This function
124 // also assigns nbulk_value from the required_nvalue of the bulk element
125 element_pt->build_face_element(face_index, this);
126
127 // Zero traction
129
130#ifdef PARANOID
131 {
132 // Check that the bulk element is not a refineable 3d element
133 if (!called_from_refineable_constructor)
134 {
135 if (element_pt->dim() == 3)
136 {
137 // Is it refineable
138 RefineableElement* ref_el_pt =
139 dynamic_cast<RefineableElement*>(element_pt);
140 if (ref_el_pt != 0)
141 {
142 if (this->has_hanging_nodes())
143 {
144 throw OomphLibError(
145 "This face element will not work correctly if nodes are "
146 "hanging.\nUse the refineable version instead. ",
147 OOMPH_CURRENT_FUNCTION,
148 OOMPH_EXCEPTION_LOCATION);
149 }
150 }
151 }
152 }
153 }
154#endif
155 }
156
157
158 /// Reference to the traction function pointer
159 void (*&traction_fct_pt())(const Vector<double>& xi,
160 const Vector<double>& x,
161 const Vector<double>& n,
163 {
164 return Traction_fct_pt;
165 }
166
167
168 /// Return the residuals
170 {
172 }
173
174
175 /// Fill in contribution from Jacobian
177 DenseMatrix<double>& jacobian)
178 {
179 // Call the residuals
181
182 // Call the generic FD jacobian calculation
184 jacobian);
185
186 // Derivs w.r.t. to any external data (e.g. during displacement control)
187 this->fill_in_jacobian_from_external_by_fd(residuals, jacobian);
188 }
189
190 /// Output function
191 void output(std::ostream& outfile)
192 {
193 unsigned n_plot = 5;
194 output(outfile, n_plot);
195 }
196
197
198 /// Output function
199 void output(std::ostream& outfile, const unsigned& n_plot)
200 {
201 unsigned n_dim = this->nodal_dimension();
202
203 Vector<double> x(n_dim);
204 Vector<double> xi(n_dim);
205 Vector<double> s(n_dim - 1);
206
207 // Tecplot header info
208 outfile << this->tecplot_zone_string(n_plot);
209
210 // Loop over plot points
211 unsigned num_plot_points = this->nplot_points(n_plot);
212 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
213 {
214 // Get local coordinates of plot point
215 this->get_s_plot(iplot, n_plot, s);
216
217 // Get Eulerian and Lagrangian coordinates
218 this->interpolated_x(s, x);
219 this->interpolated_xi(s, xi);
220
221 // Outer unit normal
222 Vector<double> unit_normal(n_dim);
223 outer_unit_normal(s, unit_normal);
224
225 // Dummy
226 unsigned ipt = 0;
227
228 // Traction vector
230 get_traction(ipt, xi, x, unit_normal, traction);
231
232 // Output the x,y,..
233 for (unsigned i = 0; i < n_dim; i++)
234 {
235 outfile << x[i] << " ";
236 }
237
238 // Output traction
239 for (unsigned i = 0; i < n_dim; i++)
240 {
241 outfile << traction[i] << " ";
242 }
243
244 // Output normal
245 for (unsigned i = 0; i < n_dim; i++)
246 {
247 outfile << unit_normal[i] << " ";
248 }
249
250 outfile << std::endl;
251 }
252
253 // Write tecplot footer (e.g. FE connectivity lists)
254 this->write_tecplot_zone_footer(outfile, n_plot);
255 }
256
257 /// C_style output function
258 void output(FILE* file_pt)
259 {
260 FiniteElement::output(file_pt);
261 }
262
263 /// C-style output function
264 void output(FILE* file_pt, const unsigned& n_plot)
265 {
266 FiniteElement::output(file_pt, n_plot);
267 }
268
269
270 /// Compute traction vector at specified local coordinate
271 /// Should only be used for post-processing; ignores dependence
272 /// on integration point!
274 };
275
276 /// ////////////////////////////////////////////////////////////////////
277 /// ////////////////////////////////////////////////////////////////////
278 /// ////////////////////////////////////////////////////////////////////
279
280 //=====================================================================
281 /// Compute traction vector at specified local coordinate
282 /// Should only be used for post-processing; ignores dependence
283 /// on integration point!
284 //=====================================================================
285 template<class ELEMENT>
287 Vector<double>& traction)
288 {
289 unsigned n_dim = this->nodal_dimension();
290
291 // Position vector
292 Vector<double> x(n_dim);
293 interpolated_x(s, x);
294
295
296 // Lagrangian coordinate
297 Vector<double> xi(n_dim);
298 this->interpolated_xi(s, xi);
299
300 // Outer unit normal
301 Vector<double> unit_normal(n_dim);
302 outer_unit_normal(s, unit_normal);
303
304 // Dummy
305 unsigned ipt = 0;
306
307 // Traction vector
308 get_traction(ipt, xi, x, unit_normal, traction);
309 }
310
311
312 //=====================================================================
313 /// Return the residuals for the SolidTractionElement equations
314 //=====================================================================
315 template<class ELEMENT>
318 {
319 // Find out how many nodes there are
320 unsigned n_node = nnode();
321
322 // Find out how many positional dofs there are
323 unsigned n_position_type = this->nnodal_position_type();
324
325 // Find out the dimension of the node
326 unsigned n_dim = this->nodal_dimension();
327
328 // Integer to hold the local equation number
329 int local_eqn = 0;
330
331 // Set up memory for the shape functions
332 // Note that in this case, the number of lagrangian coordinates is always
333 // equal to the dimension of the nodes
334 Shape psi(n_node, n_position_type);
335 DShape dpsids(n_node, n_position_type, n_dim - 1);
336
337 // Set the value of n_intpt
338 unsigned n_intpt = integral_pt()->nweight();
339
340 // Loop over the integration points
341 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
342 {
343 // Get the integral weight
344 double w = integral_pt()->weight(ipt);
345
346 // Only need to call the local derivatives
347 dshape_local_at_knot(ipt, psi, dpsids);
348
349 // Calculate the Eulerian and Lagrangian coordinates
350 Vector<double> interpolated_x(n_dim, 0.0);
351 Vector<double> interpolated_xi(n_dim, 0.0);
352
353 // Also calculate the surface Vectors (derivatives wrt local coordinates)
354 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
355
356 // Calculate displacements and derivatives
357 for (unsigned l = 0; l < n_node; l++)
358 {
359 // Loop over positional dofs
360 for (unsigned k = 0; k < n_position_type; k++)
361 {
362 // Loop over displacement components (deformed position)
363 for (unsigned i = 0; i < n_dim; i++)
364 {
365 // Calculate the Eulerian and Lagrangian positions
366 interpolated_x[i] +=
367 nodal_position_gen(l, bulk_position_type(k), i) * psi(l, k);
368
369 interpolated_xi[i] +=
370 this->lagrangian_position_gen(l, bulk_position_type(k), i) *
371 psi(l, k);
372
373 // Loop over LOCAL derivative directions, to calculate the
374 // tangent(s)
375 for (unsigned j = 0; j < n_dim - 1; j++)
376 {
377 interpolated_A(j, i) +=
378 nodal_position_gen(l, bulk_position_type(k), i) *
379 dpsids(l, k, j);
380 }
381 }
382 }
383 }
384
385 // Now find the local deformed metric tensor from the tangent Vectors
386 DenseMatrix<double> A(n_dim - 1);
387 for (unsigned i = 0; i < n_dim - 1; i++)
388 {
389 for (unsigned j = 0; j < n_dim - 1; j++)
390 {
391 // Initialise surface metric tensor to zero
392 A(i, j) = 0.0;
393 // Take the dot product
394 for (unsigned k = 0; k < n_dim; k++)
395 {
396 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
397 }
398 }
399 }
400
401 // Get the outer unit normal
402 Vector<double> interpolated_normal(n_dim);
403 outer_unit_normal(ipt, interpolated_normal);
404
405 // Find the determinant of the metric tensor
406 double Adet = 0.0;
407 switch (n_dim)
408 {
409 case 2:
410 Adet = A(0, 0);
411 break;
412 case 3:
413 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
414 break;
415 default:
416 throw OomphLibError(
417 "Wrong dimension in SolidTractionElement",
418 "SolidTractionElement::fill_in_contribution_to_residuals()",
419 OOMPH_EXCEPTION_LOCATION);
420 }
421
422 // Premultiply the weights and the square-root of the determinant of
423 // the metric tensor
424 double W = w * sqrt(Adet);
425
426 // Now calculate the load
427 Vector<double> traction(n_dim);
428 get_traction(
429 ipt, interpolated_xi, interpolated_x, interpolated_normal, traction);
430
431 //=====LOAD TERMS FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
432
433 // Loop over the test functions, nodes of the element
434 for (unsigned l = 0; l < n_node; l++)
435 {
436 // Loop of types of dofs
437 for (unsigned k = 0; k < n_position_type; k++)
438 {
439 // Loop over the displacement components
440 for (unsigned i = 0; i < n_dim; i++)
441 {
442 local_eqn = this->position_local_eqn(l, bulk_position_type(k), i);
443 /*IF it's not a boundary condition*/
444 if (local_eqn >= 0)
445 {
446 // Add the loading terms to the residuals
447 residuals[local_eqn] -= traction[i] * psi(l, k) * W;
448 }
449 }
450 } // End of if not boundary condition
451 } // End of loop over shape functions
452 } // End of loop over integration points
453 }
454
455
456 /// //////////////////////////////////////////////////////////////////////
457 /// //////////////////////////////////////////////////////////////////////
458 /// //////////////////////////////////////////////////////////////////////
459
460
461 //======================================================================
462 /// A class for elements that allow the imposition of an applied traction
463 /// in the principle of virtual displacements.
464 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
465 /// class and and thus, we can be generic enough without the need to have
466 /// a separate equations class.
467 ///
468 /// THIS IS THE REFINEABLE VERSION.
469 //======================================================================
470 template<class ELEMENT>
472 : public virtual SolidTractionElement<ELEMENT>,
474 {
475 public:
476 /// Constructor, which takes a "bulk" element and the face index
478 const int& face_index)
479 : // we're calling this from the constructor of the refineable version.
480 SolidTractionElement<ELEMENT>(element_pt, face_index, true)
481 {
482 }
483
484 /// Destructor should not delete anything
486
487 /// Number of continuously interpolated values are the
488 /// same as those in the bulk element.
490 {
491 return dynamic_cast<ELEMENT*>(this->bulk_element_pt())
493 }
494
495 /// This function returns just the residuals
497 {
498 // Call the generic residuals function
500 }
501
502
503 /// This function returns the residuals and the Jacobian
505 DenseMatrix<double>& jacobian)
506 {
507 // Get the residuals
509
510 // Jacobian only contains derivatives w.r.t. to positions which we
511 // evaluate by finite differencing
513 jacobian);
514
515 // Derivs w.r.t. to any external data (e.g. during displacement control)
516 this->fill_in_jacobian_from_external_by_fd(residuals, jacobian);
517 }
518
519
520 protected:
521 /// This function returns the residuals for the
522 /// traction function.
524 Vector<double>& residuals);
525 };
526
527
528 /// ////////////////////////////////////////////////////////////////////
529 /// ////////////////////////////////////////////////////////////////////
530 /// ////////////////////////////////////////////////////////////////////
531
532
533 //============================================================================
534 /// Function that returns the residuals for the imposed traction solid
535 /// equations
536 //============================================================================
537 template<class ELEMENT>
540 Vector<double>& residuals)
541 {
542 // Find out how many nodes there are
543 unsigned n_node = nnode();
544
545 // Find out how many positional dofs there are
546 unsigned n_position_type = this->nnodal_position_type();
547
548
549#ifdef PARANOID
550 if (n_position_type != 1)
551 {
552 throw OomphLibError(
553 "RefineableSolidTractionElement only works for n_position_type=1",
554 OOMPH_CURRENT_FUNCTION,
555 OOMPH_EXCEPTION_LOCATION);
556 }
557#endif
558
559
560 // Find out the dimension of the node
561 unsigned n_dim = this->nodal_dimension();
562
563 // Integer to hold the local equation number
564 int local_eqn = 0;
565
566 // Set up memory for the shape functions
567 // Note that in this case, the number of lagrangian coordinates is always
568 // equal to the dimension of the nodes
569 Shape psi(n_node, n_position_type);
570 DShape dpsids(n_node, n_position_type, n_dim - 1);
571
572 // Set the value of n_intpt
573 unsigned n_intpt = integral_pt()->nweight();
574
575 // Loop over the integration points
576 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
577 {
578 // Get the integral weight
579 double w = integral_pt()->weight(ipt);
580
581 // Only need to call the local derivatives
582 dshape_local_at_knot(ipt, psi, dpsids);
583
584 // Calculate the Eulerian and Lagrangian coordinates
585 Vector<double> interpolated_x(n_dim, 0.0);
586 Vector<double> interpolated_xi(n_dim, 0.0);
587
588 // Also calculate the surface Vectors (derivatives wrt local coordinates)
589 DenseMatrix<double> interpolated_A(n_dim - 1, n_dim, 0.0);
590
591 // Calculate displacements and derivatives
592 for (unsigned l = 0; l < n_node; l++)
593 {
594 // Loop over positional dofs
595 for (unsigned k = 0; k < n_position_type; k++)
596 {
597 // Loop over displacement components (deformed position)
598 for (unsigned i = 0; i < n_dim; i++)
599 {
600 // Calculate the Eulerian and Lagrangian positions
601 interpolated_x[i] +=
602 nodal_position_gen(l, this->bulk_position_type(k), i) * psi(l, k);
603
604 interpolated_xi[i] +=
605 this->lagrangian_position_gen(l, this->bulk_position_type(k), i) *
606 psi(l, k);
607
608 // Loop over LOCAL derivative directions, to calculate the
609 // tangent(s)
610 for (unsigned j = 0; j < n_dim - 1; j++)
611 {
612 interpolated_A(j, i) +=
613 nodal_position_gen(l, this->bulk_position_type(k), i) *
614 dpsids(l, k, j);
615 }
616 }
617 }
618 }
619
620 // Now find the local deformed metric tensor from the tangent Vectors
621 DenseMatrix<double> A(n_dim - 1);
622 for (unsigned i = 0; i < n_dim - 1; i++)
623 {
624 for (unsigned j = 0; j < n_dim - 1; j++)
625 {
626 // Initialise surface metric tensor to zero
627 A(i, j) = 0.0;
628 // Take the dot product
629 for (unsigned k = 0; k < n_dim; k++)
630 {
631 A(i, j) += interpolated_A(i, k) * interpolated_A(j, k);
632 }
633 }
634 }
635
636 // Get the outer unit normal
637 Vector<double> interpolated_normal(n_dim);
638 this->outer_unit_normal(ipt, interpolated_normal);
639
640 // Find the determinant of the metric tensor
641 double Adet = 0.0;
642 switch (n_dim)
643 {
644 case 2:
645 Adet = A(0, 0);
646 break;
647 case 3:
648 Adet = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
649 break;
650 default:
651 throw OomphLibError(
652 "Wrong dimension in RefineableSolidTractionElement",
653 OOMPH_CURRENT_FUNCTION,
654 OOMPH_EXCEPTION_LOCATION);
655 }
656
657 // Premultiply the weights and the square-root of the determinant of
658 // the metric tensor
659 double W = w * sqrt(Adet);
660
661 // Now calculate the load
662 Vector<double> traction(n_dim);
663 this->get_traction(
664 ipt, interpolated_xi, interpolated_x, interpolated_normal, traction);
665
666 // Number of master nodes and storage for the weight of the shape function
667 unsigned n_master = 1;
668 double hang_weight = 1.0;
669
670 //=====LOAD TERMS FROM PRINCIPLE OF VIRTUAL DISPLACEMENTS========
671
672 // Loop over the test functions, nodes of the element
673 for (unsigned l = 0; l < n_node; l++)
674 {
675 // Get pointer to local node l
676 Node* local_node_pt = node_pt(l);
677
678 // Cache hang status
679 bool is_hanging = local_node_pt->is_hanging();
680
681 // If the node is a hanging node
682 if (is_hanging)
683 {
684 n_master = local_node_pt->hanging_pt()->nmaster();
685 }
686 // Otherwise the node is its own master
687 else
688 {
689 n_master = 1;
690 }
691
692
693 // Storage for local equation numbers at node indexed by
694 // type and direction
695 DenseMatrix<int> position_local_eqn_at_node(n_position_type, n_dim);
696
697 // Loop over the master nodes
698 for (unsigned m = 0; m < n_master; m++)
699 {
700 if (is_hanging)
701 {
702 // Find the equation numbers
703 position_local_eqn_at_node = local_position_hang_eqn(
704 local_node_pt->hanging_pt()->master_node_pt(m));
705
706 // Find the hanging node weight
707 hang_weight = local_node_pt->hanging_pt()->master_weight(m);
708 }
709 else
710 {
711 // Loop of types of dofs
712 for (unsigned k = 0; k < n_position_type; k++)
713 {
714 // Loop over the displacement components
715 for (unsigned i = 0; i < n_dim; i++)
716 {
717 position_local_eqn_at_node(k, i) =
718 this->position_local_eqn(l, k, i);
719 }
720 }
721
722 // Hang weight is one
723 hang_weight = 1.0;
724 }
725
726 // Loop of types of dofs
727 for (unsigned k = 0; k < n_position_type; k++)
728 {
729 // Loop over the displacement components
730 for (unsigned i = 0; i < n_dim; i++)
731 {
732 local_eqn = position_local_eqn_at_node(k, i);
733
734 /*IF it's not a boundary condition*/
735 if (local_eqn >= 0)
736 {
737 /* //Loop over the test functions, nodes of the element */
738 /* for(unsigned l=0;l<n_node;l++) */
739 /* { */
740 /* //Loop of types of dofs */
741 /* for(unsigned k=0;k<n_position_type;k++) */
742 /* { */
743 /* //Loop over the displacement components */
744 /* for(unsigned i=0;i<n_dim;i++) */
745 /* { */
746 /* local_eqn =
747 * this->position_local_eqn(l,bulk_position_type(k),i); */
748 /* /\*IF it's not a boundary condition*\/ */
749 /* if(local_eqn >= 0) */
750 /* { */
751
752
753 // Add the loading terms to the residuals
754 residuals[local_eqn] -=
755 traction[i] * psi(l, k) * W * hang_weight;
756
757 } // End of if not boundary condition
758 }
759 }
760 }
761 }
762 } // End of loop over integration points
763 }
764
765
766 /// //////////////////////////////////////////////////////////////////////
767 /// //////////////////////////////////////////////////////////////////////
768 /// //////////////////////////////////////////////////////////////////////
769
770
771 //=========================================================================
772 /// SolidTractionElement "upgraded" to a FSIWallElement (and thus,
773 /// by inheritance, a GeomObject), so it can be used in FSI.
774 /// The element is templated by the bulk solid element and
775 /// the spatial (Eulerian) dimension of the bulk element.
776 //=========================================================================
777 template<class ELEMENT, unsigned DIM>
778 class FSISolidTractionElement : public virtual SolidTractionElement<ELEMENT>,
779 public virtual FSIWallElement
780 {
781 private:
782 /// Boolean flag to indicate whether the normal is directed into the fluid
784
785 public:
786 /// Constructor: Create element as FSIWallElement (and thus,
787 /// by inheritance, a GeomObject) with DIM-1 Lagrangian and DIM
788 /// Eulerian coordinates. By default, we assume that the
789 /// normal vector computed by the underlying FaceElement
790 /// points into the fluid. If this is not the case, overwrite this
791 /// with the access function
792 /// FSISolidTractionElement::set_normal_pointing_out_of_fluid()
793 /// Constructor for GeomObject is called explicitly because
794 /// of virtual inheritance!
796 FiniteElement* const& element_pt,
797 const int& face_index,
798 const bool& called_from_refineable_constructor = false)
799 : SolidTractionElement<ELEMENT>(
800 element_pt, face_index, called_from_refineable_constructor),
802 {
803 unsigned n_lagr = DIM - 1;
804 unsigned n_dim = DIM;
805 setup_fsi_wall_element(n_lagr, n_dim);
806 }
807
808 /// Destructor: empty
810
811
812 /// Set the normal computed by underlying FaceElement
813 /// to point into the fluid
815 {
817 }
818
819 /// Set the normal computed by underlying FaceElement
820 /// to point out of the fluid
822 {
824 }
825
826 /// Derivative of position vector w.r.t. the SolidFiniteElement's
827 /// Lagrangian coordinates; evaluated at current time.
829 const Vector<double>& s, DenseMatrix<double>& drdxi) const
830 {
831 throw OomphLibError("Broken -- who calls this? \n",
832 OOMPH_CURRENT_FUNCTION,
833 OOMPH_EXCEPTION_LOCATION);
834 }
835
836
837 /// Final overload... Forwards to the version in the FSIWallElement
839 DenseMatrix<double>& jacobian)
840 {
841 // Call the underlying element's jacobian function
843 jacobian);
844 // Add the contribution of the external load data by finite differences
846 }
847
848
849 /// Get the load vector: Pass number of the integration point,
850 /// Lagr. coordinate, Eulerian coordinate (neither of the last two
851 /// are used in the FSI implementation of this function!) and normal vector
852 /// and return the load vector, taking
853 /// the sign of the normal into account.
854 virtual void get_traction(const unsigned& intpt,
855 const Vector<double>& xi,
856 const Vector<double>& x,
857 const Vector<double>& n,
859 {
860 // Get the fluid load on the wall stress scale, i.e. this
861 // includes the ratio of stresses represented by Q.
862 fluid_load_vector(intpt, n, traction);
863
864 // If the normal is outer to the fluid switch the direction
866 {
867 for (unsigned i = 0; i < DIM; i++)
868 {
869 traction[i] *= -1.0;
870 }
871 }
872
873 } // end of get_traction
874
875
876 /// Output function: Note we can only output the traction
877 /// at Gauss points so n_plot is actually ignored.
878 void output(std::ostream& outfile, const unsigned& n_plot)
879 {
880 // Tecplot header info
881 outfile << "ZONE" << std::endl;
882
883 // Find the number of Gauss points of the element
884 unsigned n_intpt = this->integral_pt()->nweight();
885
886 // Find the dimension of the element (i.e. its number of local
887 // coordinates)
888 unsigned el_dim = this->dim();
889
890 // Set storage for the local coordinates of the Gauss points
891 // in the solid
892 Vector<double> s(el_dim);
893
894 // Loop over the integration points
895 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
896 {
897 // Loop over the dimension of the solid element and find the local
898 // coordinates of the Gauss points, then the global
899 for (unsigned i = 0; i < el_dim; i++)
900 {
901 s[i] = integral_pt()->knot(ipt, i);
902 }
903
904 // Eulerian position of Gauss point
905 Vector<double> r(DIM);
906 interpolated_x(s, r);
907
908 // Outer unit normal
909 Vector<double> unit_normal(DIM);
910 this->outer_unit_normal(s, unit_normal);
911
912 // Compute traction with dummy arguments for xi
913 Vector<double> xi(el_dim, 0.0);
915 // Don't get the traction if you are a halo element because
916 // the external elements will not have been set up
917#ifdef OOMPH_HAS_MPI
918 if (!this->is_halo())
919#endif
920 {
921 get_traction(ipt, xi, r, unit_normal, traction);
922 }
923
924
925 for (unsigned i = 0; i < DIM; i++)
926 {
927 outfile << r[i] << " ";
928 }
929 for (unsigned i = 0; i < DIM; i++)
930 {
931 outfile << traction[i] << " ";
932 }
933 for (unsigned i = 0; i < DIM; i++)
934 {
935 outfile << unit_normal[i] << " ";
936 }
937 outfile << std::endl;
938 }
939 }
940
941 /// Broken overloaded reference to the traction function pointer.
942 /// It doesn't make sense to specify an external
943 /// traction.
944 virtual void (*&traction_fct_pt())(const Vector<double>& xi,
945 const Vector<double>& x,
946 const Vector<double>& n,
948 {
949 throw OomphLibError("It doesn't make sense to specify an external "
950 "traction in an FSI context",
951 OOMPH_CURRENT_FUNCTION,
952 OOMPH_EXCEPTION_LOCATION);
953
954 // Dummy return to shut up the compiler
955 return this->Traction_fct_pt;
956 }
957
958 /// The number of "DOF types" that degrees of freedom in this element
959 /// are sub-divided into: Just the solid degrees of freedom themselves.
960 unsigned ndof_types() const
961 {
962 return 1;
963 }
964
965 /// Create a list of pairs for all unknowns in this element,
966 /// so that the first entry in each pair contains the global equation
967 /// number of the unknown, while the second one contains the number
968 /// of the "DOF types" that this unknown is associated with.
969 /// (Function can obviously only be called if the equation numbering
970 /// scheme has been set up.)
972 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const;
973 };
974
975
976 //=============================================================================
977 /// Create a list of pairs for all unknowns in this element,
978 /// so that the first entry in each pair contains the global equation
979 /// number of the unknown, while the second one contains the number
980 /// of the "DOF types" that this unknown is associated with.
981 /// (Function can obviously only be called if the equation numbering
982 /// scheme has been set up.) This element is only in charge of the solid dofs.
983 //=============================================================================
984 template<class ELEMENT, unsigned DIM>
986 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
987 {
988 // temporary pair (used to store dof lookup prior to being added to list)
989 std::pair<unsigned, unsigned> dof_lookup;
990
991 // number of nodes
992 const unsigned n_node = this->nnode();
993
994 // Get the number of position dofs and dimensions at the node
995 const unsigned n_position_type = nnodal_position_type();
996 const unsigned nodal_dim = nodal_dimension();
997
998 // Integer storage for local unknown
999 int local_unknown = 0;
1000
1001 // Loop over the nodes
1002 for (unsigned n = 0; n < n_node; n++)
1003 {
1004 // Loop over position dofs
1005 for (unsigned k = 0; k < n_position_type; k++)
1006 {
1007 // Loop over dimension
1008 for (unsigned i = 0; i < nodal_dim; i++)
1009 {
1010 // If the variable is free
1011 local_unknown = position_local_eqn(n, k, i);
1012
1013 // ignore pinned values
1014 if (local_unknown >= 0)
1015 {
1016 // store dof lookup in temporary pair: First entry in pair
1017 // is global equation number; second entry is dof type
1018 dof_lookup.first = this->eqn_number(local_unknown);
1019 dof_lookup.second = 0;
1020
1021 // add to list
1022 dof_lookup_list.push_front(dof_lookup);
1023 }
1024 }
1025 }
1026 }
1027 }
1028
1029
1030 /// //////////////////////////////////////////////////////////////////////
1031 /// //////////////////////////////////////////////////////////////////////
1032 /// //////////////////////////////////////////////////////////////////////
1033
1034
1035 //=========================================================================
1036 /// RefineableSolidTractionElement "upgraded" to a FSIWallElement (and thus,
1037 /// by inheritance, a GeomObject), so it can be used in FSI.
1038 /// The element is templated by the bulk solid element and
1039 /// the spatial (Eulerian) dimension of the bulk element.
1040 //=========================================================================
1041 template<class ELEMENT, unsigned DIM>
1043 : public virtual RefineableSolidTractionElement<ELEMENT>,
1044 public virtual FSISolidTractionElement<ELEMENT, DIM>,
1045 public virtual FSIWallElement
1046 {
1047 public:
1048 /// Constructor: Create element as FSIWallElement (and thus,
1049 /// by inheritance, a GeomObject) with DIM-1 Lagrangian and DIM
1050 /// Eulerian coordinates. By default, we assume that the
1051 /// normal vector computed by the underlying FaceElement
1052 /// points into the fluid. If this is not the case, overwrite this
1053 /// with the access function
1054 /// FSISolidTractionElement::set_normal_pointing_out_of_fluid()
1055 /// Constructor for GeomObject is called explicitly because
1056 /// of virtual inheritance!
1058 const int& face_index)
1059 : SolidTractionElement<ELEMENT>(element_pt, face_index, true),
1060 RefineableSolidTractionElement<ELEMENT>(element_pt, face_index),
1061 FSISolidTractionElement<ELEMENT, DIM>(element_pt, face_index, true)
1062 {
1063 }
1064
1065 /// Destructor: empty
1067
1068
1069 /// Final overload. Get contributions from refineable solid
1070 /// traction element and derivatives from external data
1072 DenseMatrix<double>& jacobian)
1073 {
1074 // Call the underlying element's jacobian function
1076 residuals, jacobian);
1077 // Add the contribution of the external load data by finite differences
1079 }
1080 };
1081
1082
1083 /// /////////////////////////////////////////////////////////////////////////
1084 /// /////////////////////////////////////////////////////////////////////////
1085 /// /////////////////////////////////////////////////////////////////////////
1086
1087
1088 //======================================================================
1089 /// A class for elements that allow the imposition of a displacement
1090 /// constraint for "bulk" solid elements via a Lagrange multiplier.
1091 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
1092 /// class and and thus, we can be generic enough without the need to have
1093 /// a separate equations class.
1094 /// \b NOTE: Currently (and for the foreseeable future) this
1095 /// element only works with bulk elements that do not have
1096 /// generalised degrees of freedom (so it won't work with
1097 /// Hermite-type elements, say). The additional functionality
1098 /// to deal with such elements could easily be added (once a
1099 /// a suitable test case is written). For now we simply throw
1100 /// errors if an attempt is made to use the element with an unsuitable
1101 /// bulk element.
1102 //======================================================================
1103 template<class ELEMENT>
1105 : public virtual FaceGeometry<ELEMENT>,
1106 public virtual SolidFaceElement
1107 {
1108 public:
1109 /// Constructor takes a "bulk" element and the
1110 /// index that identifies which face the FaceElement is supposed
1111 /// to be attached to. The optional identifier can be used
1112 /// to distinguish the additional nodal values created by
1113 /// this element from thos created by other FaceElements.
1115 FiniteElement* const& element_pt,
1116 const int& face_index,
1117 const unsigned& id = 0,
1118 const bool& called_from_refineable_constructor = false)
1120 {
1121 // Store the ID of the FaceElement -- this is used to distinguish
1122 // it from any others
1123 Id = id;
1124
1125 // By default sparsify, i.e. check if the GeomObject that
1126 // defines the boundary contains sub-GeomObjects. If so,
1127 // only use their GeomData as the external Data that affects
1128 // this element's residuals.
1129 Sparsify = true;
1130
1131 // Build the face element
1132 element_pt->build_face_element(face_index, this);
1133
1134#ifdef PARANOID
1135 {
1136 // Initialise number of assigned geom Data.
1138
1139 // Check that the bulk element is not a refineable 3d element
1140 if (!called_from_refineable_constructor)
1141 {
1142 if (element_pt->dim() == 3)
1143 {
1144 // Is it refineable
1145 RefineableElement* ref_el_pt =
1146 dynamic_cast<RefineableElement*>(element_pt);
1147 if (ref_el_pt != 0)
1148 {
1149 if (this->has_hanging_nodes())
1150 {
1151 throw OomphLibError(
1152 "This face element will not work correctly if nodes are "
1153 "hanging\nUse the refineable version instead. ",
1154 OOMPH_CURRENT_FUNCTION,
1155 OOMPH_EXCEPTION_LOCATION);
1156 }
1157 }
1158 }
1159 }
1160 }
1161
1162 {
1163 // Check that the bulk element does not require generalised positional
1164 // degrees of freedom
1165 if (element_pt->nnodal_position_type() != 1)
1166 {
1167 throw OomphLibError(
1168 "ImposeDisplacementByLagrangeMultiplierElement cannot (currently) "
1169 "be used with elements that have generalised positional dofs",
1170 OOMPH_CURRENT_FUNCTION,
1171 OOMPH_EXCEPTION_LOCATION);
1172 }
1173 }
1174#endif
1175
1176 // Dimension of the bulk element
1177 unsigned dim = element_pt->dim();
1178
1179 // We need dim additional values for each FaceElement node
1180 // to store the dim Lagrange multipliers.
1181 Vector<unsigned> n_additional_values(nnode(), dim);
1182
1183 // Now add storage for Lagrange multipliers and set the map containing
1184 // the position of the first entry of this face element's
1185 // additional values.
1186 add_additional_values(n_additional_values, id);
1187 }
1188
1189
1190 /// Access to GeomObject that specifies the prescribed
1191 /// boundary displacement; GeomObject is assumed to be
1192 /// parametrised by the same coordinate that is used as
1193 /// the boundary coordinate in the bulk solid mesh to which
1194 /// this element is attached.
1196 {
1198 }
1199
1200
1201 /// Set GeomObject that specifies the prescribed
1202 /// boundary displacement; GeomObject is assumed to be
1203 /// parametrised by the same coordinate that is used as
1204 /// the boundary coordinate in the bulk solid mesh to which
1205 /// this element is attached. GeomData of GeomObject
1206 /// is added to this element's external Data. Also specify
1207 /// the boundary number in the bulk mesh to which this element is
1208 /// attached.
1211 const unsigned& boundary_number_in_bulk_mesh)
1212 {
1213 // Record boundary number
1214#ifdef PARANOID
1216#endif
1218
1219
1220 // Store (possibly compound) GeomObject that specifies the
1221 // the desired boundary shape.
1223
1224
1225 // Don't sparsify: Use all the geometric Data associated with
1226 // the (possibly compound) GeomObject that specifies the
1227 // boundary shape as external data for this element.
1228 if (!Sparsify)
1229 {
1230 unsigned n_geom_data = boundary_shape_geom_object_pt->ngeom_data();
1231
1232#ifdef PARANOID
1233 if ((this->nexternal_data() > 0) &&
1235 {
1236 std::ostringstream error_message;
1237 error_message << "About to wipe external data for "
1238 << "ImposeDisplacementByLagrangeMultiplierElement.\n"
1239 << "I noted that N_assigned_geom_data = "
1241 << " != nexternal_data() = " << this->nexternal_data()
1242 << " \n"
1243 << "so we're going to wipe some external data that\n"
1244 << "is not geometric Data of the GeomObject that\n"
1245 << "specifies the desired boundary shape.\n"
1246 << std::endl;
1247 throw OomphLibError(error_message.str(),
1248 OOMPH_CURRENT_FUNCTION,
1249 OOMPH_EXCEPTION_LOCATION);
1250 }
1251#endif
1252 this->flush_external_data();
1253 for (unsigned i = 0; i < n_geom_data; i++)
1254 {
1256 }
1257#ifdef PARANOID
1258 N_assigned_geom_data = n_geom_data;
1259#endif
1260 }
1261 // Sparsify: Use locate_zeta to determine the sub-GeomObjects that
1262 // make up the (possibly compound) GeomObject that specifies the
1263 // boundary shape. Use their geometric Data as external data for
1264 // this element.
1265 else
1266 {
1267 // Find out how many nodes there are
1268 unsigned n_node = nnode();
1269
1270 // Get the number of position dofs and dimensions at the node
1271 const unsigned n_position_type = nnodal_position_type();
1272
1273 // Dimension of element
1274 unsigned dim_el = dim();
1275
1276 // Set up memory for the shape functions
1277 Shape psi(n_node);
1278
1279
1280#ifdef PARANOID
1281 if ((this->nexternal_data() > 0) &&
1283 {
1284 std::ostringstream error_message;
1285 error_message << "About to wipe external data for "
1286 << "ImposeDisplacementByLagrangeMultiplierElement.\n"
1287 << "I noted that N_assigned_geom_data = "
1289 << " != nexternal_data() = " << this->nexternal_data()
1290 << " \n"
1291 << "so we're going to wipe some external data that\n"
1292 << "is not geometric Data of the GeomObject that\n"
1293 << "specifies the desired boundary shape.\n"
1294 << std::endl;
1295 throw OomphLibError(error_message.str(),
1296 OOMPH_CURRENT_FUNCTION,
1297 OOMPH_EXCEPTION_LOCATION);
1298 }
1299#endif
1300
1301 // Flush the data
1302 this->flush_external_data();
1303
1304#ifdef PARANOID
1306#endif
1307
1308 // Prepare local storage
1309 unsigned n_intpt = integral_pt()->nweight();
1310 Sub_geom_object_pt.resize(n_intpt);
1311 Zeta_sub_geom_object.resize(n_intpt);
1312
1313 // Loop over the integration points
1314 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1315 {
1316 // Get shape function
1317 shape_at_knot(ipt, psi);
1318
1319 // Calculate the intrinsic coordinates
1320 Vector<double> zeta(dim_el, 0.0);
1321 Vector<double> s(dim_el);
1322
1323 // Loop over nodes
1324 for (unsigned j = 0; j < n_node; j++)
1325 {
1326 for (unsigned k = 0; k < n_position_type; k++)
1327 {
1328 // Assemble the intrinsic coordinate
1329 for (unsigned i = 0; i < dim_el; i++)
1330 {
1331 zeta[i] += zeta_nodal(j, k, i) * psi(j, k);
1332 }
1333 }
1334 }
1335
1336 // Find sub-GeomObject and local coordinate within it
1337 // at integration point
1338 Zeta_sub_geom_object[ipt].resize(dim_el);
1340 zeta, Sub_geom_object_pt[ipt], Zeta_sub_geom_object[ipt]);
1341
1342 unsigned n_geom_data = Sub_geom_object_pt[ipt]->ngeom_data();
1343 for (unsigned i = 0; i < n_geom_data; i++)
1344 {
1346 }
1347#ifdef PARANOID
1348 N_assigned_geom_data += n_geom_data;
1349#endif
1350 }
1351 }
1352 }
1353
1354 /// Fill in the residuals
1356 {
1357 // Call the generic routine with the flag set to 0
1359 residuals, GeneralisedElement::Dummy_matrix, 0);
1360 }
1361
1362
1363 /// Fill in contribution from Jacobian
1365 DenseMatrix<double>& jacobian)
1366 {
1367 // Call the generic routine with the flag set to 1
1369 residuals, jacobian, 1);
1370
1371 // Fill in the derivatives w.r.t. external data by FD, re-using
1372 // the pre-computed residual vector
1373 fill_in_jacobian_from_external_by_fd(residuals, jacobian);
1374 }
1375
1376
1377 /// Fill in contribution to Mass matrix and
1378 /// Jacobian. There is no contributiont to mass matrix
1379 /// so simply call the fill_in_contribution_to_jacobian term
1380 /// Note that the Jacobian is multiplied by minus one to
1381 /// ensure that the mass matrix is positive semi-definite.
1383 Vector<double>& residuals,
1384 DenseMatrix<double>& jacobian,
1385 DenseMatrix<double>& mass_matrix)
1386 {
1387 // Just call the jacobian calculation
1388 fill_in_contribution_to_jacobian(residuals, jacobian);
1389
1390 // Multiply the residuals and jacobian by minus one
1391 const unsigned n_dof = this->ndof();
1392 for (unsigned i = 0; i < n_dof; i++)
1393 {
1394 residuals[i] *= -1.0;
1395 for (unsigned j = 0; j < n_dof; j++)
1396 {
1397 jacobian(i, j) *= -1.0;
1398 }
1399 }
1400 }
1401
1402
1403 /// Output function
1404 void output(std::ostream& outfile, const unsigned& n_plot)
1405 {
1406 // Elemental dimension
1407 unsigned dim_el = dim();
1408
1409 // Find the number of positional types
1410 unsigned n_position_type = this->nnodal_position_type();
1411
1412#ifdef PARANOID
1413 if (n_position_type != 1)
1414 {
1415 throw OomphLibError(
1416 "ImposeDisplacementByLagrangeMultiplierElement cannot (currently) be "
1417 "used with elements that have generalised positional dofs",
1418 OOMPH_CURRENT_FUNCTION,
1419 OOMPH_EXCEPTION_LOCATION);
1420 }
1421#endif
1422
1423
1424 // Local coord
1425 Vector<double> s(dim_el);
1426
1427 // # of nodes,
1428 unsigned n_node = nnode();
1429 Shape psi(n_node, n_position_type);
1430
1431 // Tecplot header info
1432 outfile << this->tecplot_zone_string(n_plot);
1433
1434 // Loop over plot points
1435 unsigned num_plot_points = this->nplot_points(n_plot);
1436 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
1437 {
1438 // Get local coordinates of plot point
1439 this->get_s_plot(iplot, n_plot, s);
1440
1441 // Get shape function
1442 shape(s, psi);
1443
1444 // Calculate the Eulerian coordinates and Lagrange multiplier
1445 Vector<double> x(dim_el + 1, 0.0);
1446 Vector<double> lambda(dim_el + 1, 0.0);
1447 Vector<double> zeta(dim_el, 0.0);
1448 for (unsigned j = 0; j < n_node; j++)
1449 {
1450 // Cast to a boundary node
1451 BoundaryNodeBase* bnod_pt =
1452 dynamic_cast<BoundaryNodeBase*>(node_pt(j));
1453
1454 // get the node pt
1455 Node* nod_pt = node_pt(j);
1456
1457 // Get the index of the first nodal value associated with
1458 // this FaceElement
1459 unsigned first_index =
1461
1462 // higher dimensional quantities
1463 for (unsigned i = 0; i < dim_el + 1; i++)
1464 {
1465 x[i] += nodal_position(j, i) * psi(j, 0); // need to sort
1466 // this out properly
1467 // for generalised dofs
1468 lambda[i] += nod_pt->value(first_index + i) * psi(j, 0);
1469 }
1470 // In-element quantities
1471 for (unsigned i = 0; i < dim_el; i++)
1472 {
1473 // Loop over positional types
1474 for (unsigned k = 0; k < n_position_type; k++)
1475 {
1476 zeta[i] += zeta_nodal(j, k, i) * psi(j, k);
1477 }
1478 }
1479 }
1480
1481 // Get prescribed wall shape
1482 Vector<double> r_prescribed(dim_el + 1);
1483 Boundary_shape_geom_object_pt->position(zeta, r_prescribed);
1484
1485 // Output stuff
1486 for (unsigned i = 0; i < dim_el + 1; i++)
1487 {
1488 outfile << x[i] << " ";
1489 }
1490 for (unsigned i = 0; i < dim_el + 1; i++)
1491 {
1492 outfile << -lambda[i] << " ";
1493 }
1494 for (unsigned i = 0; i < dim_el + 1; i++)
1495 {
1496 outfile << r_prescribed[i] << " ";
1497 }
1498 /* for(unsigned i=0;i<dim_el;i++) */
1499 /* { */
1500 /* outfile << zeta[i] << " "; */
1501 /* } */
1502 outfile << std::endl;
1503 }
1504 }
1505
1506
1507 /// Output function
1508 void output(std::ostream& outfile)
1509 {
1510 unsigned n_plot = 5;
1511 output(outfile, n_plot);
1512 }
1513
1514
1515 /// Compute square of L2 norm of error between
1516 /// prescribed and actual displacement
1518 {
1519 // Find out how many positional dofs there are
1520 unsigned n_position_type = this->nnodal_position_type();
1521
1522#ifdef PARANOID
1523 if (n_position_type != 1)
1524 {
1525 throw OomphLibError(
1526 "ImposeDisplacementByLagrangeMultiplierElement cannot (currently) be "
1527 "used with elements that have generalised positional dofs",
1528 OOMPH_CURRENT_FUNCTION,
1529 OOMPH_EXCEPTION_LOCATION);
1530 }
1531#endif
1532
1533 // Find out how many nodes there are
1534 unsigned n_node = nnode();
1535
1536 // Dimension of element
1537 unsigned dim_el = dim();
1538
1539 // Set up memory for the shape functions
1540 Shape psi(n_node);
1541 DShape dpsids(n_node, dim_el);
1542
1543 // Set the value of n_intpt
1544 unsigned n_intpt = integral_pt()->nweight();
1545
1546
1547 // Initialise error
1548 double squared_error = 0.0;
1549
1550 // Loop over the integration points
1551 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1552 {
1553 // Get the integral weight
1554 double w = integral_pt()->weight(ipt);
1555
1556 // Only need to call the local derivatives
1557 dshape_local_at_knot(ipt, psi, dpsids);
1558
1559 // Calculate the Eulerian coordinates and Lagrange multiplier
1560 Vector<double> x(dim_el + 1, 0.0);
1561 Vector<double> lambda(dim_el + 1, 0.0);
1562 Vector<double> zeta(dim_el, 0.0);
1563 DenseMatrix<double> interpolated_a(dim_el, dim_el + 1, 0.0);
1564
1565 // Loop over nodes
1566 for (unsigned j = 0; j < n_node; j++)
1567 {
1568 Node* nod_pt = node_pt(j);
1569
1570 // Cast to a boundary node
1571 BoundaryNodeBase* bnod_pt =
1572 dynamic_cast<BoundaryNodeBase*>(node_pt(j));
1573
1574 // Get the index of the first nodal value associated with
1575 // this FaceElement
1576 unsigned first_index =
1578
1579 // Assemble higher-dimensional quantities
1580 for (unsigned i = 0; i < dim_el + 1; i++)
1581 {
1582 x[i] += nodal_position(j, i) * psi(j);
1583 lambda[i] += nod_pt->value(first_index + i) * psi(j);
1584 for (unsigned ii = 0; ii < dim_el; ii++)
1585 {
1586 interpolated_a(ii, i) +=
1587 lagrangian_position(j, i) * dpsids(j, ii);
1588 }
1589 }
1590 if (!Sparsify)
1591 {
1592 for (unsigned k = 0; k < n_position_type; k++)
1593 {
1594 // Assemble in-element quantities: boundary coordinate
1595 for (unsigned i = 0; i < dim_el; i++)
1596 {
1597 zeta[i] += zeta_nodal(j, k, i) * psi(j, k);
1598 }
1599 }
1600 }
1601 }
1602
1603 if (Sparsify) zeta = Zeta_sub_geom_object[ipt];
1604
1605
1606 // Now find the local undeformed metric tensor from the tangent Vectors
1607 DenseMatrix<double> a(dim_el);
1608 for (unsigned i = 0; i < dim_el; i++)
1609 {
1610 for (unsigned j = 0; j < dim_el; j++)
1611 {
1612 // Initialise surface metric tensor to zero
1613 a(i, j) = 0.0;
1614 // Take the dot product
1615 for (unsigned k = 0; k < dim_el + 1; k++)
1616 {
1617 a(i, j) += interpolated_a(i, k) * interpolated_a(j, k);
1618 }
1619 }
1620 }
1621
1622
1623 // Find the determinant of the metric tensor
1624 double adet = 0.0;
1625 switch (dim_el + 1)
1626 {
1627 case 2:
1628 adet = a(0, 0);
1629 break;
1630
1631 case 3:
1632 adet = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
1633 break;
1634
1635 default:
1636 throw OomphLibError(
1637 "Wrong dimension "
1638 "fill_in_generic_contribution_to_residuals_displ_lagr_multiplier",
1639 "ImposeDisplacementByLagrangeMultiplierElement::fill_in_generic_"
1640 "contribution_to_residuals_displ_lagr_multiplier()",
1641 OOMPH_EXCEPTION_LOCATION);
1642 }
1643
1644 // Get prescribed wall shape
1645 Vector<double> r_prescribed(dim_el + 1);
1646 if (!Sparsify)
1647 {
1648 Boundary_shape_geom_object_pt->position(zeta, r_prescribed);
1649 }
1650 else
1651 {
1652 Sub_geom_object_pt[ipt]->position(zeta, r_prescribed);
1653 }
1654
1655 // Premultiply the weights and the square-root of the determinant of
1656 // the metric tensor
1657 double W = w * sqrt(adet);
1658
1659 // Assemble error
1660
1661 // Loop over directions
1662 for (unsigned i = 0; i < dim_el + 1; i++)
1663 {
1664 squared_error +=
1665 (x[i] - r_prescribed[i]) * (x[i] - r_prescribed[i]) * W;
1666 }
1667 } // End of loop over the integration points
1668
1669 return squared_error;
1670 }
1671
1672
1673 protected:
1674 /// Helper function to compute the residuals and, if flag==1, the
1675 /// Jacobian
1677 Vector<double>& residuals,
1678 DenseMatrix<double>& jacobian,
1679 const unsigned& flag)
1680 {
1681 // Find out how many positional dofs there are
1682 unsigned n_position_type = this->nnodal_position_type();
1683
1684#ifdef PARANOID
1685 if (n_position_type != 1)
1686 {
1687 throw OomphLibError(
1688 "ImposeDisplacementByLagrangeMultiplierElement cannot (currently) be "
1689 "used with elements that have generalised positional dofs",
1690 OOMPH_CURRENT_FUNCTION,
1691 OOMPH_EXCEPTION_LOCATION);
1692 }
1693#endif
1694
1695 // Find out how many nodes there are
1696 unsigned n_node = nnode();
1697
1698 // Dimension of element
1699 unsigned dim_el = dim();
1700
1701 // Set up memory for the shape functions
1702 Shape psi(n_node);
1703 DShape dpsids(n_node, dim_el);
1704
1705 // Set the value of n_intpt
1706 unsigned n_intpt = integral_pt()->nweight();
1707
1708 // Loop over the integration points
1709 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
1710 {
1711 // Get the integral weight
1712 double w = integral_pt()->weight(ipt);
1713
1714 // Only need to call the local derivatives
1715 dshape_local_at_knot(ipt, psi, dpsids);
1716
1717 // Calculate the Eulerian coordinates and Lagrange multiplier
1718 Vector<double> x(dim_el + 1, 0.0);
1719 Vector<double> lambda(dim_el + 1, 0.0);
1720 Vector<double> zeta(dim_el, 0.0);
1721 DenseMatrix<double> interpolated_a(dim_el, dim_el + 1, 0.0);
1722
1723 // Loop over nodes
1724 for (unsigned j = 0; j < n_node; j++)
1725 {
1726 Node* nod_pt = node_pt(j);
1727
1728 // Cast to a boundary node
1729 BoundaryNodeBase* bnod_pt =
1730 dynamic_cast<BoundaryNodeBase*>(node_pt(j));
1731
1732 // Get the index of the first nodal value associated with
1733 // this FaceElement
1734 unsigned first_index =
1736
1737 // Assemble higher-dimensional quantities
1738 for (unsigned i = 0; i < dim_el + 1; i++)
1739 {
1740 x[i] += nodal_position(j, i) * psi(j);
1741 lambda[i] += nod_pt->value(first_index + i) * psi(j);
1742 for (unsigned ii = 0; ii < dim_el; ii++)
1743 {
1744 interpolated_a(ii, i) +=
1745 lagrangian_position(j, i) * dpsids(j, ii);
1746 }
1747 }
1748 if (!Sparsify)
1749 {
1750 for (unsigned k = 0; k < n_position_type; k++)
1751 {
1752 // Assemble in-element quantities: boundary coordinate
1753 for (unsigned i = 0; i < dim_el; i++)
1754 {
1755 zeta[i] += zeta_nodal(j, k, i) * psi(j, k);
1756 }
1757 }
1758 }
1759 }
1760
1761 if (Sparsify) zeta = Zeta_sub_geom_object[ipt];
1762
1763
1764 // Now find the local undeformed metric tensor from the tangent Vectors
1765 DenseMatrix<double> a(dim_el);
1766 for (unsigned i = 0; i < dim_el; i++)
1767 {
1768 for (unsigned j = 0; j < dim_el; j++)
1769 {
1770 // Initialise surface metric tensor to zero
1771 a(i, j) = 0.0;
1772 // Take the dot product
1773 for (unsigned k = 0; k < dim_el + 1; k++)
1774 {
1775 a(i, j) += interpolated_a(i, k) * interpolated_a(j, k);
1776 }
1777 }
1778 }
1779
1780
1781 // Find the determinant of the metric tensor
1782 double adet = 0.0;
1783 switch (dim_el + 1)
1784 {
1785 case 2:
1786 adet = a(0, 0);
1787 break;
1788
1789 case 3:
1790 adet = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
1791 break;
1792
1793 default:
1794 throw OomphLibError(
1795 "Wrong dimension "
1796 "fill_in_generic_contribution_to_residuals_displ_lagr_multiplier",
1797 "ImposeDisplacementByLagrangeMultiplierElement::fill_in_generic_"
1798 "contribution_to_residuals_displ_lagr_multiplier()",
1799 OOMPH_EXCEPTION_LOCATION);
1800 }
1801
1802 // Get prescribed wall shape
1803 Vector<double> r_prescribed(dim_el + 1);
1804 if (!Sparsify)
1805 {
1806 Boundary_shape_geom_object_pt->position(zeta, r_prescribed);
1807 }
1808 else
1809 {
1810 Sub_geom_object_pt[ipt]->position(zeta, r_prescribed);
1811 }
1812
1813 // Premultiply the weights and the square-root of the determinant of
1814 // the metric tensor
1815 double W = w * sqrt(adet);
1816
1817 // Assemble residuals and jacobian
1818
1819 // Loop over directions
1820 for (unsigned i = 0; i < dim_el + 1; i++)
1821 {
1822 // Loop over the nodes
1823 for (unsigned j = 0; j < n_node; j++)
1824 {
1825 // Assemble residual for Lagrange multiplier:
1826
1827 // Cast to a boundary node
1828 BoundaryNodeBase* bnod_pt =
1829 dynamic_cast<BoundaryNodeBase*>(node_pt(j));
1830
1831 // Local eqn number:
1832 int local_eqn = nodal_local_eqn(
1833 j,
1835
1836
1837 if (local_eqn >= 0)
1838 {
1839 residuals[local_eqn] += (x[i] - r_prescribed[i]) * psi(j) * W;
1840
1841 // Do Jacobian too?
1842 if (flag == 1)
1843 {
1844 // Loop over the nodes again for unknowns (only diagonal
1845 // contribution to direction!).
1846 for (unsigned jj = 0; jj < n_node; jj++)
1847 {
1848 int local_unknown = position_local_eqn(jj, 0, i);
1849 if (local_unknown >= 0)
1850 {
1851 jacobian(local_eqn, local_unknown) += psi(jj) * psi(j) * W;
1852 }
1853 }
1854 }
1855 }
1856
1857
1858 // Add Lagrange multiplier contribution to bulk equations
1859
1860 // Local eqn number: Node, type, direction
1861 local_eqn = position_local_eqn(j, 0, i);
1862 if (local_eqn >= 0)
1863 {
1864 // Add to residual
1865 residuals[local_eqn] += lambda[i] * psi(j) * W;
1866
1867 // Do Jacobian too?
1868 if (flag == 1)
1869 {
1870 // Loop over the nodes again for unknowns (only diagonal
1871 // contribution to direction!).
1872 for (unsigned jj = 0; jj < n_node; jj++)
1873 {
1874 // Cast to a boundary node
1875 BoundaryNodeBase* bnode_pt =
1876 dynamic_cast<BoundaryNodeBase*>(node_pt(jj));
1877
1878 int local_unknown = nodal_local_eqn(
1879 jj,
1881 Id) +
1882 i);
1883 if (local_unknown >= 0)
1884 {
1885 jacobian(local_eqn, local_unknown) += psi(jj) * psi(j) * W;
1886 }
1887 }
1888 }
1889 }
1890 }
1891 }
1892
1893
1894 } // End of loop over the integration points
1895 }
1896
1897
1898 /// The number of "DOF types" that degrees of freedom in this element
1899 /// are sub-divided into: We only label the
1900 /// Lagrange multiplier degrees of freedom (one for each spatial dimension)
1901 unsigned ndof_types() const
1902 {
1903 return this->dim() + 1;
1904 }
1905
1906
1907 /// Create a list of pairs for all unknowns in this element,
1908 /// so that the first entry in each pair contains the global equation
1909 /// number of the unknown, while the second one contains the number
1910 /// of the dof that this unknown is associated with.
1911 /// (Function can obviously only be called if the equation numbering
1912 /// scheme has been set up.)
1914 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
1915 {
1916 // temporary pair (used to store dof lookup prior to being added to list)
1917 std::pair<unsigned, unsigned> dof_lookup;
1918
1919 // number of nodes
1920 const unsigned n_node = this->nnode();
1921
1922 // Loop over directions
1923 unsigned dim_el = this->dim();
1924 for (unsigned i = 0; i < dim_el + 1; i++)
1925 {
1926 // Loop over the nodes
1927 for (unsigned j = 0; j < n_node; j++)
1928 {
1929 // Cast to a boundary node
1930 BoundaryNodeBase* bnod_pt =
1931 dynamic_cast<BoundaryNodeBase*>(node_pt(j));
1932
1933 // Local eqn number:
1934 int local_eqn = nodal_local_eqn(
1936 if (local_eqn >= 0)
1937 {
1938 // store dof lookup in temporary pair: First entry in pair
1939 // is global equation number; second entry is dof type
1940 dof_lookup.first = this->eqn_number(local_eqn);
1941 dof_lookup.second = i;
1942
1943 // add to list
1944 dof_lookup_list.push_front(dof_lookup);
1945 }
1946 }
1947 }
1948 }
1949
1950
1951 /// Lagrange Id
1952 unsigned Id;
1953
1954
1955#ifdef PARANOID
1956
1957 /// Bool to record the number of geom Data that has been
1958 /// assigned to external data (we're keeping a record to make
1959 /// sure we're not accidentally wiping more than we assigned). Only
1960 /// included if compiled with PARANOID switched on.
1962
1963#endif
1964
1965 /// GeomObject that specifies the prescribed
1966 /// boundary displacement; GeomObject is assumed to be
1967 /// parametrised by the same coordinate the is used as
1968 /// the boundary coordinate in the bulk solid mesh to which
1969 /// this element is attached.
1971
1972 /// Storage for sub-GeomObject at the integration points
1974
1975 /// Storage for local coordinates in sub-GeomObjects at integration
1976 /// points
1978
1979 /// Boolean flag to indicate if we're using geometric Data of
1980 /// sub-GeomObjects that make up the (possibly compound) GeomObject
1981 /// that specifies the boundary shape. Defaults to true.
1983 };
1984
1985
1986 /// /////////////////////////////////////////////////////////////////////////
1987 /// /////////////////////////////////////////////////////////////////////////
1988 /// /////////////////////////////////////////////////////////////////////////
1989
1990
1991 //======================================================================
1992 /// A class for elements that allow the imposition of a displacement
1993 /// constraint for "bulk" solid elements via a Lagrange multiplier.
1994 /// The geometrical information can be read from the FaceGeometry<ELEMENT>
1995 /// class and and thus, we can be generic enough without the need to have
1996 /// a separate equations class.
1997 /// \b NOTE: Currently (and for the foreseeable future) this
1998 /// element only works with bulk elements that do not have
1999 /// generalised degrees of freedom (so it won't work with
2000 /// Hermite-type elements, say). The additional functionality
2001 /// to deal with such elements could easily be added (once a
2002 /// a suitable test case is written). For now we simply throw
2003 /// errors if an attempt is made to use the element with an unsuitable
2004 /// bulk element.
2005 ///
2006 /// REFINEABLE VERSION
2007 //======================================================================
2008 template<class ELEMENT>
2010 : public virtual ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>,
2012
2013 {
2014 public:
2015 /// Constructor takes a "bulk" element and the
2016 /// index that identifies which face the FaceElement is supposed
2017 /// to be attached to. The optional identifier can be used
2018 /// to distinguish the additional nodal values created by
2019 /// this element from thos created by other FaceElements.
2021 FiniteElement* const& element_pt,
2022 const int& face_index,
2023 const unsigned& id = 0)
2025 element_pt, face_index, id, true)
2026 {
2027 }
2028
2029
2030 /// Number of continuously interpolated values: Same for
2031 /// all nodes since it's a refineable element
2033 {
2034 return node_pt(0)->nvalue();
2035 }
2036
2037 /// Fill in the residuals
2039 {
2040 // Call the generic routine with the flag set to 0
2042 residuals, GeneralisedElement::Dummy_matrix, 0);
2043 }
2044
2045
2046 /// Fill in contribution from Jacobian
2048 DenseMatrix<double>& jacobian)
2049 {
2050 // Call the generic routine with the flag set to 1
2052 residuals, jacobian, 1);
2053
2054 // Fill in the derivatives w.r.t. external data by FD, re-using
2055 // the pre-computed residual vector
2056 this->fill_in_jacobian_from_external_by_fd(residuals, jacobian);
2057 }
2058
2059
2060 protected:
2061 /// Helper function to compute the residuals and, if flag==1, the
2062 /// Jacobian
2064 Vector<double>& residuals,
2065 DenseMatrix<double>& jacobian,
2066 const unsigned& flag)
2067 {
2068 // Find out how many positional dofs there are
2069 unsigned n_position_type = this->nnodal_position_type();
2070
2071#ifdef PARANOID
2072 if (n_position_type != 1)
2073 {
2074 throw OomphLibError(
2075 "RefineableImposeDisplacementByLagrangeMultiplierElement \n cannot "
2076 "(currently) be used with elements that have generalised\n "
2077 "positional dofs. Upgrade should be straightforward though the code "
2078 "is in a \n bit of state, with generalised degrees of freedom "
2079 "sometimes half taken into \n account, sometimes completely "
2080 "ignored...",
2081 OOMPH_CURRENT_FUNCTION,
2082 OOMPH_EXCEPTION_LOCATION);
2083 }
2084#endif
2085
2086 // Find out how many nodes there are
2087 unsigned n_node = this->nnode();
2088
2089 // Dimension of element
2090 unsigned dim_el = this->dim();
2091
2092 // Set up memory for the shape functions
2093 Shape psi(n_node);
2094 DShape dpsids(n_node, dim_el);
2095
2096 // Set the value of n_intpt
2097 unsigned n_intpt = this->integral_pt()->nweight();
2098
2099
2100 // Integers to store local equation number
2101 int local_eqn = 0;
2102 int local_unknown = 0;
2103
2104 // Loop over the integration points
2105 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2106 {
2107 // Get the integral weight
2108 double w = this->integral_pt()->weight(ipt);
2109
2110 // Only need to call the local derivatives
2111 this->dshape_local_at_knot(ipt, psi, dpsids);
2112
2113 // Calculate the Eulerian coordinates and Lagrange multiplier
2114 Vector<double> x(dim_el + 1, 0.0);
2115 Vector<double> lambda(dim_el + 1, 0.0);
2116 Vector<double> zeta(dim_el, 0.0);
2117 DenseMatrix<double> interpolated_a(dim_el, dim_el + 1, 0.0);
2118
2119 // Loop over nodes -- note in these calculations hang-ness is
2120 // automatically taken into account because of calls to position(...)
2121 // etc
2122 for (unsigned j = 0; j < n_node; j++)
2123 {
2124 Node* nod_pt = this->node_pt(j);
2125
2126 // Cast to a boundary node
2127 BoundaryNodeBase* bnod_pt = dynamic_cast<BoundaryNodeBase*>(nod_pt);
2128
2129 // Get the index of the first nodal value associated with
2130 // this FaceElement
2131 unsigned first_index =
2133
2134 // Assemble higher-dimensional quantities
2135 for (unsigned i = 0; i < dim_el + 1; i++)
2136 {
2137 x[i] += this->nodal_position(j, i) * psi(j);
2138 lambda[i] += nod_pt->value(first_index + i) * psi(j);
2139 for (unsigned ii = 0; ii < dim_el; ii++)
2140 {
2141 interpolated_a(ii, i) +=
2142 this->lagrangian_position(j, i) * dpsids(j, ii);
2143 }
2144 }
2145 if (!this->Sparsify)
2146 {
2147 for (unsigned k = 0; k < n_position_type; k++)
2148 {
2149 // Assemble in-element quantities: boundary coordinate
2150 for (unsigned i = 0; i < dim_el; i++)
2151 {
2152 zeta[i] += this->zeta_nodal(j, k, i) * psi(j, k);
2153 }
2154 }
2155 }
2156 }
2157
2158 if (this->Sparsify) zeta = this->Zeta_sub_geom_object[ipt];
2159
2160
2161 // Now find the local undeformed metric tensor from the tangent Vectors
2162 DenseMatrix<double> a(dim_el);
2163 for (unsigned i = 0; i < dim_el; i++)
2164 {
2165 for (unsigned j = 0; j < dim_el; j++)
2166 {
2167 // Initialise surface metric tensor to zero
2168 a(i, j) = 0.0;
2169 // Take the dot product
2170 for (unsigned k = 0; k < dim_el + 1; k++)
2171 {
2172 a(i, j) += interpolated_a(i, k) * interpolated_a(j, k);
2173 }
2174 }
2175 }
2176
2177
2178 // Find the determinant of the metric tensor
2179 double adet = 0.0;
2180 switch (dim_el + 1)
2181 {
2182 case 2:
2183 adet = a(0, 0);
2184 break;
2185
2186 case 3:
2187 adet = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
2188 break;
2189
2190 default:
2191 throw OomphLibError(
2192 "Wrong dimension "
2193 "refineable_fill_in_generic_contribution_to_residuals_displ_lagr_"
2194 "multiplier",
2195 "RefineableImposeDisplacementByLagrangeMultiplierElement::"
2196 "refineablefill_in_generic_contribution_to_residuals_displ_lagr_"
2197 "multiplier()",
2198 OOMPH_EXCEPTION_LOCATION);
2199 }
2200
2201 // Get prescribed wall shape
2202 Vector<double> r_prescribed(dim_el + 1);
2203 if (!this->Sparsify)
2204 {
2205 this->Boundary_shape_geom_object_pt->position(zeta, r_prescribed);
2206 }
2207 else
2208 {
2209 this->Sub_geom_object_pt[ipt]->position(zeta, r_prescribed);
2210 }
2211
2212 // Premultiply the weights and the square-root of the determinant of
2213 // the metric tensor
2214 double W = w * sqrt(adet);
2215
2216 // Assemble residuals and jacobian
2217
2218
2219 // Number of master nodes and storage for the weight of the shape
2220 // function
2221 unsigned n_master = 1;
2222 unsigned n_master2 = 1;
2223 double hang_weight = 1.0;
2224 double hang_weight2 = 1.0;
2225
2226 // Pointer to hang info object
2227 HangInfo* hang_info_pt = 0;
2228 HangInfo* hang_info2_pt = 0;
2229
2230
2231 // Loop over nodes
2232 for (unsigned j = 0; j < n_node; j++)
2233 {
2234 // Local node itself (hanging or not)
2235 Node* local_node_pt = this->node_pt(j);
2236
2237 // Local boolean to indicate whether the node is hanging
2238 bool is_node_hanging = local_node_pt->is_hanging();
2239
2240 // If the node is hanging
2241 if (is_node_hanging)
2242 {
2243 hang_info_pt = local_node_pt->hanging_pt();
2244
2245 // Read out number of master nodes from hanging data
2246 n_master = hang_info_pt->nmaster();
2247 }
2248 // Otherwise the node is its own master
2249 else
2250 {
2251 n_master = 1;
2252 }
2253
2254 // Loop over the master nodes
2255 for (unsigned m = 0; m < n_master; m++)
2256 {
2257 // Loop over velocity components for equations
2258 for (unsigned i = 0; i < dim_el + 1; i++)
2259 {
2260 // Get the equation number for Lagrange multiplier eqn
2261
2262 // If the node is hanging
2263 if (is_node_hanging)
2264 {
2265 // Cast to a boundary node
2266 BoundaryNodeBase* bnod_pt = dynamic_cast<BoundaryNodeBase*>(
2267 hang_info_pt->master_node_pt(m));
2268
2269 // Get the equation number from the master node
2270 local_eqn = this->local_hang_eqn(
2271 hang_info_pt->master_node_pt(m),
2273 this->Id) +
2274 i);
2275
2276 // Get the hang weight from the master node
2277 hang_weight = hang_info_pt->master_weight(m);
2278 }
2279 // If the node is not hanging
2280 else
2281 {
2282 // Cast to a boundary node
2283 BoundaryNodeBase* bnod_pt =
2284 dynamic_cast<BoundaryNodeBase*>(local_node_pt);
2285
2286 // Local equation number
2287 local_eqn = this->nodal_local_eqn(
2288 j,
2290 this->Id) +
2291 i);
2292
2293 // Node contributes with full weight
2294 hang_weight = 1.0;
2295 }
2296
2297 // If it's not a boundary condition...
2298 if (local_eqn >= 0)
2299 {
2300 residuals[local_eqn] +=
2301 (x[i] - r_prescribed[i]) * psi(j) * W * hang_weight;
2302
2303 // Do Jacobian too?
2304 if (flag == 1)
2305 {
2306 // Loop over the nodes again for unknowns (only diagonal
2307 // contribution to direction!).
2308 for (unsigned jj = 0; jj < n_node; jj++)
2309 {
2310 // Local node itself (hanging or not)
2311 Node* local_node2_pt = this->node_pt(jj);
2312
2313 // Local boolean to indicate whether the node is hanging
2314 bool is_node_hanging2 = local_node2_pt->is_hanging();
2315
2316 // If the node is hanging
2317 if (is_node_hanging2)
2318 {
2319 hang_info2_pt = local_node2_pt->hanging_pt();
2320
2321 // Read out number of master nodes from hanging data
2322 n_master2 = hang_info2_pt->nmaster();
2323 }
2324 // Otherwise the node is its own master
2325 else
2326 {
2327 n_master2 = 1;
2328 }
2329
2330 // Loop over the master nodes
2331 for (unsigned m2 = 0; m2 < n_master2; m2++)
2332 {
2333 // Storage for local equation numbers at node indexed by
2334 // type and direction
2335 DenseMatrix<int> position_local_eqn_at_node2(
2336 n_position_type, dim_el + 1);
2337
2338 if (is_node_hanging2)
2339 {
2340 // Find the equation numbers
2341 position_local_eqn_at_node2 = local_position_hang_eqn(
2342 local_node2_pt->hanging_pt()->master_node_pt(m2));
2343
2344 // Find the hanging node weight
2345 hang_weight2 =
2346 local_node2_pt->hanging_pt()->master_weight(m2);
2347 }
2348 else
2349 {
2350 // Non-loop of types of dofs
2351 // for(unsigned k2=0;k2<n_position_type;k2++)
2352 // {
2353 unsigned k2 = 0;
2354
2355 // Loop over the displacement components
2356 // for(unsigned i2=0;i2<dim_el+1;i2++)
2357 unsigned i2 =
2358 i; // only need that one, but need to store
2359 // information in this container because
2360 // it's required for hanging case.
2361 {
2362 position_local_eqn_at_node2(k2, i2) =
2363 this->position_local_eqn(jj, k2, i2);
2364 }
2365
2366 // Hang weight is one
2367 hang_weight2 = 1.0;
2368 }
2369
2370 unsigned k2 = 0;
2371 local_unknown = position_local_eqn_at_node2(k2, i);
2372 if (local_unknown >= 0)
2373 {
2374 jacobian(local_eqn, local_unknown) +=
2375 psi(jj) * psi(j) * W * hang_weight * hang_weight2;
2376 }
2377 }
2378 }
2379 }
2380
2381
2382 // Add Lagrange multiplier contribution to bulk equations
2383
2384
2385 // Storage for local equation numbers at node indexed by
2386 // type and direction
2387 DenseMatrix<int> position_local_eqn_at_node(n_position_type,
2388 dim_el + 1);
2389
2390 if (is_node_hanging)
2391 {
2392 // Find the equation numbers
2393 position_local_eqn_at_node = local_position_hang_eqn(
2394 local_node_pt->hanging_pt()->master_node_pt(m));
2395 }
2396 else
2397 {
2398 // Non-loop of types of dofs
2399 // for(unsigned k2=0;k2<n_position_type;k2++)
2400 // {
2401 unsigned k2 = 0;
2402
2403 // Loop over the displacement components
2404 // for(unsigned i2=0;i2<dim_el+1;i2++)
2405 unsigned i2 = i; // only need that one, but need to store
2406 // information in this container because
2407 // it's required for hanging case.
2408 {
2409 position_local_eqn_at_node(k2, i2) =
2410 this->position_local_eqn(j, k2, i2);
2411 }
2412 }
2413 unsigned k = 0;
2414 local_eqn = position_local_eqn_at_node(k, i);
2415
2416 /*IF it's not a boundary condition*/
2417 if (local_eqn >= 0)
2418 {
2419 // Add to residual
2420 residuals[local_eqn] += lambda[i] * psi(j) * W * hang_weight;
2421
2422 // Do Jacobian too?
2423 if (flag == 1)
2424 {
2425 // Loop over the nodes again for unknowns (only diagonal
2426 // contribution to direction!).
2427 for (unsigned jj = 0; jj < n_node; jj++)
2428 {
2429 // Local node itself (hanging or not)
2430 Node* local_node2_pt = this->node_pt(jj);
2431
2432 // Local boolean to indicate whether the node is hanging
2433 bool is_node_hanging2 = local_node2_pt->is_hanging();
2434
2435 // If the node is hanging
2436 if (is_node_hanging2)
2437 {
2438 hang_info2_pt = local_node2_pt->hanging_pt();
2439
2440 // Read out number of master nodes from hanging data
2441 n_master2 = hang_info2_pt->nmaster();
2442 }
2443 // Otherwise the node is its own master
2444 else
2445 {
2446 n_master2 = 1;
2447 }
2448
2449 // Loop over the master nodes
2450 for (unsigned m2 = 0; m2 < n_master2; m2++)
2451 {
2452 // Get the equation number for Lagrange multiplier eqn
2453
2454 // If the node is hanging
2455 if (is_node_hanging2)
2456 {
2457 // Cast to a boundary node
2458 BoundaryNodeBase* bnod2_pt =
2459 dynamic_cast<BoundaryNodeBase*>(
2460 hang_info2_pt->master_node_pt(m2));
2461
2462 // Get the equation number from the master node
2463 local_unknown = this->local_hang_eqn(
2464 hang_info2_pt->master_node_pt(m2),
2465 bnod2_pt
2467 this->Id) +
2468 i);
2469
2470 // Get the hang weight from the master node
2471 hang_weight2 = hang_info2_pt->master_weight(m2);
2472 }
2473 // If the node is not hanging
2474 else
2475 {
2476 // Cast to a boundary node
2477 BoundaryNodeBase* bnod2_pt =
2478 dynamic_cast<BoundaryNodeBase*>(local_node2_pt);
2479
2480 // Local equation number
2481 local_unknown = this->nodal_local_eqn(
2482 jj,
2483 bnod2_pt
2484 ->index_of_first_value_assigned_by_face_element(
2485 this->Id) +
2486 i);
2487
2488 // Node contributes with full weight
2489 hang_weight2 = 1.0;
2490 }
2491
2492 if (local_unknown >= 0)
2493 {
2494 jacobian(local_eqn, local_unknown) +=
2495 psi(jj) * psi(j) * W * hang_weight * hang_weight2;
2496 }
2497 }
2498 }
2499 }
2500 }
2501 }
2502 }
2503 }
2504 }
2505
2506 } // End of loop over the integration points
2507 }
2508 };
2509
2510
2511 /// /////////////////////////////////////////////////////////////////////////
2512 /// /////////////////////////////////////////////////////////////////////////
2513 /// /////////////////////////////////////////////////////////////////////////
2514
2515
2516 //======================================================================
2517 /// A class for elements that allow the imposition of a displacement
2518 /// constraint for bulk solid elements via a Lagrange multiplier.
2519 /// Prescribed displaced is obtained from an adjacent bulk solid
2520 /// element (rather than from a lower-dimensional GeomObject
2521 /// as in the corresponding ImposeDisplacementByLagrangeMultiplierElement
2522 /// class. The present class is particularly suited for parallel
2523 /// FSI computations.
2524 /// \b NOTE: Currently (and for the foreseeable future) this
2525 /// element only works with bulk elements that do not have
2526 /// generalised degrees of freedom (so it won't work with
2527 /// Hermite-type elements, say). The additional functionality
2528 /// to deal with such elements could easily be added (once a
2529 /// a suitable test case is written). For now we simply throw
2530 /// errors if an attempt is made to use the element with an unsuitable
2531 /// bulk element.
2532 //======================================================================
2533 template<class ELEMENT>
2535 : public virtual FaceGeometry<ELEMENT>,
2536 public virtual SolidFaceElement,
2537 public virtual ElementWithExternalElement
2538 {
2539 public:
2540 /// Function to describe the local dofs of the elements. The ostream
2541 /// specifies the output stream to which the description
2542 /// is written; the string stores the currently
2543 /// assembled output that is ultimately written to the
2544 /// output stream by Data::describe_dofs(...); it is typically
2545 /// built up incrementally as we descend through the
2546 /// call hierarchy of this function when called from
2547 /// Problem::describe_dofs(...)
2548 void describe_local_dofs(std::ostream& out,
2549 const std::string& current_string) const
2550 {
2552 describe_nodal_local_dofs(out, current_string);
2553 }
2554
2555 /// Constructor takes a "bulk" element and the
2556 /// index that identifies which face the FaceElement is supposed
2557 /// to be attached to. The optional identifier can be used
2558 /// to distinguish the additional nodal values created by
2559 /// this element from thos created by other FaceElements.
2561 FiniteElement* const& element_pt,
2562 const int& face_index,
2563 const unsigned& id = 0,
2564 const bool& called_from_refineable_constructor = false)
2566 {
2567 // Set external element storage - one interaction
2568 this->set_ninteraction(1);
2569
2570 // Store the ID of the FaceElement -- this is used to distinguish
2571 // it from any others
2572 Id = id;
2573
2574 // Build the face element
2575 element_pt->build_face_element(face_index, this);
2576
2577#ifdef PARANOID
2578 {
2579 // Check that the bulk element is not a refineable 3d element
2580 if (!called_from_refineable_constructor)
2581 {
2582 if (element_pt->dim() == 3)
2583 {
2584 // Is it refineable
2585 RefineableElement* ref_el_pt =
2586 dynamic_cast<RefineableElement*>(element_pt);
2587 if (ref_el_pt != 0)
2588 {
2589 if (this->has_hanging_nodes())
2590 {
2591 throw OomphLibError(
2592 "This face element will not work correctly if nodes are "
2593 "hanging\nUse the refineable version instead. ",
2594 OOMPH_CURRENT_FUNCTION,
2595 OOMPH_EXCEPTION_LOCATION);
2596 }
2597 }
2598 }
2599 }
2600 }
2601
2602 {
2603 // Check that the bulk element does not require generalised positional
2604 // degrees of freedom
2605 if (element_pt->nnodal_position_type() != 1)
2606 {
2607 throw OomphLibError("FSIImposeDisplacementByLagrangeMultiplierElement"
2608 " cannot (currently) be used with elements that "
2609 "have generalised positional dofs",
2610 OOMPH_CURRENT_FUNCTION,
2611 OOMPH_EXCEPTION_LOCATION);
2612 }
2613 }
2614#endif
2615
2616 // Dimension of the bulk element
2617 unsigned dim = element_pt->dim();
2618
2619 // We need dim additional values for each FaceElement node
2620 // to store the dim Lagrange multipliers.
2621 Vector<unsigned> n_additional_values(nnode(), dim);
2622
2623 // Now add storage for Lagrange multipliers and set the map containing
2624 // the position of the first entry of this face element's
2625 // additional values.
2626 add_additional_values(n_additional_values, id);
2627 }
2628
2629
2630 /// Fill in the residuals
2632 {
2633 // Call the generic routine with the flag set to 0
2635 residuals, GeneralisedElement::Dummy_matrix, 0);
2636 }
2637
2638
2639 /// Fill in contribution from Jacobian
2641 DenseMatrix<double>& jacobian)
2642 {
2643 // Call the generic routine with the flag set to 1
2645 residuals, jacobian, 1);
2646
2647 // Add the contribution of the external interaction by finite differences
2649 }
2650
2651 /// Fill in contribution to Mass matrix and
2652 /// Jacobian. There is no contributiont to mass matrix
2653 /// so simply call the fill_in_contribution_to_jacobian term
2654 /// Note that the Jacobian is multiplied by minus one to
2655 /// ensure that the mass matrix is positive semi-definite.
2657 Vector<double>& residuals,
2658 DenseMatrix<double>& jacobian,
2659 DenseMatrix<double>& mass_matrix)
2660 {
2661 // Just call the jacobian calculation
2662 fill_in_contribution_to_jacobian(residuals, jacobian);
2663
2664 // Multiply the residuals and jacobian by minus one
2665 const unsigned n_dof = this->ndof();
2666 for (unsigned i = 0; i < n_dof; i++)
2667 {
2668 residuals[i] *= -1.0;
2669 for (unsigned j = 0; j < n_dof; j++)
2670 {
2671 jacobian(i, j) *= -1.0;
2672 }
2673 }
2674 }
2675
2676
2677 /// Output function
2678 void output(std::ostream& outfile, const unsigned& n_plot)
2679 {
2680 // Elemental dimension
2681 unsigned dim_el = dim();
2682
2683 // Find the number of positional types
2684 unsigned n_position_type = this->nnodal_position_type();
2685
2686#ifdef PARANOID
2687 if (n_position_type != 1)
2688 {
2689 throw OomphLibError(
2690 "FSIImposeDisplacementByLagrangeMultiplierElement cannot (currently) "
2691 "be used with elements that have generalised positional dofs",
2692 OOMPH_CURRENT_FUNCTION,
2693 OOMPH_EXCEPTION_LOCATION);
2694 }
2695#endif
2696
2697
2698 // Local coord
2699 Vector<double> s(dim_el);
2700
2701 // # of nodes,
2702 unsigned n_node = nnode();
2703 Shape psi(n_node, n_position_type);
2704
2705 // Tecplot header info
2706 outfile << this->tecplot_zone_string(n_plot);
2707
2708 // Loop over plot points
2709 unsigned num_plot_points = this->nplot_points(n_plot);
2710 for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
2711 {
2712 // Get local coordinates of plot point
2713 this->get_s_plot(iplot, n_plot, s);
2714
2715 // Get shape function
2716 shape(s, psi);
2717
2718 // Calculate the Eulerian coordinates and Lagrange multiplier
2719 Vector<double> x(dim_el + 1, 0.0);
2720 Vector<double> lambda(dim_el + 1, 0.0);
2721 Vector<double> zeta(dim_el, 0.0);
2722 for (unsigned j = 0; j < n_node; j++)
2723 {
2724 // Cast to a boundary node
2725 BoundaryNodeBase* bnod_pt =
2726 dynamic_cast<BoundaryNodeBase*>(node_pt(j));
2727
2728 // get the node pt
2729 Node* nod_pt = node_pt(j);
2730
2731 // Get the index of the first nodal value associated with
2732 // this FaceElement
2733 unsigned first_index =
2735
2736 // higher dimensional quantities
2737 for (unsigned i = 0; i < dim_el + 1; i++)
2738 {
2739 x[i] += nodal_position(j, i) * psi(j, 0); // need to sort
2740 // this out properly
2741 // for generalised dofs
2742 lambda[i] += nod_pt->value(first_index + i) * psi(j, 0);
2743 }
2744 // In-element quantities
2745 for (unsigned i = 0; i < dim_el; i++)
2746 {
2747 // Loop over positional types
2748 for (unsigned k = 0; k < n_position_type; k++)
2749 {
2750 zeta[i] += zeta_nodal(j, k, i) * psi(j, k);
2751 }
2752 }
2753 }
2754
2755 // Output stuff
2756 for (unsigned i = 0; i < dim_el + 1; i++)
2757 {
2758 outfile << x[i] << " ";
2759 }
2760 for (unsigned i = 0; i < dim_el + 1; i++)
2761 {
2762 outfile << -lambda[i] << " ";
2763 }
2764 outfile << std::endl;
2765 }
2766 }
2767
2768
2769 /// Output function
2770 void output(std::ostream& outfile)
2771 {
2772 unsigned n_plot = 5;
2773 output(outfile, n_plot);
2774 }
2775
2776
2777 protected:
2778 /// Helper function to compute the residuals and, if flag==1, the
2779 /// Jacobian
2781 Vector<double>& residuals,
2782 DenseMatrix<double>& jacobian,
2783 const unsigned& flag)
2784 {
2785#ifdef PARANOID
2786 // Find out how many positional dofs there are
2787 unsigned n_position_type = this->nnodal_position_type();
2788
2789 if (n_position_type != 1)
2790 {
2791 throw OomphLibError(
2792 "FSIImposeDisplacementByLagrangeMultiplierElement cannot (currently) "
2793 "be used with elements that have generalised positional dofs",
2794 OOMPH_CURRENT_FUNCTION,
2795 OOMPH_EXCEPTION_LOCATION);
2796 }
2797#endif
2798
2799 // Find out how many nodes there are
2800 unsigned n_node = nnode();
2801
2802 // Dimension of element
2803 unsigned dim_el = dim();
2804
2805 // Set up memory for the shape functions
2806 Shape psi(n_node);
2807 DShape dpsids(n_node, dim_el);
2808
2809 // Set the value of n_intpt
2810 unsigned n_intpt = integral_pt()->nweight();
2811
2812 // Loop over the integration points
2813 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
2814 {
2815 // Get the integral weight
2816 double w = integral_pt()->weight(ipt);
2817
2818 // Only need to call the local derivatives
2819 dshape_local_at_knot(ipt, psi, dpsids);
2820
2821 // Calculate the Eulerian coordinates and Lagrange multiplier
2822 Vector<double> x(dim_el + 1, 0.0);
2823 Vector<double> lambda(dim_el + 1, 0.0);
2824 Vector<double> zeta(dim_el, 0.0);
2825 DenseMatrix<double> interpolated_a(dim_el, dim_el + 1, 0.0);
2826
2827 // Loop over nodes
2828 for (unsigned j = 0; j < n_node; j++)
2829 {
2830 Node* nod_pt = node_pt(j);
2831
2832 // Cast to a boundary node
2833 BoundaryNodeBase* bnod_pt =
2834 dynamic_cast<BoundaryNodeBase*>(node_pt(j));
2835
2836 // Get the index of the first nodal value associated with
2837 // this FaceElement
2838 unsigned first_index =
2840
2841 // Assemble higher-dimensional quantities
2842 for (unsigned i = 0; i < dim_el + 1; i++)
2843 {
2844 x[i] += nodal_position(j, i) * psi(j);
2845 lambda[i] += nod_pt->value(first_index + i) * psi(j);
2846 for (unsigned ii = 0; ii < dim_el; ii++)
2847 {
2848 interpolated_a(ii, i) +=
2849 lagrangian_position(j, i) * dpsids(j, ii);
2850 }
2851 }
2852 }
2853
2854 // Now find the local undeformed metric tensor from the tangent Vectors
2855 DenseMatrix<double> a(dim_el);
2856 for (unsigned i = 0; i < dim_el; i++)
2857 {
2858 for (unsigned j = 0; j < dim_el; j++)
2859 {
2860 // Initialise surface metric tensor to zero
2861 a(i, j) = 0.0;
2862 // Take the dot product
2863 for (unsigned k = 0; k < dim_el + 1; k++)
2864 {
2865 a(i, j) += interpolated_a(i, k) * interpolated_a(j, k);
2866 }
2867 }
2868 }
2869
2870
2871 // Find the determinant of the metric tensor
2872 double adet = 0.0;
2873 switch (dim_el + 1)
2874 {
2875 case 2:
2876 adet = a(0, 0);
2877 break;
2878
2879 case 3:
2880 adet = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
2881 break;
2882
2883 default:
2884 throw OomphLibError(
2885 "Wrong dimension "
2886 "fill_in_generic_contribution_to_residuals_displ_lagr_multiplier",
2887 "FSIImposeDisplacementByLagrangeMultiplierElement::fill_in_"
2888 "generic_contribution_to_residuals_fsi_displ_lagr_multiplier()",
2889 OOMPH_EXCEPTION_LOCATION);
2890 }
2891
2892 // Get prescribed wall shape from adjacent bulk element
2893 Vector<double> r_prescribed(dim_el + 1);
2894
2895 // Get the local coordinate in the solid element (copy
2896 // operation for Vector)
2897 Vector<double> s_adjacent(external_element_local_coord(0, ipt));
2898
2899 // Get the position in the adjacent element
2900 FiniteElement* bulk_el_pt = external_element_pt(0, ipt);
2901 bulk_el_pt->interpolated_x(s_adjacent, r_prescribed);
2902
2903 // Premultiply the weights and the square-root of the determinant of
2904 // the metric tensor
2905 double W = w * sqrt(adet);
2906
2907 // Assemble residuals and jacobian
2908
2909 // Loop over directions
2910 for (unsigned i = 0; i < dim_el + 1; i++)
2911 {
2912 // Loop over the nodes
2913 for (unsigned j = 0; j < n_node; j++)
2914 {
2915 // Assemble residual for Lagrange multiplier:
2916
2917 // Cast to a boundary node
2918 BoundaryNodeBase* bnod_pt =
2919 dynamic_cast<BoundaryNodeBase*>(node_pt(j));
2920
2921 // Local eqn number:
2922 int local_eqn = nodal_local_eqn(
2923 j,
2925
2926
2927 if (local_eqn >= 0)
2928 {
2929 residuals[local_eqn] += (x[i] - r_prescribed[i]) * psi(j) * W;
2930
2931 // Do Jacobian too?
2932 if (flag == 1)
2933 {
2934 // Loop over the nodes again for unknowns (only diagonal
2935 // contribution to direction!).
2936 for (unsigned jj = 0; jj < n_node; jj++)
2937 {
2938 int local_unknown = position_local_eqn(jj, 0, i);
2939 if (local_unknown >= 0)
2940 {
2941 jacobian(local_eqn, local_unknown) += psi(jj) * psi(j) * W;
2942 }
2943 }
2944 }
2945 }
2946
2947
2948 // Add Lagrange multiplier contribution to bulk equations
2949
2950 // Local eqn number: Node, type, direction
2951 local_eqn = position_local_eqn(j, 0, i);
2952 if (local_eqn >= 0)
2953 {
2954 // Add to residual
2955 residuals[local_eqn] += lambda[i] * psi(j) * W;
2956
2957 // Do Jacobian too?
2958 if (flag == 1)
2959 {
2960 // Loop over the nodes again for unknowns (only diagonal
2961 // contribution to direction!).
2962 for (unsigned jj = 0; jj < n_node; jj++)
2963 {
2964 // Cast to a boundary node
2965 BoundaryNodeBase* bnode_pt =
2966 dynamic_cast<BoundaryNodeBase*>(node_pt(jj));
2967
2968 int local_unknown = nodal_local_eqn(
2969 jj,
2971 Id) +
2972 i);
2973 if (local_unknown >= 0)
2974 {
2975 jacobian(local_eqn, local_unknown) += psi(jj) * psi(j) * W;
2976 }
2977 }
2978 }
2979 }
2980 }
2981 }
2982
2983
2984 } // End of loop over the integration points
2985 }
2986
2987
2988 /// The number of "DOF types" that degrees of freedom in this element
2989 /// are sub-divided into: Just the solid degrees of freedom themselves.
2990 unsigned ndof_types() const
2991 {
2992 return this->dim() + 1;
2993 };
2994
2995
2996 /// Create a list of pairs for all unknowns in this element,
2997 /// so that the first entry in each pair contains the global equation
2998 /// number of the unknown, while the second one contains the number
2999 /// of the dof that this unknown is associated with.
3000 /// (Function can obviously only be called if the equation numbering
3001 /// scheme has been set up.)
3003 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list) const
3004 {
3005 // temporary pair (used to store dof lookup prior to being added to list)
3006 std::pair<unsigned, unsigned> dof_lookup;
3007
3008 // number of nodes
3009 const unsigned n_node = this->nnode();
3010
3011 // Loop over directions
3012 unsigned dim_el = this->dim();
3013 for (unsigned i = 0; i < dim_el + 1; i++)
3014 {
3015 // Loop over the nodes
3016 for (unsigned j = 0; j < n_node; j++)
3017 {
3018 // Cast to a boundary node
3019 BoundaryNodeBase* bnod_pt =
3020 dynamic_cast<BoundaryNodeBase*>(node_pt(j));
3021
3022 // Local eqn number:
3023 int local_eqn = nodal_local_eqn(
3025 if (local_eqn >= 0)
3026 {
3027 // store dof lookup in temporary pair: First entry in pair
3028 // is global equation number; second entry is dof type
3029 dof_lookup.first = this->eqn_number(local_eqn);
3030 dof_lookup.second = i;
3031
3032 // add to list
3033 dof_lookup_list.push_front(dof_lookup);
3034 }
3035 }
3036 }
3037 }
3038
3039 /// Lagrange Id
3040 unsigned Id;
3041 };
3042
3043
3044 /// /////////////////////////////////////////////////////////////////////////
3045 /// /////////////////////////////////////////////////////////////////////////
3046 /// /////////////////////////////////////////////////////////////////////////
3047
3048
3049 //======================================================================
3050 /// A class for elements that allow the imposition of a displacement
3051 /// constraint for bulk solid elements via a Lagrange multiplier.
3052 /// Prescribed displaced is obtained from an adjacent bulk solid
3053 /// element (rather than from a lower-dimensional GeomObject
3054 /// as in the corresponding ImposeDisplacementByLagrangeMultiplierElement
3055 /// class. The present class is particularly suited for parallel
3056 /// FSI computations.
3057 /// \b NOTE: Currently (and for the foreseeable future) this
3058 /// element only works with bulk elements that do not have
3059 /// generalised degrees of freedom (so it won't work with
3060 /// Hermite-type elements, say). The additional functionality
3061 /// to deal with such elements could easily be added (once a
3062 /// a suitable test case is written). For now we simply throw
3063 /// errors if an attempt is made to use the element with an unsuitable
3064 /// bulk element.
3065 ///
3066 /// REFINEABLE VERSION
3067 //======================================================================
3068 template<class ELEMENT>
3070 : public virtual FSIImposeDisplacementByLagrangeMultiplierElement<ELEMENT>,
3072
3073 {
3074 public:
3075 /// Function to describe the local dofs of the element. The ostream
3076 /// specifies the output stream to which the description
3077 /// is written; the string stores the currently
3078 /// assembled output that is ultimately written to the
3079 /// output stream by Data::describe_dofs(...); it is typically
3080 /// built up incrementally as we descend through the
3081 /// call hierarchy of this function when called from
3082 /// Problem::describe_dofs(...)
3084 ELEMENT>::describe_local_dofs;
3085
3086 /// Constructor takes a "bulk" element and the
3087 /// index that identifies which face the FaceElement is supposed
3088 /// to be attached to. The optional identifier can be used
3089 /// to distinguish the additional nodal values created by
3090 /// this element from thos created by other FaceElements.
3092 FiniteElement* const& element_pt,
3093 const int& face_index,
3094 const unsigned& id = 0)
3096 element_pt, face_index, id, true)
3097 {
3098 }
3099
3100
3101 /// Number of continuously interpolated values: Same for
3102 /// all nodes since it's a refineable element
3104 {
3105 return node_pt(0)->nvalue();
3106 }
3107
3108 /// Fill in the residuals
3110 {
3111 // Call the generic routine with the flag set to 0
3113 residuals, GeneralisedElement::Dummy_matrix, 0);
3114 }
3115
3116
3117 /// Fill in contribution from Jacobian
3119 DenseMatrix<double>& jacobian)
3120 {
3121 // Call the generic routine with the flag set to 1
3123 residuals, jacobian, 1);
3124
3125 // Add the contribution of the external interaction by finite differences
3127 }
3128
3129
3130 protected:
3131 /// Helper function to compute the residuals and, if flag==1, the
3132 /// Jacobian
3134 Vector<double>& residuals,
3135 DenseMatrix<double>& jacobian,
3136 const unsigned& flag)
3137 {
3138 // Find out how many positional dofs there are
3139 unsigned n_position_type = this->nnodal_position_type();
3140
3141#ifdef PARANOID
3142 if (n_position_type != 1)
3143 {
3144 throw OomphLibError(
3145 "RefineableImposeDisplacementByLagrangeMultiplierElement \n cannot "
3146 "(currently) be used with elements that have generalised\n "
3147 "positional dofs. Upgrade should be straightforward though the code "
3148 "is in a \n bit of state, with generalised degrees of freedom "
3149 "sometimes half taken into \n account, sometimes completely "
3150 "ignored...",
3151 OOMPH_CURRENT_FUNCTION,
3152 OOMPH_EXCEPTION_LOCATION);
3153 }
3154#endif
3155
3156 // Find out how many nodes there are
3157 unsigned n_node = this->nnode();
3158
3159 // Dimension of element
3160 unsigned dim_el = this->dim();
3161
3162 // Set up memory for the shape functions
3163 Shape psi(n_node);
3164 DShape dpsids(n_node, dim_el);
3165
3166 // Set the value of n_intpt
3167 unsigned n_intpt = this->integral_pt()->nweight();
3168
3169
3170 // Integers to store local equation number
3171 int local_eqn = 0;
3172 int local_unknown = 0;
3173
3174 // Loop over the integration points
3175 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
3176 {
3177 // Get the integral weight
3178 double w = this->integral_pt()->weight(ipt);
3179
3180 // Only need to call the local derivatives
3181 this->dshape_local_at_knot(ipt, psi, dpsids);
3182
3183 // Calculate the Eulerian coordinates and Lagrange multiplier
3184 Vector<double> x(dim_el + 1, 0.0);
3185 Vector<double> lambda(dim_el + 1, 0.0);
3186 Vector<double> zeta(dim_el, 0.0);
3187 DenseMatrix<double> interpolated_a(dim_el, dim_el + 1, 0.0);
3188
3189 // Loop over nodes -- note in these calculations hang-ness is
3190 // automatically taken into account because of calls to position(...)
3191 // etc
3192 for (unsigned j = 0; j < n_node; j++)
3193 {
3194 Node* nod_pt = this->node_pt(j);
3195
3196 // Cast to a boundary node
3197 BoundaryNodeBase* bnod_pt = dynamic_cast<BoundaryNodeBase*>(nod_pt);
3198
3199 // Get the index of the first nodal value associated with
3200 // this FaceElement
3201 unsigned first_index =
3203
3204 // Assemble higher-dimensional quantities
3205 for (unsigned i = 0; i < dim_el + 1; i++)
3206 {
3207 x[i] += this->nodal_position(j, i) * psi(j);
3208 lambda[i] += nod_pt->value(first_index + i) * psi(j);
3209 for (unsigned ii = 0; ii < dim_el; ii++)
3210 {
3211 interpolated_a(ii, i) +=
3212 this->lagrangian_position(j, i) * dpsids(j, ii);
3213 }
3214 }
3215 }
3216
3217
3218 // Now find the local undeformed metric tensor from the tangent Vectors
3219 DenseMatrix<double> a(dim_el);
3220 for (unsigned i = 0; i < dim_el; i++)
3221 {
3222 for (unsigned j = 0; j < dim_el; j++)
3223 {
3224 // Initialise surface metric tensor to zero
3225 a(i, j) = 0.0;
3226 // Take the dot product
3227 for (unsigned k = 0; k < dim_el + 1; k++)
3228 {
3229 a(i, j) += interpolated_a(i, k) * interpolated_a(j, k);
3230 }
3231 }
3232 }
3233
3234
3235 // Find the determinant of the metric tensor
3236 double adet = 0.0;
3237 switch (dim_el + 1)
3238 {
3239 case 2:
3240 adet = a(0, 0);
3241 break;
3242
3243 case 3:
3244 adet = a(0, 0) * a(1, 1) - a(0, 1) * a(1, 0);
3245 break;
3246
3247 default:
3248 throw OomphLibError(
3249 "Wrong dimension "
3250 "refineable_fill_in_generic_contribution_to_residuals_displ_lagr_"
3251 "multiplier",
3252 "RefineableImposeDisplacementByLagrangeMultiplierElement::"
3253 "refineablefill_in_generic_contribution_to_residuals_displ_lagr_"
3254 "multiplier()",
3255 OOMPH_EXCEPTION_LOCATION);
3256 }
3257
3258 // Get prescribed wall shape from adjacent bulk element
3259 Vector<double> r_prescribed(dim_el + 1);
3260
3261 // Get the local coordinate in the solid element (copy
3262 // operation for Vector)
3263 Vector<double> s_adjacent(this->external_element_local_coord(0, ipt));
3264
3265 // Get the position in the adjacent element
3266 FiniteElement* bulk_el_pt = this->external_element_pt(0, ipt);
3267 bulk_el_pt->interpolated_x(s_adjacent, r_prescribed);
3268
3269
3270 // Premultiply the weights and the square-root of the determinant of
3271 // the metric tensor
3272 double W = w * sqrt(adet);
3273
3274 // Assemble residuals and jacobian
3275
3276
3277 // Number of master nodes and storage for the weight of the shape
3278 // function
3279 unsigned n_master = 1;
3280 unsigned n_master2 = 1;
3281 double hang_weight = 1.0;
3282 double hang_weight2 = 1.0;
3283
3284 // Pointer to hang info object
3285 HangInfo* hang_info_pt = 0;
3286 HangInfo* hang_info2_pt = 0;
3287
3288
3289 // Loop over nodes
3290 for (unsigned j = 0; j < n_node; j++)
3291 {
3292 // Local node itself (hanging or not)
3293 Node* local_node_pt = this->node_pt(j);
3294
3295 // Local boolean to indicate whether the node is hanging
3296 bool is_node_hanging = local_node_pt->is_hanging();
3297
3298 // If the node is hanging
3299 if (is_node_hanging)
3300 {
3301 hang_info_pt = local_node_pt->hanging_pt();
3302
3303 // Read out number of master nodes from hanging data
3304 n_master = hang_info_pt->nmaster();
3305 }
3306 // Otherwise the node is its own master
3307 else
3308 {
3309 n_master = 1;
3310 }
3311
3312 // Loop over the master nodes
3313 for (unsigned m = 0; m < n_master; m++)
3314 {
3315 // Loop over velocity components for equations
3316 for (unsigned i = 0; i < dim_el + 1; i++)
3317 {
3318 // Get the equation number for Lagrange multiplier eqn
3319
3320 // If the node is hanging
3321 if (is_node_hanging)
3322 {
3323 // Cast to a boundary node
3324 BoundaryNodeBase* bnod_pt = dynamic_cast<BoundaryNodeBase*>(
3325 hang_info_pt->master_node_pt(m));
3326
3327 // Get the equation number from the master node
3328 local_eqn = this->local_hang_eqn(
3329 hang_info_pt->master_node_pt(m),
3331 this->Id) +
3332 i);
3333
3334 // Get the hang weight from the master node
3335 hang_weight = hang_info_pt->master_weight(m);
3336 }
3337 // If the node is not hanging
3338 else
3339 {
3340 // Cast to a boundary node
3341 BoundaryNodeBase* bnod_pt =
3342 dynamic_cast<BoundaryNodeBase*>(local_node_pt);
3343
3344 // Local equation number
3345 local_eqn = this->nodal_local_eqn(
3346 j,
3348 this->Id) +
3349 i);
3350
3351 // Node contributes with full weight
3352 hang_weight = 1.0;
3353 }
3354
3355 // If it's not a boundary condition...
3356 if (local_eqn >= 0)
3357 {
3358 residuals[local_eqn] +=
3359 (x[i] - r_prescribed[i]) * psi(j) * W * hang_weight;
3360
3361 // Do Jacobian too?
3362 if (flag == 1)
3363 {
3364 // Loop over the nodes again for unknowns (only diagonal
3365 // contribution to direction!).
3366 for (unsigned jj = 0; jj < n_node; jj++)
3367 {
3368 // Local node itself (hanging or not)
3369 Node* local_node2_pt = this->node_pt(jj);
3370
3371 // Local boolean to indicate whether the node is hanging
3372 bool is_node_hanging2 = local_node2_pt->is_hanging();
3373
3374 // If the node is hanging
3375 if (is_node_hanging2)
3376 {
3377 hang_info2_pt = local_node2_pt->hanging_pt();
3378
3379 // Read out number of master nodes from hanging data
3380 n_master2 = hang_info2_pt->nmaster();
3381 }
3382 // Otherwise the node is its own master
3383 else
3384 {
3385 n_master2 = 1;
3386 }
3387
3388 // Loop over the master nodes
3389 for (unsigned m2 = 0; m2 < n_master2; m2++)
3390 {
3391 // Storage for local equation numbers at node indexed by
3392 // type and direction
3393 DenseMatrix<int> position_local_eqn_at_node2(
3394 n_position_type, dim_el + 1);
3395
3396 if (is_node_hanging2)
3397 {
3398 // Find the equation numbers
3399 position_local_eqn_at_node2 = local_position_hang_eqn(
3400 local_node2_pt->hanging_pt()->master_node_pt(m2));
3401
3402 // Find the hanging node weight
3403 hang_weight2 =
3404 local_node2_pt->hanging_pt()->master_weight(m2);
3405 }
3406 else
3407 {
3408 // Non-loop of types of dofs
3409 // for(unsigned k2=0;k2<n_position_type;k2++)
3410 // {
3411 unsigned k2 = 0;
3412
3413 // Loop over the displacement components
3414 // for(unsigned i2=0;i2<dim_el+1;i2++)
3415 unsigned i2 =
3416 i; // only need that one, but need to store
3417 // information in this container because
3418 // it's required for hanging case.
3419 {
3420 position_local_eqn_at_node2(k2, i2) =
3421 this->position_local_eqn(jj, k2, i2);
3422 }
3423
3424 // Hang weight is one
3425 hang_weight2 = 1.0;
3426 }
3427
3428 unsigned k2 = 0;
3429 local_unknown = position_local_eqn_at_node2(k2, i);
3430 if (local_unknown >= 0)
3431 {
3432 jacobian(local_eqn, local_unknown) +=
3433 psi(jj) * psi(j) * W * hang_weight * hang_weight2;
3434 }
3435 }
3436 }
3437 }
3438
3439
3440 // Add Lagrange multiplier contribution to bulk equations
3441
3442
3443 // Storage for local equation numbers at node indexed by
3444 // type and direction
3445 DenseMatrix<int> position_local_eqn_at_node(n_position_type,
3446 dim_el + 1);
3447
3448 if (is_node_hanging)
3449 {
3450 // Find the equation numbers
3451 position_local_eqn_at_node = local_position_hang_eqn(
3452 local_node_pt->hanging_pt()->master_node_pt(m));
3453 }
3454 else
3455 {
3456 // Non-loop of types of dofs
3457 // for(unsigned k2=0;k2<n_position_type;k2++)
3458 // {
3459 unsigned k2 = 0;
3460
3461 // Loop over the displacement components
3462 // for(unsigned i2=0;i2<dim_el+1;i2++)
3463 unsigned i2 = i; // only need that one, but need to store
3464 // information in this container because
3465 // it's required for hanging case.
3466 {
3467 position_local_eqn_at_node(k2, i2) =
3468 this->position_local_eqn(j, k2, i2);
3469 }
3470 }
3471 unsigned k = 0;
3472 local_eqn = position_local_eqn_at_node(k, i);
3473
3474 /*IF it's not a boundary condition*/
3475 if (local_eqn >= 0)
3476 {
3477 // Add to residual
3478 residuals[local_eqn] += lambda[i] * psi(j) * W * hang_weight;
3479
3480 // Do Jacobian too?
3481 if (flag == 1)
3482 {
3483 // Loop over the nodes again for unknowns (only diagonal
3484 // contribution to direction!).
3485 for (unsigned jj = 0; jj < n_node; jj++)
3486 {
3487 // Local node itself (hanging or not)
3488 Node* local_node2_pt = this->node_pt(jj);
3489
3490 // Local boolean to indicate whether the node is hanging
3491 bool is_node_hanging2 = local_node2_pt->is_hanging();
3492
3493 // If the node is hanging
3494 if (is_node_hanging2)
3495 {
3496 hang_info2_pt = local_node2_pt->hanging_pt();
3497
3498 // Read out number of master nodes from hanging data
3499 n_master2 = hang_info2_pt->nmaster();
3500 }
3501 // Otherwise the node is its own master
3502 else
3503 {
3504 n_master2 = 1;
3505 }
3506
3507 // Loop over the master nodes
3508 for (unsigned m2 = 0; m2 < n_master2; m2++)
3509 {
3510 // Get the equation number for Lagrange multiplier eqn
3511
3512 // If the node is hanging
3513 if (is_node_hanging2)
3514 {
3515 // Cast to a boundary node
3516 BoundaryNodeBase* bnod2_pt =
3517 dynamic_cast<BoundaryNodeBase*>(
3518 hang_info2_pt->master_node_pt(m2));
3519
3520 // Get the equation number from the master node
3521 local_unknown = this->local_hang_eqn(
3522 hang_info2_pt->master_node_pt(m2),
3523 bnod2_pt
3525 this->Id) +
3526 i);
3527
3528 // Get the hang weight from the master node
3529 hang_weight2 = hang_info2_pt->master_weight(m2);
3530 }
3531 // If the node is not hanging
3532 else
3533 {
3534 // Cast to a boundary node
3535 BoundaryNodeBase* bnod2_pt =
3536 dynamic_cast<BoundaryNodeBase*>(local_node2_pt);
3537
3538 // Local equation number
3539 local_unknown = this->nodal_local_eqn(
3540 jj,
3541 bnod2_pt
3542 ->index_of_first_value_assigned_by_face_element(
3543 this->Id) +
3544 i);
3545
3546 // Node contributes with full weight
3547 hang_weight2 = 1.0;
3548 }
3549
3550 if (local_unknown >= 0)
3551 {
3552 jacobian(local_eqn, local_unknown) +=
3553 psi(jj) * psi(j) * W * hang_weight * hang_weight2;
3554 }
3555 }
3556 }
3557 }
3558 }
3559 }
3560 }
3561 }
3562 }
3563
3564 } // End of loop over the integration points
3565 }
3566 };
3567
3568
3569 /// /////////////////////////////////////////////////////////////////////////
3570 /// /////////////////////////////////////////////////////////////////////////
3571 /// /////////////////////////////////////////////////////////////////////////
3572
3573
3574} // namespace oomph
3575
3576#endif
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A class that contains the information required by Nodes that are located on Mesh boundaries....
Definition: nodes.h:1996
unsigned index_of_first_value_assigned_by_face_element(const unsigned &face_id=0) const
Return the index of the first value associated with the i-th face element value. If no argument is sp...
Definition: nodes.h:2061
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:483
This is a base class for all elements that require external sources (e.g. FSI, multi-domain problems ...
Vector< double > & external_element_local_coord(const unsigned &interaction_index, const unsigned &ipt)
Access function to get source element's local coords for specified interaction index at specified int...
void set_ninteraction(const unsigned &n_interaction)
Set the number of interactions in the element This function is usually called in the specific element...
void describe_local_dofs(std::ostream &out, const std::string &curr_string) const
Function to describe the local dofs of the element. The ostream specifies the output stream to which ...
void fill_in_jacobian_from_external_interaction_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Calculate the contributions to the jacobian from all external interaction degrees of freedom (geometr...
FiniteElement *& external_element_pt(const unsigned &interaction_index, const unsigned &ipt)
Access function to source element for specified interaction index at specified integration point.
///////////////////////////////////////////////////////////////////////// ///////////////////////////...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
FSIImposeDisplacementByLagrangeMultiplierElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0, const bool &called_from_refineable_constructor=false)
Constructor takes a "bulk" element and the index that identifies which face the FaceElement is suppos...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Fill in contribution to Mass matrix and Jacobian. There is no contributiont to mass matrix so simply ...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals.
void fill_in_generic_contribution_to_residuals_fsi_displ_lagr_multiplier(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Helper function to compute the residuals and, if flag==1, the Jacobian.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Just the soli...
void describe_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the elements. The ostream specifies the output stream to which...
void output(std::ostream &outfile)
Output function.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
void set_normal_pointing_into_fluid()
Set the normal computed by underlying FaceElement to point into the fluid.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: Just the soli...
void dposition_dlagrangian_at_local_coordinate(const Vector< double > &s, DenseMatrix< double > &drdxi) const
Derivative of position vector w.r.t. the SolidFiniteElement's Lagrangian coordinates; evaluated at cu...
FSISolidTractionElement(FiniteElement *const &element_pt, const int &face_index, const bool &called_from_refineable_constructor=false)
Constructor: Create element as FSIWallElement (and thus, by inheritance, a GeomObject) with DIM-1 Lag...
virtual void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Final overload... Forwards to the version in the FSIWallElement.
virtual void get_traction(const unsigned &intpt, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Get the load vector: Pass number of the integration point, Lagr. coordinate, Eulerian coordinate (nei...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function: Note we can only output the traction at Gauss points so n_plot is actually ignored.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
void set_normal_pointing_out_of_fluid()
Set the normal computed by underlying FaceElement to point out of the fluid.
bool Normal_points_into_fluid
Boolean flag to indicate whether the normal is directed into the fluid.
virtual void(*&)(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) traction_fct_pt()
Broken overloaded reference to the traction function pointer. It doesn't make sense to specify an ext...
////////////////////////////////////////////////////////////////////////// //////////////////////////...
Definition: fsi.h:194
void setup_fsi_wall_element(const unsigned &nlagr_solid, const unsigned &ndim_fluid)
Setup: Assign storage – pass the Eulerian dimension of the "adjacent" fluid elements and the number o...
Definition: fsi.cc:121
void fluid_load_vector(const unsigned &intpt, const Vector< double > &N, Vector< double > &load)
Get FE Jacobian by systematic finite differencing w.r.t. nodal positition Data, internal and external...
Definition: fsi.cc:162
FaceElements are elements that coincide with the faces of higher-dimensional "bulk" elements....
Definition: elements.h:4338
bool Boundary_number_in_bulk_mesh_has_been_set
Has the Boundary_number_in_bulk_mesh been set? Only included if compiled with PARANOID switched on.
Definition: elements.h:4394
unsigned Boundary_number_in_bulk_mesh
The boundary number in the bulk mesh to which this element is attached.
Definition: elements.h:4388
int & face_index()
Index of the face (a number that uniquely identifies the face in the element)
Definition: elements.h:4626
const unsigned & boundary_number_in_bulk_mesh() const
Broken assignment operator.
Definition: elements.h:4475
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:6006
void add_additional_values(const Vector< unsigned > &nadditional_values, const unsigned &id)
Helper function adding additional values for the unknowns associated with the FaceElement....
Definition: elements.h:4428
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
FiniteElement *& bulk_element_pt()
Pointer to higher-dimensional "bulk" element.
Definition: elements.h:4735
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:4998
A general Finite Element class.
Definition: elements.h:1313
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual void dshape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids) const
Return the geometric shape function and its derivative w.r.t. the local coordinates at the ipt-th int...
Definition: elements.cc:3239
virtual void output(std::ostream &outfile)
Output the element data — typically the values at the nodes in a format suitable for post-processing.
Definition: elements.h:3050
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
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
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
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
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement....
Definition: elements.cc:5132
virtual void describe_nodal_local_dofs(std::ostream &out, const std::string &current_string) const
Function to describe the local dofs of the element[s]. The ostream specifies the output stream to whi...
Definition: elements.cc:1727
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
virtual void shape_at_knot(const unsigned &ipt, Shape &psi) const
Return the geometric shape function at the ipt-th integration point.
Definition: elements.cc:3220
bool has_hanging_nodes() const
Return boolean to indicate if any of the element's nodes are geometrically hanging.
Definition: elements.h:2470
void fill_in_jacobian_from_external_by_fd(Vector< double > &residuals, DenseMatrix< double > &jacobian, const bool &fd_all_data=false)
Calculate the contributions to the jacobian from the external degrees of freedom using finite differe...
Definition: elements.cc:1199
unsigned nexternal_data() const
Return the number of external data objects.
Definition: elements.h:829
bool is_halo() const
Is this element a halo?
Definition: elements.h:1163
unsigned ndof() const
Return the number of equations/dofs in the element.
Definition: elements.h:835
unsigned long eqn_number(const unsigned &ieqn_local) const
Return the global equation number corresponding to the ieqn_local-th local equation number.
Definition: elements.h:704
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
void flush_external_data()
Flush all external data.
Definition: elements.cc:387
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) external data object to the element and return its index (i....
Definition: elements.cc:307
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
virtual unsigned ngeom_data() const
How many items of Data does the shape of the object depend on? This is implemented as a broken virtua...
Definition: geom_objects.h:209
virtual Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on. This is implemented as a bro...
Definition: geom_objects.h:233
virtual void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
A geometric object may be composed of may sub-objects (e.g. a finite-element representation of a boun...
Definition: geom_objects.h:381
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
///////////////////////////////////////////////////////////////////////// ///////////////////////////...
Vector< Vector< double > > Zeta_sub_geom_object
Storage for local coordinates in sub-GeomObjects at integration points.
unsigned N_assigned_geom_data
Bool to record the number of geom Data that has been assigned to external data (we're keeping a recor...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals.
Vector< GeomObject * > Sub_geom_object_pt
Storage for sub-GeomObject at the integration points.
unsigned ndof_types() const
The number of "DOF types" that degrees of freedom in this element are sub-divided into: We only label...
void output(std::ostream &outfile)
Output function.
void set_boundary_shape_geom_object_pt(GeomObject *boundary_shape_geom_object_pt, const unsigned &boundary_number_in_bulk_mesh)
Set GeomObject that specifies the prescribed boundary displacement; GeomObject is assumed to be param...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
double square_of_l2_norm_of_error()
Compute square of L2 norm of error between prescribed and actual displacement.
GeomObject * Boundary_shape_geom_object_pt
GeomObject that specifies the prescribed boundary displacement; GeomObject is assumed to be parametri...
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned > > &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
GeomObject * boundary_shape_geom_object_pt() const
Access to GeomObject that specifies the prescribed boundary displacement; GeomObject is assumed to be...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
bool Sparsify
Boolean flag to indicate if we're using geometric Data of sub-GeomObjects that make up the (possibly ...
ImposeDisplacementByLagrangeMultiplierElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0, const bool &called_from_refineable_constructor=false)
Constructor takes a "bulk" element and the index that identifies which face the FaceElement is suppos...
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
Fill in contribution to Mass matrix and Jacobian. There is no contributiont to mass matrix so simply ...
void fill_in_generic_contribution_to_residuals_displ_lagr_multiplier(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Helper function to compute the residuals and, if flag==1, the Jacobian.
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.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1285
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation....
Definition: nodes.cc:2408
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1228
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
An OomphLibError object which should be thrown when an run-time error is encountered....
RefineableElements are FiniteElements that may be subdivided into children to provide a better local ...
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
///////////////////////////////////////////////////////////////////////// ///////////////////////////...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
RefineableFSIImposeDisplacementByLagrangeMultiplierElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0)
Constructor takes a "bulk" element and the index that identifies which face the FaceElement is suppos...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals.
void refineable_fill_in_generic_contribution_to_residuals_fsi_displ_lagr_multiplier(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Helper function to compute the residuals and, if flag==1, the Jacobian.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: Same for all nodes since it's a refineable element.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
virtual void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Final overload. Get contributions from refineable solid traction element and derivatives from externa...
RefineableFSISolidTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor: Create element as FSIWallElement (and thus, by inheritance, a GeomObject) with DIM-1 Lag...
///////////////////////////////////////////////////////////////////////// ///////////////////////////...
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals.
void refineable_fill_in_generic_contribution_to_residuals_displ_lagr_multiplier(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Helper function to compute the residuals and, if flag==1, the Jacobian.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values: Same for all nodes since it's a refineable element.
RefineableImposeDisplacementByLagrangeMultiplierElement(FiniteElement *const &element_pt, const int &face_index, const unsigned &id=0)
Constructor takes a "bulk" element and the index that identifies which face the FaceElement is suppos...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
DenseMatrix< int > & local_position_hang_eqn(Node *const &node_pt)
Access the local equation number of of hanging node variables associated with nodal positions....
////////////////////////////////////////////////////////////////////// //////////////////////////////...
~RefineableSolidTractionElement()
Destructor should not delete anything.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
This function returns just the residuals.
RefineableSolidTractionElement(FiniteElement *const &element_pt, const int &face_index)
Constructor, which takes a "bulk" element and the face index.
unsigned ncont_interpolated_values() const
Number of continuously interpolated values are the same as those in the bulk element.
void refineable_fill_in_contribution_to_residuals_solid_traction(Vector< double > &residuals)
This function returns the residuals for the traction function.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
This function returns the residuals and the Jacobian.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
SolidFaceElements combine FaceElements and SolidFiniteElements and overload various functions so they...
Definition: elements.h:4914
double zeta_nodal(const unsigned &n, const unsigned &k, const unsigned &i) const
The "global" intrinsic coordinate of the element when viewed as part of a geometric object should be ...
Definition: elements.h:4921
double interpolated_xi(const Vector< double > &s, const unsigned &i) const
Return i-th FE-interpolated Lagrangian coordinate xi[i] at local coordinate s. Overloaded from SolidF...
Definition: elements.h:4937
Data * geom_data_pt(const unsigned &j)
Return pointer to the j-th Data item that the object's shape depends on. (Redirects to the nodes' pos...
Definition: elements.h:3615
double lagrangian_position(const unsigned &n, const unsigned &i) const
Return i-th Lagrangian coordinate at local node n.
Definition: elements.h:3905
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
A class for elements that allow the imposition of an applied traction in the principle of virtual dis...
void output(std::ostream &outfile, const unsigned &n_plot)
Output function.
void fill_in_contribution_to_residuals_solid_traction(Vector< double > &residuals)
Helper function that actually calculates the residuals.
void(*&)(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction) traction_fct_pt()
Reference to the traction function pointer.
SolidTractionElement(FiniteElement *const &element_pt, const int &face_index, const bool &called_from_refineable_constructor=false)
Constructor, which takes a "bulk" element and the value of the index and its limit.
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in contribution from Jacobian.
void output(std::ostream &outfile)
Output function.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Return the residuals.
void traction(const Vector< double > &s, Vector< double > &traction)
Compute traction vector at specified local coordinate Should only be used for post-processing; ignore...
virtual void get_traction(const unsigned &intpt, const Vector< double > &xi, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Get the traction vector: Pass number of integration point (dummy), Lagr. coordinate and normal vector...
void output(FILE *file_pt)
C_style output function.
void output(FILE *file_pt, const unsigned &n_plot)
C-style output function.
void(* Traction_fct_pt)(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &n, Vector< double > &result)
Pointer to an imposed traction function. Arguments: Lagrangian coordinate; Eulerian coordinate; outer...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void Zero_traction_fct(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Default load function (zero traction)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...