oomph::VorticitySmoother< ELEMENT > Class Template Reference

Smoother for vorticity in 2D. More...

#include <vorticity_smoother.h>

Public Member Functions
	VorticitySmoother (const unsigned &recovery_order)
	Constructor: Set order of recovery shape functions. More...
	VorticitySmoother (const VorticitySmoother &)=delete
	Broken copy constructor. More...
void	operator= (const VorticitySmoother &)=delete
	Broken assignment operator. More...
virtual	~VorticitySmoother ()
	Empty virtual destructor. More...
unsigned &	recovery_order ()
	Access function for order of recovery polynomials. More...
void	shape_rec (const Vector< double > &x, Vector< double > &psi_r)
Integral *	integral_rec (const bool &is_q_mesh)
void	setup_patches (Mesh *&mesh_pt, std::map< Node *, Vector< ELEMENT * > * > &adjacent_elements_pt, Vector< Node * > &vertex_node_pt)
void	get_recovered_vorticity_in_patch (const Vector< ELEMENT * > &patch_el_pt, const unsigned &num_recovery_terms, Vector< double > *&recovered_vorticity_coefficient_pt, unsigned &n_deriv)
unsigned	nrecovery_order () const
	Get the recovery order. More...
void	recover_vorticity (Mesh *mesh_pt)
	Recover vorticity from patches. More...
void	recover_vorticity (Mesh *mesh_pt, DocInfo &doc_info)

Private Attributes
unsigned	Recovery_order
	Order of recovery polynomials. More...

Detailed Description

template<class ELEMENT>
class oomph::VorticitySmoother< ELEMENT >

Smoother for vorticity in 2D.

/////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////

Constructor & Destructor Documentation

VorticitySmoother() [1/2]

template<class ELEMENT >

oomph::VorticitySmoother< ELEMENT >::VorticitySmoother ( const unsigned &  recovery_order )

inline

Constructor: Set order of recovery shape functions.

      : Recovery_order(recovery_order)
    {
    }

VorticitySmoother() [2/2]

template<class ELEMENT >

oomph::VorticitySmoother< ELEMENT >::VorticitySmoother ( const VorticitySmoother< ELEMENT > &  )

delete

Broken copy constructor.

~VorticitySmoother()

template<class ELEMENT >

virtual oomph::VorticitySmoother< ELEMENT >::~VorticitySmoother ( )

inlinevirtual

Empty virtual destructor.

{}

Member Function Documentation

get_recovered_vorticity_in_patch()

template<class ELEMENT >

void oomph::VorticitySmoother< ELEMENT >::get_recovered_vorticity_in_patch ( const Vector< ELEMENT * > &  patch_el_pt, const unsigned &  num_recovery_terms, Vector< double > *&  recovered_vorticity_coefficient_pt, unsigned &  n_deriv  )

inline

