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