biharmonic_elements.cc
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 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32  #include <oomph-lib-config.h>
33 #endif
34 
35 //oomph-lib includes
36 #include "biharmonic_elements.h"
37 
38 
39 namespace oomph
40 {
41 
42  /// \short Compute element residual Vector only (if JFLAG=and/or element
43  /// Jacobian matrix
44 template<unsigned DIM>
47  Vector<double> &residuals,
48  DenseMatrix<double>& jacobian,
49  unsigned JFLAG)
50 {
51 
52  //Find out how many nodes there are
53  unsigned n_node = this->nnode();
54 
55  //Find out how many values there
56  unsigned n_value = this->node_pt(0)->nvalue();
57 
58  //Find out how many degrees of freedom are associated with 2nd derivative
59  unsigned d2_dof = this->get_d2_dof();
60 
61  //set up memory for shape and test functions
62  Shape psi(n_node,n_value);
63  DShape dpsidx(n_node,n_value,DIM);
64  DShape d2psidx(n_node,n_value, d2_dof);
65 // Shape psi(n_node,n_value);
66  DShape dpsids(n_node,n_value,DIM);
67  DShape d2psids(n_node,n_value, d2_dof);
68 
69  // storage of single local coordinate
70  Vector<double> s(DIM);
71 
72  //number of integration points
73  unsigned n_intpt = this->integral_pt()->nweight();
74 
75  // loop over integration points
76  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
77  {
78 
79  // set local coordinates to be intergration scheme knot
80  for (unsigned i = 0; i < DIM; i++) s[i] = this->integral_pt()->knot(ipt,i);
81 
82  // find weight at knot
83  double w = this->integral_pt()->weight(ipt);
84 
85  // find shape fn and derivate fns at knot point,
86  //and Jacobian of local to global mapping
87  double J = this->d2shape_eulerian(s,psi,dpsidx,d2psidx);
88  this->d2shape_local(s,psi,dpsids,d2psids);
89 
90  // premultiply weight by jacobian
91  double W = w*J;
92 
93  // initialise storage for d2u_interpolated
94  Vector<double> d2u_interpolated(DIM);
95  for (unsigned i = 0; i < DIM; i++)
96  {
97  d2u_interpolated[i] = 0.0;
98  }
99 
100  // loop over nodes, degrees of freedom, and dimension to calculate
101  // interpolated_d2u
102  for (unsigned n = 0; n < n_node; n++)
103  {
104  for (unsigned k = 0; k < n_value; k++)
105  {
106  for (unsigned d = 0; d < DIM; d++)
107  {
108  d2u_interpolated[d] += this->node_pt(n)->value(k)*d2psidx(n,k,d);
109  }
110  }
111  }
112 
113  // create vector for interpolated position and calculate
114  Vector<double> interpolated_position(DIM);
115  this->interpolated_x(s,interpolated_position);
116 
117  // evaluate source function at knot position
118  double source = 0.0;
119  get_source(ipt, interpolated_position, source);
120 
121  // loop over nodes
122  for (unsigned n1 = 0; n1 < n_node; n1++)
123  {
124 
125  // loop over types of degrees of freedom
126  for (unsigned k1=0; k1 < n_value; k1++)
127  {
128 
129  // determine the local equation number
130  int local_eqn_number = this->nodal_local_eqn(n1,k1);
131 
132  // if local equation number equal to -1 then its a boundary node(pinned)
133  if (local_eqn_number >= 0)
134  {
135 
136  // add source contribution to residual
137  residuals[local_eqn_number] -= source*psi(n1,k1)*W;
138 
139  // loop over dimension of d2u interpolated
140  for (unsigned d1 = 0; d1 < DIM; d1++)
141  {
142 
143  // loop over derivatives of d2psidx
144  for (unsigned d2 = 0; d2 < DIM; d2++)
145  {
146 
147  // add residual contribution
148  residuals[local_eqn_number] +=
149  d2u_interpolated[d1]*d2psidx(n1,k1,d2)*W;
150 
151  }
152 
153  //calculate the jacobian
154  if(JFLAG)
155  {
156 
157  // loop over nodes
158  for (unsigned n2 = 0; n2 < n_node; n2++)
159  {
160 
161  // loop over types of degrees of freedom
162  for (unsigned k2=0; k2 < n_value; k2++)
163  {
164 
165  // get local dof, -if -1 then node pinned
166  int local_dof_number = this->nodal_local_eqn(n2,k2);
167 
168  // if local dof # = -1, then its pinned
169  if (local_dof_number >= 0)
170  {
171 
172  // loop over derivatives
173  for (unsigned d2 = 0; d2 < DIM; d2++)
174  {
175 
176  // add contribution to jacobian
177  jacobian(local_eqn_number,local_dof_number) +=
178  d2psidx(n1,k1,d1)*d2psidx(n2,k2,d2)*W;
179 
180  }
181  }
182  }
183  }
184  }
185  }
186  }
187  }
188  }
189  }
190 }
191 
192  template class BiharmonicElement<2>;
193  template class BiharmonicElement<1>;
194 
195 }
196 
197 
198 
199 
200 
201 
202 
203 
cstr elem_len * i
Definition: cfortran.h:607
static char t char * s
Definition: cfortran.h:572
biharmonic element class
virtual void fill_in_generic_residual_contribution_biharmonic(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned JFLAG)
Compute element residual Vector only (if JFLAG=and/or element Jacobian matrix.