Given the vector of elements that make up a patch, compute the vector of recovered vorticity coefficients and return a pointer to it. n_deriv indicates which derivative of the vorticity is supposed to be smoothed: 0: zeroth (i.e. the vorticity itself) 1: d/dx; 2: d/dy; 3: d^2/dx^2; 4: d^2/dxdy 5: d^2/dy^2 6: d^3/dx^3, 7: d^3/dx^2dy, 8: d^3/dxdy^2, 9: d^3/dy^3, 10: du/dx, 11: du/dy, 12: dv/dx, 13: dv/dy

    {
      // Find the number of elements in the patch
      unsigned nelem = patch_el_pt.size();
 
      // Get a pointer to any element
      ELEMENT* el_pt = patch_el_pt[0];
 
#ifdef PARANOID
      // If there's at least one element
      if (nelem > 0)
      {
        // Get the number of vorticity derivatives to recover
        int n_vort_derivs = el_pt->get_maximum_order_of_vorticity_derivative();
 
        // Get the number of vorticity derivatives to recover
        int n_veloc_derivs = el_pt->get_maximum_order_of_velocity_derivative();
 
        // If we're not recovering anything, we shouldn't be here
        if (n_vort_derivs + n_veloc_derivs == -1)
        {
          // Create an ostringstream object to create an error message
          std::ostringstream error_stream;
 
          // Create the error message
          error_stream << "Not recovering anything. Change the maximum number "
                       << "of derivatives to recover.";
 
          // Throw an error
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
      } // if (nelem>0)
#endif
 
      // Find the container indices associated with n_deriv
      std::pair<unsigned, unsigned> container_id =
        el_pt->vorticity_dof_to_container_id(n_deriv);
 
      // Maximum vorticity derivative order we can recover
      unsigned max_recoverable_vort_order =
        el_pt->get_maximum_order_of_recoverable_vorticity_derivative();
 
      // Maximum velocity derivative order we can recover
      unsigned max_recoverable_veloc_order =
        el_pt->get_maximum_order_of_recoverable_velocity_derivative();
 
      // Make a counter
      unsigned counter = 0;
 
      // Calculate the case value (initialise to -1 so we know if it's set
      // later)
      int case_value = -1;
 
      // Loop over the derivatives
      for (unsigned i = 0; i < max_recoverable_vort_order + 1; i++)
      {
        // Increment by the number of partial derivatives of order i
        counter += el_pt->npartial_derivative(i);
 
        // If we've exceeded the value of n_deriv then we know which vorticity
        // derivative to recover
        if (n_deriv < counter)
        {
          // We need to recover the i-th order of derivative of the vorticity
          case_value = i;
 
          // We're done here
          break;
        }
      } // for (unsigned i=0;i<max_recoverable_order+1;i++)
 
      // If we haven't set the case value yet then we must be recovering a
      // velocity derivative
      if (case_value == -1)
      {
        // Loop over the velocity order
        for (unsigned i = 1; i < max_recoverable_veloc_order + 1; i++)
        {
          // Increment by the number of velocity partial derivatives of order i
          counter += 2 * el_pt->npartial_derivative(i);
 
          // If we've exceeded the value of n_deriv then we know which vorticity
          // derivative to recover
          if (n_deriv < counter)
          {
            // We need to recover the i-th order of derivative of the vorticity
            case_value = max_recoverable_vort_order + i;
 
            // We're done here
            break;
          }
        } // for (unsigned i=1;i<max_recoverable_veloc_order+1;i++)
      } // if (case_value==-1)
 
#ifdef PARANOID
      // Sanity check: if the case value hasn't been set then something's wrong
      if (case_value == -1)
      {
        // Create a ostringstream object to create an error message
        std::ostringstream error_message_stream;
 
        // Create an error message
        error_message_stream
          << "Case order has not been set. Something's wrong!";
 
        // Throw an error
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Create space for the recovered quantity
      double recovered_quantity = 0.0;
 
      // Create/initialise matrix for linear system
      DenseDoubleMatrix recovery_mat(
        num_recovery_terms, num_recovery_terms, 0.0);
 
      // Create/initialise vector for RHS
      Vector<double> rhs(num_recovery_terms, 0.0);
 
      // Create a new integration scheme based on the recovery order
      // in the elements. Need to find the type of the element, default
      // is to assume a quad
      bool is_q_mesh = true;
 
      // If we can dynamic cast to the TElementBase, then it's a triangle/tet
      // Note that I'm assuming that all elements are of the same geometry, but
      // if they weren't we could adapt...
      if (dynamic_cast<TElementBase*>(patch_el_pt[0]))
      {
        // We're dealing with a triangle-based mesh so change the bool value
        is_q_mesh = false;
      }
 
      // Get a pointer to the appropriate integration type
      Integral* const integ_pt = this->integral_rec(is_q_mesh);
 
      // Loop over all elements in patch to assemble linear system
      for (unsigned e = 0; e < nelem; e++)
      {
        // Get pointer to element
        ELEMENT* const el_pt = patch_el_pt[e];
 
        // Create storage for the recovery shape function values
        Vector<double> psi_r(num_recovery_terms);
 
        // Create vector to hold local coordinates
        Vector<double> s(2, 0.0);
 
        // Find the number of integration points
        unsigned n_intpt = integ_pt->nweight();
 
        // Loop over the integration points
        for (unsigned ipt = 0; ipt < n_intpt; ipt++)
        {
          // Assign values of s, the local coordinate
          for (unsigned i = 0; i < 2; i++)
          {
            // Get the i-th entry of the local coordinate
            s[i] = integ_pt->knot(ipt, i);
          }
 
          // Get the integral weight
          double w = integ_pt->weight(ipt);
 
          // Jacobian of mapping
          double J = el_pt->J_eulerian(s);
 
          // Allocate space for the global coordinates
          Vector<double> x(2, 0.0);
 
          // Interpolate the global (Eulerian) coordinate
          el_pt->interpolated_x(s, x);
 
          // Premultiply the weights and the Jacobian and the geometric
          // Jacobian weight (used in axisymmetric and spherical coordinate
          // systems) -- hierher really fct of x? probably yes, actually).
          double W = w * J * (el_pt->geometric_jacobian(x));
 
          // Recovery shape functions at global (Eulerian) coordinate
          shape_rec(x, psi_r);
 
          // Use a switch statement to decide on which function to call
          switch (case_value)
          {
            case 0:
            {
              Vector<double> vorticity(1, 0.0);
              el_pt->get_vorticity(s, vorticity);
              recovered_quantity = vorticity[0];
              break;
            }
            case 1:
            {
              el_pt->get_raw_vorticity_deriv(
                s, recovered_quantity, container_id.second);
              break;
            }
            case 2:
            {
              el_pt->get_raw_vorticity_second_deriv(
                s, recovered_quantity, container_id.second);
              break;
            }
            case 3:
            {
              el_pt->get_raw_vorticity_third_deriv(
                s, recovered_quantity, container_id.second);
              break;
            }
            case 4:
            {
              el_pt->get_raw_velocity_deriv(
                s, recovered_quantity, container_id.second);
              break;
            }
            default:
            {
              oomph_info << "Never get here." << std::endl;
              abort();
            }
          }
 
          // Add elemental RHSs and recovery matrix to global versions
          //----------------------------------------------------------
          // RHS: Loop over the nodes for the test functions
          for (unsigned l = 0; l < num_recovery_terms; l++)
          {
            // Update the RHS entry ()
            rhs[l] += recovered_quantity * psi_r[l] * W;
          }
 
          // Loop over the nodes for the test functions
          for (unsigned l = 0; l < num_recovery_terms; l++)
          {
            // Loop over the nodes for the variables
            for (unsigned l2 = 0; l2 < num_recovery_terms; l2++)
            {
              // Add contribution to recovery matrix
              recovery_mat(l, l2) += psi_r[l] * psi_r[l2] * W;
            }
          } // for (unsigned l=0;l<num_recovery_terms;l++)
        } // for (unsigned ipt=0;ipt<n_intpt;ipt++)
      } // End of loop over elements that make up patch.
 
      // Delete the integration scheme
      delete integ_pt;
 
      // Linear system is now assembled: Solve recovery system
      //------------------------------------------------------
      // LU decompose the recovery matrix
      recovery_mat.ludecompose();
 
      // Back-substitute
      recovery_mat.lubksub(rhs);
 
      // Now create a matrix to store the vorticity recovery coefficients.
      // Pointer to this matrix will be returned.
      recovered_vorticity_coefficient_pt =
        new Vector<double>(num_recovery_terms);
 
      // Loop over the number of recovered terms
      for (unsigned icoeff = 0; icoeff < num_recovery_terms; icoeff++)
      {
        // Copy the RHS value over
        (*recovered_vorticity_coefficient_pt)[icoeff] = rhs[icoeff];
      }
    } // End of get_recovered_vorticity_in_patch

References e(), i, oomph::VorticitySmoother< ELEMENT >::integral_rec(), J, oomph::Integral::knot(), oomph::DenseDoubleMatrix::lubksub(), oomph::DenseDoubleMatrix::ludecompose(), oomph::Integral::nweight(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, s, oomph::VorticitySmoother< ELEMENT >::shape_rec(), w, oomph::QuadTreeNames::W, oomph::Integral::weight(), and plotDoE::x.

Referenced by oomph::VorticitySmoother< ELEMENT >::recover_vorticity().

integral_rec()

template<class ELEMENT >

Integral* oomph::VorticitySmoother< ELEMENT >::integral_rec ( const bool &  is_q_mesh )

inline

Integation scheme associated with the recovery shape functions must be of sufficiently high order to integrate the mass matrix associated with the recovery shape functions. The argument is the dimension of the elements. The integration is performed locally over the elements, so the integration scheme does depend on the geometry of the element. The type of element is specified by the boolean which is true if elements in the patch are QElements and false if they are TElements (will need change if we ever have other element types)

Find order of recovery shape functions

    {
      // Create an ostringstream object to create a string
      std::ostringstream error_stream;
 
      //----
      // 2D:
      //----
      switch (recovery_order())
      {
        case 1:
          // Complete linear polynomial in 2D:
          if (is_q_mesh)
          {
            // Return the appropriate Gauss integration scheme
            return (new Gauss<2, 2>);
          }
          else
          {
            // Return the appropriate Gauss integration scheme
            return (new TGauss<2, 2>);
          }
          break;
 
        case 2:
          // Complete quadratic polynomial in 2D:
          if (is_q_mesh)
          {
            // Return the appropriate Gauss integration scheme
            return (new Gauss<2, 3>);
          }
          else
          {
            // Return the appropriate Gauss integration scheme
            return (new TGauss<2, 3>);
          }
          break;
 
        case 3:
          // Complete cubic polynomial in 2D:
          if (is_q_mesh)
          {
            // Return the appropriate Gauss integration scheme
            return (new Gauss<2, 4>);
          }
          else
          {
            // Return the appropriate Gauss integration scheme
            return (new TGauss<2, 4>);
          }
          break;
 
        default:
          // Create an error messaage
          error_stream << "Recovery shape functions for recovery order "
                       << recovery_order()
                       << " haven't yet been implemented for 2D" << std::endl;
 
          // Throw an error
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
      }
 
      // Dummy return (never get here)
      return 0;
    } // End of integral_rec

References OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::VorticitySmoother< ELEMENT >::recovery_order().

Referenced by oomph::VorticitySmoother< ELEMENT >::get_recovered_vorticity_in_patch().

nrecovery_order()

template<class ELEMENT >

unsigned oomph::VorticitySmoother< ELEMENT >::nrecovery_order ( ) const

inline

Get the recovery order.

    {
      // Use a switch statement
      switch (Recovery_order)
      {
        case 1:
          // Linear recovery shape functions
          //--------------------------------
          return 3; // 1, x, y
          break;
 
        case 2:
          // Quadratic recovery shape functions
          //-----------------------------------
          return 6; // 1, x, y, x^2, xy, y^2
          break;
 
        case 3:
          // Cubic recovery shape functions
          //--------------------------------
          return 10; // 1, x, y, x^2, xy, y^2, x^3, x^2 y, x y^2, y^3
          break;
 
        default:
          // Any other recovery order?
          //--------------------------
          // Use an ostringstream object to create an error message
          std::ostringstream error_stream;
 
          // Create an error message
          error_stream << "Wrong Recovery_order " << Recovery_order
                       << std::endl;
 
          // Throw an error
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
      }
    } // End of nrecovery_order

References OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, and oomph::VorticitySmoother< ELEMENT >::Recovery_order.

Referenced by oomph::VorticitySmoother< ELEMENT >::recover_vorticity().

operator=()

template<class ELEMENT >

void oomph::VorticitySmoother< ELEMENT >::operator= ( const VorticitySmoother< ELEMENT > &  )

delete

Broken assignment operator.

recover_vorticity() [1/2]

template<class ELEMENT >

void oomph::VorticitySmoother< ELEMENT >::recover_vorticity ( Mesh *  mesh_pt )

inline

Recover vorticity from patches.

    {
      // Create a DocInfo object (used as a dummy argument)
      DocInfo doc_info;
 
      // Disable any documentation
      doc_info.disable_doc();
 
      // Recover the vorticity
      recover_vorticity(mesh_pt, doc_info);
    }

References oomph::DocInfo::disable_doc().

recover_vorticity() [2/2]

template<class ELEMENT >

void oomph::VorticitySmoother< ELEMENT >::recover_vorticity ( Mesh *  mesh_pt, DocInfo &  doc_info  )

inline

Recover vorticity from patches – output intermediate steps to directory specified by DocInfo object

    {
      // Start the timer
      double t_start = TimingHelpers::timer();
 
      // Allocate space for the local coordinates
      Vector<double> s(2, 0.0);
 
      // Allocate space for the global coordinates
      Vector<double> x(2, 0.0);
 
      // Make patches
      //-------------
      // Allocate space for the mapping from nodes to elements
      std::map<Node*, Vector<ELEMENT*>*> adjacent_elements_pt;
 
      // Allocate space for the vertex nodes
      Vector<Node*> vertex_node_pt;
 
      // Set up the patches information
      setup_patches(mesh_pt, adjacent_elements_pt, vertex_node_pt);
 
      // Grab any element (this shouldn't be a null pointer)
      ELEMENT* const el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(0));
 
      // Get the index of the vorticity
      unsigned smoothed_vorticity_index = el_pt->smoothed_vorticity_index();
 
      // Maximum order of vorticity derivative (that can be recovered)
      unsigned max_vort_order =
        el_pt->get_maximum_order_of_recoverable_vorticity_derivative();
 
      // Maximum order of velocity derivative (that can be recovered)
      unsigned max_veloc_order =
        el_pt->get_maximum_order_of_recoverable_velocity_derivative();
 
      // Maximum number of recoverable vorticity terms
      unsigned max_vort_recov = 0;
 
      // Maximum number of recoverable velocity terms
      unsigned max_veloc_recov = 0;
 
      // Loop over the entries of the vector
      for (unsigned i = 0; i < max_vort_order + 1; i++)
      {
        // Get the number of partial derivatives of the vorticity
        max_vort_recov += el_pt->npartial_derivative(i);
      }
 
      // Loop over the entries of the vector
      for (unsigned i = 1; i < max_veloc_order + 1; i++)
      {
        // Get the number of partial derivatives of the vorticity
        max_veloc_recov += 2 * el_pt->npartial_derivative(i);
      }
 
      // Number of recovered vorticity derivatives
      unsigned n_recovered_vort_derivs =
        el_pt->nvorticity_derivatives_to_recover();
 
      // Number of recovered velocity derivatives
      unsigned n_recovered_veloc_derivs =
        el_pt->nvelocity_derivatives_to_recover();
 
      // Determine number of coefficients for expansion of recovered vorticity.
      // Use complete polynomial of given order for recovery
      unsigned num_recovery_terms = nrecovery_order();
 
      // Counter for averaging of recovered vorticity and its derivatives
      std::map<Node*, unsigned> count;
 
      // Counter for which nodal value we're assigning
      unsigned nodal_dof = 0;
 
      // Loop over derivatives
      for (unsigned deriv = 0; deriv < max_vort_recov + max_veloc_recov;
           deriv++)
      {
        // If we're not recovering this vorticity derivative. Note, we cast
        // to an int because n_recovered_vort_derivs can be zero (so subtracting
        // any positive integer can cause trouble)
        if ((int(deriv) > int(n_recovered_vort_derivs - 1)) &&
            (deriv < max_vort_recov))
        {
          // We're done here
          continue;
        }
        // If we're not recovering any of the velocity derivatives and we're
        // finished with the vorticity derivatives
        else if ((n_recovered_veloc_derivs == 0) && (deriv >= max_vort_recov))
        {
          // We're done here
          continue;
        }
 
        // Storage for accumulated nodal vorticity (used to compute nodal
        // averages)
        std::map<Node*, double> averaged_recovered_vort;
 
        // Calculation of vorticity
        //-------------------------
        // Do patch recovery
        // unsigned counter=0;
        for (typename std::map<Node*, Vector<ELEMENT*>*>::iterator it =
               adjacent_elements_pt.begin();
             it != adjacent_elements_pt.end();
             it++)
        {
          // Pointer to the recovered vorticity coefficients
          Vector<double>* recovered_vorticity_coefficient_pt;
 
          // Setup smoothed vorticity field for patches
          get_recovered_vorticity_in_patch(*(it->second),
                                           num_recovery_terms,
                                           recovered_vorticity_coefficient_pt,
                                           deriv);
 
          // Now get the nodal average of the recovered vorticity (nodes are
          // generally part of multiple patches):
 
          // Get the number of elements in adjacent_elements_pt
          unsigned nelem = (*(it->second)).size();
 
          // Loop over all elements to get recovered vorticity
          for (unsigned e = 0; e < nelem; e++)
          {
            // Get pointer to element
            ELEMENT* const el_pt = (*(it->second))[e];
 
            // Get the number of nodes by element
            unsigned nnode_el = el_pt->nnode();
 
            // Loop over the nodes in the element
            for (unsigned j = 0; j < nnode_el; j++)
            {
              // Get a pointer to the j-th node in this element
              Node* nod_pt = el_pt->node_pt(j);
 
              // Get the local coordinates of the node
              el_pt->local_coordinate_of_node(j, s);
 
              // Interpolate the global (Eulerian) coordinate
              el_pt->interpolated_x(s, x);
 
              // Recovery shape functions at global (Eulerian) coordinate
              Vector<double> psi_r(num_recovery_terms);
 
              // Recover the shape function values at the position x
              shape_rec(x, psi_r);
 
              // Initialise the value of the recovered quantity
              double recovered_vort = 0.0;
 
              // Loop over the recovery terms
              for (unsigned i = 0; i < num_recovery_terms; i++)
              {
                // Assemble recovered vorticity
                recovered_vort +=
                  (*recovered_vorticity_coefficient_pt)[i] * psi_r[i];
              }
 
              // Keep adding
              averaged_recovered_vort[nod_pt] += recovered_vort;
 
              // Increment the counter
              count[nod_pt]++;
            } // for (unsigned j=0;j<nnode_el;j++)
          } // for (unsigned e=0;e<nelem;e++)
 
          // Delete the recovered coefficient data
          delete recovered_vorticity_coefficient_pt;
 
          // Make it a null pointer
          recovered_vorticity_coefficient_pt = 0;
        } // for (typename std::map<Node*,Vector<ELEMENT*>*>::iterator it=...
 
        // Find out how many nodes there are in the mesh
        unsigned nnod = mesh_pt->nnode();
 
        // Loop over all nodes to actually work out the average
        for (unsigned j = 0; j < nnod; j++)
        {
          // Make a pointer to the j-th node
          Node* nod_pt = mesh_pt->node_pt(j);
 
          // Calculate the values of the smoothed vorticity
          averaged_recovered_vort[nod_pt] /= count[nod_pt];
 
          // Assign smoothed vorticity to nodal values
          nod_pt->set_value(smoothed_vorticity_index + nodal_dof,
                            averaged_recovered_vort[nod_pt]);
        }
 
        // We're done with this dof so increment the counter
        nodal_dof++;
 
        // Start again
        count.clear();
      } // for (unsigned deriv=0;deriv<max_vort_recov+max_veloc_recov;deriv++)
 
      // Cleanup
      for (typename std::map<Node*, Vector<ELEMENT*>*>::iterator it =
             adjacent_elements_pt.begin();
           it != adjacent_elements_pt.end();
           it++)
      {
        // Delete the vector of element pointers
        delete it->second;
      }
 
      // Inform the user
      oomph_info << "Time for vorticity recovery [sec]: "
                 << TimingHelpers::timer() - t_start << std::endl;
    } // End of recover_vorticity

References e(), oomph::Mesh::element_pt(), oomph::VorticitySmoother< ELEMENT >::get_recovered_vorticity_in_patch(), i, j, oomph::Mesh::nnode(), oomph::Mesh::node_pt(), oomph::VorticitySmoother< ELEMENT >::nrecovery_order(), oomph::oomph_info, s, oomph::Data::set_value(), oomph::VorticitySmoother< ELEMENT >::setup_patches(), oomph::VorticitySmoother< ELEMENT >::shape_rec(), oomph::TimingHelpers::timer(), and plotDoE::x.

recovery_order()

template<class ELEMENT >

unsigned& oomph::VorticitySmoother< ELEMENT >::recovery_order ( )

inline

Access function for order of recovery polynomials.

    {
      // Return the order of recovery
      return Recovery_order;
    }

References oomph::VorticitySmoother< ELEMENT >::Recovery_order.

Referenced by oomph::VorticitySmoother< ELEMENT >::integral_rec(), and oomph::VorticitySmoother< ELEMENT >::shape_rec().

setup_patches()

template<class ELEMENT >

void oomph::VorticitySmoother< ELEMENT >::setup_patches ( Mesh *&  mesh_pt, std::map< Node *, Vector< ELEMENT * > * > &  adjacent_elements_pt, Vector< Node * > &  vertex_node_pt  )

inline

Setup patches: For each vertex node pointed to by nod_pt, adjacent_elements_pt[nod_pt] contains the pointer to the vector that contains the pointers to the elements that the node is part of. Also returns a Vector of vertex nodes for use in get_element_errors.

    {
      // Clear: hierher should we do this in Z2 as well?
      adjacent_elements_pt.clear();
 
      // Auxiliary map that contains element-adjacency for ALL nodes
      std::map<Node*, Vector<ELEMENT*>*> aux_adjacent_elements_pt;
 
#ifdef PARANOID
      // Check if all elements request the same recovery order
      unsigned ndisagree = 0;
#endif
 
      // Loop over all elements to setup adjacency for all nodes.
      // Need to do this because midside nodes can be corner nodes for
      // adjacent smaller elements! Admittedly, the inclusion of interior
      // nodes is wasteful...
      unsigned nelem = mesh_pt->nelement();
      for (unsigned e = 0; e < nelem; e++)
      {
        ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
 
#ifdef PARANOID
        // Check if all elements request the same recovery order
        if (el_pt->nrecovery_order() != Recovery_order)
        {
          ndisagree++;
        }
#endif
 
        // Loop all nodes in element
        unsigned nnod = el_pt->nnode();
        for (unsigned n = 0; n < nnod; n++)
        {
          // Make a pointer to the n-th node
          Node* nod_pt = el_pt->node_pt(n);
 
          // Has this node been considered before?
          if (aux_adjacent_elements_pt[nod_pt] == 0)
          {
            // Create Vector of pointers to its adjacent elements
            aux_adjacent_elements_pt[nod_pt] = new Vector<ELEMENT*>;
          }
 
          // Add pointer to adjacent element
          (*aux_adjacent_elements_pt[nod_pt]).push_back(el_pt);
        }
      } // end element loop
 
#ifdef PARANOID
      // Check if all elements request the same recovery order
      if (ndisagree != 0)
      {
        oomph_info
          << "\n\n======================================================\n"
          << "WARNING:\n"
          << ndisagree << " out of " << mesh_pt->nelement() << " elements"
          << "\nhave different preferences for the order of the recovery"
          << "\nshape functions. We are using: Recovery_order="
          << Recovery_order << std::endl;
        oomph_info
          << "======================================================\n\n";
      }
#endif
 
      // Loop over all elements, extract adjacency for corner nodes only
      nelem = mesh_pt->nelement();
      for (unsigned e = 0; e < nelem; e++)
      {
        ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
 
        // Loop over corner nodes
        unsigned n_node = el_pt->nvertex_node();
        for (unsigned n = 0; n < n_node; n++)
        {
          Node* nod_pt = el_pt->vertex_node_pt(n);
 
          // Has this node been considered before?
          if (adjacent_elements_pt[nod_pt] == 0)
          {
            // Add the node pointer to the vertex node container
            vertex_node_pt.push_back(nod_pt);
 
            // Create Vector of pointers to its adjacent elements
            adjacent_elements_pt[nod_pt] = new Vector<ELEMENT*>;
 
            // Copy across:
            unsigned nel = (*aux_adjacent_elements_pt[nod_pt]).size();
            for (unsigned e = 0; e < nel; e++)
            {
              (*adjacent_elements_pt[nod_pt])
                .push_back((*aux_adjacent_elements_pt[nod_pt])[e]);
            }
          }
        }
      } // End of loop over elements
 
      // Cleanup
      for (typename std::map<Node*, Vector<ELEMENT*>*>::iterator it =
             aux_adjacent_elements_pt.begin();
           it != aux_adjacent_elements_pt.end();
           it++)
      {
        delete it->second;
      }
    } // End of setup_patches

References e(), oomph::Mesh::element_pt(), n, oomph::Mesh::nelement(), oomph::oomph_info, oomph::VorticitySmoother< ELEMENT >::Recovery_order, and size.

Referenced by oomph::VorticitySmoother< ELEMENT >::recover_vorticity().

shape_rec()

template<class ELEMENT >

void oomph::VorticitySmoother< ELEMENT >::shape_rec ( const Vector< double > &  x, Vector< double > &  psi_r  )

inline

Recovery shape functions as functions of the global, Eulerian coordinate x of dimension dim. The recovery shape functions are complete polynomials of the order specified by Recovery_order.

    {
      // Create an ostringstream object to create a string
      std::ostringstream error_stream;
 
      // Find order of recovery shape functions
      switch (recovery_order())
      {
        case 1:
          // Complete linear polynomial in 2D:
          psi_r[0] = 1.0;
          psi_r[1] = x[0];
          psi_r[2] = x[1];
          break;
 
        case 2:
          // Complete quadratic polynomial in 2D:
          psi_r[0] = 1.0;
          psi_r[1] = x[0];
          psi_r[2] = x[1];
          psi_r[3] = x[0] * x[0];
          psi_r[4] = x[0] * x[1];
          psi_r[5] = x[1] * x[1];
          break;
 
        case 3:
          // Complete cubic polynomial in 2D:
          psi_r[0] = 1.0;
          psi_r[1] = x[0];
          psi_r[2] = x[1];
          psi_r[3] = x[0] * x[0];
          psi_r[4] = x[0] * x[1];
          psi_r[5] = x[1] * x[1];
          psi_r[6] = x[0] * x[0] * x[0];
          psi_r[7] = x[0] * x[0] * x[1];
          psi_r[8] = x[0] * x[1] * x[1];
          psi_r[9] = x[1] * x[1] * x[1];
          break;
 
        default:
          // Create an error message for this case
          error_stream << "Recovery shape functions for recovery order "
                       << recovery_order()
                       << " haven't yet been implemented for 2D" << std::endl;
 
          // Throw an error
          throw OomphLibError(error_stream.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
      }
    } // End of shape_rec

References OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::VorticitySmoother< ELEMENT >::recovery_order(), and plotDoE::x.

Referenced by oomph::VorticitySmoother< ELEMENT >::get_recovered_vorticity_in_patch(), and oomph::VorticitySmoother< ELEMENT >::recover_vorticity().

Member Data Documentation

Recovery_order

template<class ELEMENT >

unsigned oomph::VorticitySmoother< ELEMENT >::Recovery_order

private

Order of recovery polynomials.

Referenced by oomph::VorticitySmoother< ELEMENT >::nrecovery_order(), oomph::VorticitySmoother< ELEMENT >::recovery_order(), and oomph::VorticitySmoother< ELEMENT >::setup_patches().

---

The documentation for this class was generated from the following file:

• vorticity_smoother.h