Line3D.h

Go to the documentation of this file.

 
#pragma once
 
#include "utils/helpers.h"
#include <jem/numeric/algebra/MatmulChain.h>
#include <jem/util/Properties.h>
#include <jive/Array.h>
#include <jive/geom/ParametricLine.h>
#include <jive/geom/Shape.h>
#include <jive/geom/StdShape.h>
 
using jem::idx_t;
using jem::String;
using jem::numeric::MatmulChain;
using jem::util::Properties;
using jive::Cubix;
using jive::Matrix;
using jive::Vector;
using jive::geom::ParametricLine;
using jive::geom::Shape;
using jive_helpers::Quadix;
 
class Line3D : public Shape
{
 public:
  static const String PROP_NODES; 
  static const String PROP_INT;   
  static const String PROP_LOCAL; 
  static const idx_t GLOB_RANK;   
 
  explicit Line3D(const String &name,
                  const Properties &conf,
                  const Properties &props);
 
  inline idx_t globalRank() const override;
 
  inline idx_t localRank() const override;
 
  inline idx_t nodeCount() const override;
 
  inline idx_t vertexCount() const override;
 
  inline idx_t ipointCount() const override;
 
  inline idx_t shapeFuncCount() const override;
 
  inline Matrix getVertexCoords() const override;
 
  inline Matrix getIntegrationScheme() const override;
 
  inline bool containsLocalPoint(const Vector &u) const override;
 
  void getIntegrationWeights(const Vector &w,
                             const Matrix &c) const override;
 
  inline Matrix getVertexFunctions() const override;
 
  inline Matrix getShapeFunctions() const override;
 
  void getGlobalPoint(const Vector &x,
                      const Vector &u,
                      const Matrix &c) const override;
 
  void getShapeGradients(const Matrix &g,
                         const Vector &w,
                         const Matrix &c) const;
 
  void getRotations(const Cubix &Ri,
                    const Cubix &Rn) const;
 
  void getXi(const Quadix &Xi,
             const Vector &w,
             const Matrix &u,
             const Matrix &c) const;
 
  void getXi(const Quadix &Xi,
             const Matrix &grads,
             const Matrix &phiP) const;
 
  void getPsi(const Quadix &Psi,
              const Vector &w,
              const Matrix &c) const;
 
  void getPi(const Cubix &Pi,
             const Cubix &Ri,
             const Cubix &Rn) const;
 
  void getPi(const Cubix &Pi,
             const Cubix &Rn) const;
 
  void getRotationGradients(const Cubix &LambdaP,
                            const Vector &w,
                            const Matrix &c,
                            const Cubix &nodeLambda) const;
 
  void getRotationsAndGradients(const Cubix &Ri,
                                const Cubix &LambdaP,
                                const Matrix &grads,
                                const Cubix &nodeLambda) const;
 
  //------------------------------------------------
  // local location functions
  //------------------------------------------------
  inline void evalShapeFunctions(const Vector &h,
                                 const Vector &u) const override;
 
  inline void evalShapeGradients(const Vector &h,
                                 const Vector &g,
                                 const Vector &u) const;
 
  inline void evalShapeGradGrads(const Vector &h,
                                 const Vector &g,
                                 const Vector &gg,
                                 const Vector &u) const;
 
 private:
  void getRotations_(const Cubix &Ri,
                     const Matrix &node_psi,
                     const Matrix &Lambda_r) const;
 
  void getRotationGradients_(const Cubix &LambdaP,
                             const Matrix &grads,
                             const Matrix &ip_psi,
                             const Matrix &node_psi,
                             const Matrix &Lambda_r) const;
 
  void getRefRot_(const Matrix Lambda_r,
                  const Cubix Rn) const;
 
  void getNodeRotVecs_(const Matrix psi,
                       const Matrix Lambda_r,
                       const Cubix Rn) const;
 
 protected:
  jem::Ref<ParametricLine> intLine_; 
 
  // Work arrays and matrix product of the functions that run for every element
  // (see the note in utils/helpers.h): created once in the constructor and
  // overwritten in every call. The object is not thread safe.
  // The result of matmul_ is only valid until its next call.
  mutable Matrix nodePsi_;                
  mutable Matrix lambdaR_;                
  mutable Matrix ipPsi_;                  
  mutable MatmulChain<double, 2> matmul_; 
  mutable Matrix relRot_;                 
  mutable Matrix skewP_;                  
  mutable Vector phiIJ_;                  
 
  mutable bool gradMemoValid_ = false;
  mutable Matrix gradMemoCoords_;
  mutable Matrix gradMemoGrads_;
  mutable Vector gradMemoWeights_;
};
 
//=====================================================================
//    inline implementations
//=====================================================================
idx_t Line3D::globalRank() const
{
  return GLOB_RANK;
}
 
idx_t Line3D::localRank() const
{
  return intLine_->localRank();
}
 
idx_t Line3D::nodeCount() const
{
  return intLine_->nodeCount();
}
 
idx_t Line3D::vertexCount() const
{
  return intLine_->vertexCount();
}
 
idx_t Line3D::ipointCount() const
{
  return intLine_->ipointCount();
}
 
idx_t Line3D::shapeFuncCount() const
{
  return intLine_->shapeFuncCount();
}
 
Matrix Line3D::getVertexCoords() const
{
  return intLine_->getVertexCoords();
}
 
Matrix Line3D::getIntegrationScheme() const
{
  return intLine_->getIntegrationScheme();
}
 
bool Line3D::containsLocalPoint(const Vector &u) const
{
  return intLine_->containsLocalPoint(u);
}
 
Matrix Line3D::getVertexFunctions() const
{
  return intLine_->getVertexFunctions();
}
 
void Line3D::evalShapeFunctions(const Vector &h,
                                const Vector &u) const
{
  JEM_ASSERT2(u.size() == localRank(), "internal coordinates do not match the rank");
  intLine_->evalShapeFunctions(h, u);
}
 
void Line3D::evalShapeGradients(const Vector &h,
                                const Vector &g,
                                const Vector &u) const
{
  JEM_ASSERT2(u.size() == localRank(), "internal coordinates do not match the rank");
  Matrix g_dummy(1, g.size());
  intLine_->getSShape()->evalShapeGradients(h, g_dummy, u);
  g = g_dummy(0, ALL);
}
 
void Line3D::evalShapeGradGrads(const Vector &h,
                                const Vector &g,
                                const Vector &gg,
                                const Vector &u) const
{
  JEM_ASSERT2(u.size() == localRank(), "internal coordinates do not match the rank");
  Matrix g_dummy(1, g.size());
  Matrix gg_dummy(1, gg.size());
  intLine_->getSShape()->evalShapeGradGrads(h, g_dummy, gg_dummy, u);
  g = g_dummy(0, ALL);
  gg = gg_dummy(0, ALL);
}
 
Matrix Line3D::getShapeFunctions() const
{
  return intLine_->getShapeFunctions();
}