biharmonic_elements.cc

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// oomph-lib includes
#include "biharmonic_elements.h"
 
 
namespace oomph
{
  /// Compute element residual Vector only (if JFLAG=and/or element
  /// Jacobian matrix
  template<unsigned DIM>
  void BiharmonicEquations<DIM>::
    fill_in_generic_residual_contribution_biharmonic(
      Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned JFLAG)
  {
    // Find out how many nodes there are
    unsigned n_node = this->nnode();
 
    // Find out how many values there
    unsigned n_value = this->node_pt(0)->nvalue();
 
    // Find out how many degrees of freedom are associated with 2nd derivative
    unsigned d2_dof = this->get_d2_dof();
 
    // set up memory for shape and test functions
    Shape psi(n_node, n_value);
    DShape dpsidx(n_node, n_value, DIM);
    DShape d2psidx(n_node, n_value, d2_dof);
    // Shape psi(n_node,n_value);
    DShape dpsids(n_node, n_value, DIM);
    DShape d2psids(n_node, n_value, d2_dof);
 
    // storage of single local coordinate
    Vector<double> s(DIM);
 
    // number of integration points
    unsigned n_intpt = this->integral_pt()->nweight();
 
    // loop over integration points
    for (unsigned ipt = 0; ipt < n_intpt; ipt++)
    {
      // set local coordinates to be intergration scheme knot
      for (unsigned i = 0; i < DIM; i++)
        s[i] = this->integral_pt()->knot(ipt, i);
 
      // find weight at knot
      double w = this->integral_pt()->weight(ipt);
 
      // find shape fn and derivate fns at knot point,
      // and Jacobian of local to global mapping
      double J = this->d2shape_eulerian(s, psi, dpsidx, d2psidx);
      this->d2shape_local(s, psi, dpsids, d2psids);
 
      // premultiply weight by jacobian
      double W = w * J;
 
      // initialise storage for d2u_interpolated
      Vector<double> d2u_interpolated(DIM);
      for (unsigned i = 0; i < DIM; i++)
      {
        d2u_interpolated[i] = 0.0;
      }
 
      // loop over nodes, degrees of freedom, and dimension to calculate
      // interpolated_d2u
      for (unsigned n = 0; n < n_node; n++)
      {
        for (unsigned k = 0; k < n_value; k++)
        {
          for (unsigned d = 0; d < DIM; d++)
          {
            d2u_interpolated[d] +=
              this->node_pt(n)->value(k) * d2psidx(n, k, d);
          }
        }
      }
 
      // create vector for interpolated position and calculate
      Vector<double> interpolated_position(DIM);
      this->interpolated_x(s, interpolated_position);
 
      // evaluate source function at knot position
      double source = 0.0;
      get_source(ipt, interpolated_position, source);
 
      // loop over nodes
      for (unsigned n1 = 0; n1 < n_node; n1++)
      {
        // loop over types of degrees of freedom
        for (unsigned k1 = 0; k1 < n_value; k1++)
        {
          // determine the local equation number
          int local_eqn_number = this->nodal_local_eqn(n1, k1);
 
          // if local equation number equal to -1 then its a boundary
          // node(pinned)
          if (local_eqn_number >= 0)
          {
            // add source contribution to residual
            residuals[local_eqn_number] -= source * psi(n1, k1) * W;
 
            // loop over dimension of d2u interpolated
            for (unsigned d1 = 0; d1 < DIM; d1++)
            {
              // loop over derivatives of d2psidx
              for (unsigned d2 = 0; d2 < DIM; d2++)
              {
                // add residual contribution
                residuals[local_eqn_number] +=
                  d2u_interpolated[d1] * d2psidx(n1, k1, d2) * W;
              }
 
              // calculate the jacobian
              if (JFLAG)
              {
                // loop over nodes
                for (unsigned n2 = 0; n2 < n_node; n2++)
                {
                  // loop over types of degrees of freedom
                  for (unsigned k2 = 0; k2 < n_value; k2++)
                  {
                    // get local dof, -if -1 then node pinned
                    int local_dof_number = this->nodal_local_eqn(n2, k2);
 
                    // if local dof # = -1, then its pinned
                    if (local_dof_number >= 0)
                    {
                      // loop over derivatives
                      for (unsigned d2 = 0; d2 < DIM; d2++)
                      {
                        // add contribution to jacobian
                        jacobian(local_eqn_number, local_dof_number) +=
                          d2psidx(n1, k1, d1) * d2psidx(n2, k2, d2) * W;
                      }
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  }
 
  template class BiharmonicElement<2>;
  template class BiharmonicElement<1>;
 
} // namespace oomph