foeppl_von_karman_volume_constraint_element.h

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-2022 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//====================================================================
// Header file for FoepplvonKarman elements
#ifndef OOMPH_FVK_VOLUME_CONSTRAINT_ELEMENT_HEADER
#define OOMPH_FVK_VOLUME_CONSTRAINT_ELEMENT_HEADER
 
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "../generic/elements.h"
#include "../meshes/triangle_mesh.h"
#include "Tfoeppl_von_karman_elements.h"
 
#include <fstream>
 
namespace oomph
{
  //=============================================================
  //=============================================================
  template<class ELEMENT, template<class> class MESH>
  class FoepplvonKarmanVolumeConstraintElement
    : public virtual GeneralisedElement
  {
  public:
    FoepplvonKarmanVolumeConstraintElement(
      MESH<ELEMENT>* bounding_mesh_pt,
      const Vector<unsigned>& contributing_region,
      const double& pressure = 0.0)
      : Bounding_mesh_pt(bounding_mesh_pt),
        Contributing_region(contributing_region),
        Prescribed_volume_pt(0)
    {
      // Create instance of pressure that is traded for volume constraint
      Volume_control_pressure_pt = new Data(1);
      Volume_control_pressure_pt->set_value(0, pressure);
 
      // Add the data object as internal data for this element
      Pressure_data_index = add_internal_data(Volume_control_pressure_pt);
    }
 
    // Destructor (empty)
    ~FoepplvonKarmanVolumeConstraintElement() {}
 
    FoepplvonKarmanVolumeConstraintElement(
      const FoepplvonKarmanVolumeConstraintElement&) = delete;
 
    void operator=(const FoepplvonKarmanVolumeConstraintElement&) = delete;
 
    double get_bounded_volume()
    {
      double bounded_volume = 0.0;
 
      // Loop over the bounded regions
      unsigned n_contributing_regions = Contributing_region.size();
      for (unsigned r = 0; r < n_contributing_regions; r++)
      {
        // Loop over the elements in the bounded regions
        unsigned n_inner_el =
          Bounding_mesh_pt->nregion_element(Contributing_region[r]);
        for (unsigned e = 0; e < n_inner_el; e++)
        {
          // Add the contribution
          ELEMENT* el_pt = dynamic_cast<ELEMENT*>(
            Bounding_mesh_pt->region_element_pt(Contributing_region[r], e));
          bounded_volume += el_pt->get_bounded_volume();
        }
      }
      return bounded_volume;
    }
 
    void assign_additional_local_eqn_numbers()
    {
      Volume_control_local_eqn = internal_local_eqn(Pressure_data_index, 0);
    }
 
    void fill_in_contribution_to_residuals(Vector<double>& residuals)
    {
      if (Volume_control_local_eqn >= 0)
      {
        residuals[Volume_control_local_eqn] += -(*Prescribed_volume_pt);
      }
    }
 
 
    void fill_in_contribution_to_jacobian(Vector<double>& residuals,
                                          DenseMatrix<double>& jacobian)
    {
      if (Volume_control_local_eqn >= 0)
      {
        // Only add contribution to residual; Jacobian (derivs w.r.t. to
        // this element's external data is handled by FvK elements)
        residuals[Volume_control_local_eqn] -= (*Prescribed_volume_pt);
 
 
        /*  // We are in charge... */
        /*  else */
        /*   { */
        /*    // NOTE: This is very slow but keep code alive -- can be */
        /*    // recycled if/when we ever make the Jacobians in the */
        /*    // Foeppl von Karman elements analytical. */
        /*    double t_start=TimingHelpers::timer(); */
 
        /*    // Setup lookup scheme between local and global data */
        /*    std::map<unsigned,unsigned> local_external_eqn; */
        /*    unsigned next=nexternal_data(); */
        /*    for (unsigned j=0;j<next;j++) */
        /*     { */
        /*      Data* data_pt=external_data_pt(j); */
        /*      unsigned nval=data_pt->nvalue(); */
        /*      for (unsigned k=0;k<nval;k++) */
        /*       { */
        /*        int local_eqn=external_local_eqn(j,k); */
        /*        if (local_eqn>=0) */
        /*         { */
        /*          int global_eqn=data_pt->eqn_number(k); */
        /*          local_external_eqn[global_eqn]=local_eqn; */
        /*         } */
        /*       } */
        /*     } */
 
 
        /*    double t_end=TimingHelpers::timer(); */
        /*    oomph_info << "Time for local_external_eqn setup: "  */
        /*               << t_end-t_start << std::endl; */
        /*    t_start=TimingHelpers::timer(); */
 
 
        /*    // Add initial contribution */
        /*    residuals[Volume_control_local_eqn] -= (*Prescribed_volume_pt); */
 
        /*    // Initialise total bounded volume */
        /*    double bounded_volume = 0.0; */
 
        /*    // Loop over the bounded regions */
        /*    unsigned n_contributing_regions = Contributing_region.size(); */
        /*    for(unsigned r = 0; r < n_contributing_regions; r++) */
        /*     { */
        /*      // Loop over the elements in the bounded regions */
        /*      unsigned n_inner_el =  */
        /*       Bounding_mesh_pt->nregion_element(Contributing_region[r]); */
        /*      for(unsigned e = 0; e < n_inner_el; e++) */
        /*       { */
        /*        // Get element */
        /*        ELEMENT *el_pt = dynamic_cast<ELEMENT*> */
        /*         (Bounding_mesh_pt->region_element_pt(Contributing_region[r],e));
         */
 
        /*        // Get element's contribution to bounded volume and */
        /*        // derivative w.r.t. its unknowns */
        /*        double el_bounded_volume=0.0; */
 
        /*        std::map<unsigned,double> d_bounded_volume_d_unknown; */
        /*        el_pt->fill_in_d_bounded_volume_d_unknown(el_bounded_volume,
         */
        /*                                                  d_bounded_volume_d_unknown);
         */
 
        /*        // Add contribution to bounded volume */
        /*        bounded_volume += el_bounded_volume; */
 
 
        /*        // Add contribution to Jacobian */
        /*        for (std::map<unsigned,double>::iterator it= */
        /*              d_bounded_volume_d_unknown.begin();it!= */
        /*              d_bounded_volume_d_unknown.end();it++) */
        /*         { */
        /*          unsigned global_eqn=(*it).first; */
        /*          unsigned local_eqn=local_external_eqn[global_eqn]; */
        /*          jacobian(Volume_control_local_eqn,local_eqn)+=(*it).second;
         */
        /*         } */
        /*       } */
        /*     } */
        /*    // Add contribution to residuals */
        /*    residuals[Volume_control_local_eqn]+=bounded_volume; */
 
 
        /*    t_end=TimingHelpers::timer(); */
        /*    oomph_info << "Time for second part (actual work) of
         * fill...jacobian: "  */
        /*               << t_end-t_start << std::endl; */
 
        /*   } */
      }
    }
 
    void set_prescribed_volume(double* volume_pt)
    {
      Prescribed_volume_pt = volume_pt;
    }
 
    Data* pressure_data_pt() const
    {
      return Volume_control_pressure_pt;
    }
 
  protected:
    Data* Volume_control_pressure_pt;
 
    unsigned Pressure_data_index;
 
    int Volume_control_local_eqn;
 
    MESH<ELEMENT>* Bounding_mesh_pt;
 
    Vector<unsigned> Contributing_region;
 
    double* Prescribed_volume_pt;
  };
 
} // namespace oomph
 
#endif