apf.h

Go to the documentation of this file.

/*
 * Copyright 2011 Scientific Computation Research Center
 *
 * This work is open source software, licensed under the terms of the
 * BSD license as described in the LICENSE file in the top-level directory.
 */
 
#ifndef APF_H
#define APF_H
 
#include "apfMatrix.h"
#include "apfNew.h"
#include "apfDynamicArray.h"
 
namespace pcu{
  class PCU;
}
 
#include <vector>
#include <map>
#include <limits>
 
namespace apf {
 
class Field;
class Element;
class Mesh;
class MeshEntity;
class MeshTag;
class VectorElement;
typedef VectorElement MeshElement;
class FieldShape;
struct Sharing;
template <class T> class ReductionOp;
template <class T> class ReductionSum;
 
template <class T>
class ReductionOp
{
  public:
    /* \brief apply operation, returning a single value */
    virtual T apply(T val1, T val2) const = 0;
 
    /* \brief returns a value such that apply(val, neutral_val) == val */
    virtual T getNeutralElement() const = 0;
};
 
template <class T>
class ReductionSum : public ReductionOp<T>
{
  T apply(T val1, T val2) const { return val1 + val2; };
  T getNeutralElement() const { return 0; };
};
 
template <class T>
class ReductionMin : public ReductionOp<T>
{
  T apply(T val1, T val2) const { return ( (val1 < val2) ? val1 : val2 ); };
  T getNeutralElement() const { return std::numeric_limits<T>::max(); };
};
 
template <class T>
class ReductionMax : public ReductionOp<T>
{
  T apply(T val1, T val2) const { return ( (val1 < val2) ? val2 : val1 ); };
  T getNeutralElement() const { return std::numeric_limits<T>::min(); };
};
 
void destroyMesh(Mesh* m);
 
MeshElement* createMeshElement(Mesh* m, MeshEntity* e);
 
MeshElement* createMeshElement(Field* c, MeshEntity* e);
 
MeshEntity * getMeshEntity(MeshElement * me);
 
void destroyMeshElement(MeshElement* e);
 
enum ValueType {
  SCALAR,
  VECTOR,
  MATRIX,
  PACKED,
  VALUE_TYPES
};
 
Field* createLagrangeField(Mesh* m, const char* name, int valueType, int order);
 
Field* createStepField(Mesh* m, const char* name, int valueType);
 
Field* createIPField(Mesh* m, const char* name, int valueType, int order);
 
Field* createField(Mesh* m, const char* name, int valueType, FieldShape* shape);
 
Field* createFieldOn(Mesh* m, const char* name, int valueType);
 
Field* createPackedField(Mesh* m, const char* name, int components,
    apf::FieldShape* shape = 0);
 
Field* createGeneralField(
    Mesh* m,
    const char* name,
    int valueType,
    int components,
    FieldShape* shape);
 
Field* cloneField(Field* f, Mesh* onto);
 
Mesh* getMesh(Field* f);
 
bool hasEntity(Field* f, MeshEntity* e);
 
const char* getName(Field* f);
 
int getValueType(Field* f);
 
void destroyField(Field* f);
 
void setScalar(Field* f, MeshEntity* e, int node, double value);
 
double getScalar(Field* f, MeshEntity* e, int node);
 
void setVector(Field* f, MeshEntity* e, int node, Vector3 const& value);
 
void getVector(Field* f, MeshEntity* e, int node, Vector3& value);
 
void setMatrix(Field* f, MeshEntity* e, int node, Matrix3x3 const& value);
 
void getMatrix(Field* f, MeshEntity* e, int node, Matrix3x3& value);
 
void setComponents(Field* f, MeshEntity* e, int node, double const* components);
 
void getComponents(Field* f, MeshEntity* e, int node, double* components);
 
Element* createElement(Field* f, MeshElement* e);
 
Element* createElement(Field* f, MeshEntity* e);
 
void destroyElement(Element* e);
 
MeshElement* getMeshElement(Element* e);
 
MeshEntity* getMeshEntity(Element* e);
 
double getScalar(Element* e, Vector3 const& param);
 
void getGrad(Element* e, Vector3 const& param, Vector3& grad);
 
void getVector(Element* e, Vector3 const& param, Vector3& value);
 
double getDiv(Element* e, Vector3 const& param);
 
void getCurl(Element* e, Vector3 const& param, Vector3& curl);
 
void getVectorGrad(Element* e, Vector3 const& param, Matrix3x3& deriv);
 
void getMatrix(Element* e, Vector3 const& param, Matrix3x3& value);
 
void getMatrixGrad(Element* e, Vector3 const& param, Vector<27>& value);
 
void getComponents(Element* e, Vector3 const& param, double* components);
 
int countIntPoints(MeshElement* e, int order);
 
void getIntPoint(MeshElement* e, int order, int point, Vector3& param);
 
double getIntWeight(MeshElement* e, int order, int point);
 
void mapLocalToGlobal(MeshElement* e, Vector3 const& local, Vector3& global);
 
double getDV(MeshElement* e, Vector3 const& param);
 
class Integrator
{
  public:
    Integrator(int o);
    virtual ~Integrator();
    void process(Mesh* m, int dim=-1);
    void process(MeshElement* e);
    virtual void inElement(MeshElement*);
    virtual void outElement();
    virtual void atPoint(Vector3 const& p, double w, double dV) = 0;
    virtual void parallelReduce(pcu::PCU*);
  protected:
    int order;
    int ipnode;
};
 
double measure(MeshElement* e);
 
double measure(Mesh* m, MeshEntity* e);
 
int getOrder(MeshElement* e);
 
void getJacobian(MeshElement* e, Vector3 const& local, Matrix3x3& j);
 
void getJacobianInv(MeshElement* e, Vector3 const& local, Matrix3x3& jinv);
 
double computeCosAngle(Mesh* m, MeshEntity* pe, MeshEntity* e1, MeshEntity* e2,
    const Matrix3x3& Q);
 
double computeCosAngleInTri(Mesh* m, MeshEntity* tri,
    MeshEntity* e1, MeshEntity* e2, const Matrix3x3& Q);
 
double computeCosAngleInTet(Mesh* m, MeshEntity* tet,
    MeshEntity* e1, MeshEntity* e2, const Matrix3x3& Q);
 
Vector3 computeFaceNormalAtEdgeInTet(Mesh* m, MeshEntity* tet,
    MeshEntity* face, MeshEntity* edge, Matrix3x3 Q);
 
Vector3 computeFaceNormalAtVertex(Mesh* m, MeshEntity* face,
    MeshEntity* vert, const Matrix3x3& Q);
 
Vector3 computeEdgeTangentAtVertex(Mesh* m, MeshEntity* edge,
    MeshEntity* vert, const Matrix3x3& Q);
 
double computeShortestHeightInTet(Mesh* m, MeshEntity* tet,
    const Matrix3x3& Q = Matrix3x3(1.,0.,0.,0.,1.,0.,0.,0.,1.));
 
double computeLargestHeightInTet(Mesh* m, MeshEntity* tet,
    const Matrix3x3& Q = Matrix3x3(1.,0.,0.,0.,1.,0.,0.,0.,1.));
 
 
double computeShortestHeightInTri(Mesh* m, MeshEntity* tri,
    const Matrix3x3& Q = Matrix3x3(1.,0.,0.,0.,1.,0.,0.,0.,1.));
 
double computeLargestHeightInTri(Mesh* m, MeshEntity* tri,
    const Matrix3x3& Q = Matrix3x3(1.,0.,0.,0.,1.,0.,0.,0.,1.));
 
 
int countNodes(Element* e);
 
void getScalarNodes(Element* e, NewArray<double>& values);
 
void getVectorNodes(Element* e, NewArray<Vector3>& values);
 
void getMatrixNodes(Element* e, NewArray<Matrix3x3>& values);
 
void getShapeValues(Element* e, Vector3 const& local,
    NewArray<double>& values);
 
void getShapeGrads(Element* e, Vector3 const& local,
    NewArray<Vector3>& grads);
 
void getVectorShapeValues(Element* e, Vector3 const& local,
    NewArray<Vector3>& values);
 
void getCurlShapeValues(Element* e, Vector3 const& local,
    NewArray<Vector3>& values);
 
FieldShape* getShape(Field* f);
 
int countComponents(Field* f);
 
#define APF_ITERATE(t,w,i) \
for (t::iterator i = (w).begin(); \
     (i) != (w).end(); ++(i))
 
#define APF_CONST_ITERATE(t,w,i) \
for (t::const_iterator i = (w).begin(); \
     (i) != (w).end(); ++(i))
 
struct CGNSInfo
{
  // cgns_bc_name
  std::string cgnsBCSName;
  /* tag value
 
    Tag value holds [0, 1] as a
    marker to indicate mesh_entities 
    within bc group. 1="in group", 0="not in group"
    Tags set on vertices, edges, faces, and cells
  */
  apf::MeshTag* bcsMarkerTag = nullptr;
  // model dimension
  int mdlDim = -1;
  // model id
  int mdlId = -1;
};
 
//using CGNSBCMap = std::map<std::string, std::vector<std::tuple<std::string, apf::MeshTag *, int>>>;
using CGNSBCMap = std::map<std::string, std::vector<CGNSInfo>>;
void writeCGNS(const char *prefix, Mesh *m, const CGNSBCMap &cgnsBCMap);
 
void writeVtkFiles(const char* prefix, Mesh* m, int cellDim = -1);
 
void writeVtkFiles(const char* prefix, Mesh* m,
    std::vector<std::string> writeFields, int cellDim = -1);
 
void writeOneVtkFile(const char* prefix, Mesh* m);
 
void writeASCIIVtkFiles(const char* prefix, Mesh* m);
 
void writeASCIIVtkFiles(const char* prefix, Mesh* m,
    std::vector<std::string> writeFields);
 
void writeNedelecVtkFiles(const char* prefix, Mesh* m);
 
void getGaussPoint(int type, int order, int point, Vector3& param);
 
int countGaussPoints(int type, int order);
 
double getJacobianDeterminant(Matrix3x3 const& J, int dimension);
 
int getDimension(MeshElement* me);
 
void synchronize(Field* f, Sharing* shr = 0);
 
void accumulate(Field* f, Sharing* shr = 0, bool delete_shr = false);
 
void sharedReduction(Field* f, Sharing* shr, bool delete_shr,
           const ReductionOp<double>& sum = ReductionSum<double>());
 
bool isPrintable(Field* f);
 
void fail(const char* why) __attribute__((noreturn));
 
void freeze(Field* f);
 
void unfreeze(Field* f);
 
bool isFrozen(Field* f);
 
double* getArrayData(Field* f);
 
void zeroField(Field* f);
 
struct Function
{
  virtual ~Function();
  virtual void eval(MeshEntity* e, double* result) = 0;
};
 
/* \brief Create a Field from a user's analytic function.
 * \details This field will use no memory, has values on all
 * nodes, and calls the user Function for all value queries.
 * Writing to this field does nothing.
 * \warning if you copy the mesh with apf::convert you may want to use
 * apf::updateUserField to update function that this field refers to. This is
 * extremely important if the analytic function you use references user fields.
 */
Field* createUserField(Mesh* m, const char* name, int valueType, FieldShape* s,
    Function* f);
 
/* \brief update the analytic function from a user field
 * \details this field updates the apf::Function which the UserField refers to
 */
void updateUserField(Field* field, Function* newFunc);
 
Field* recoverGradientByVolume(Field* f);
 
void copyData(Field* to, Field* from);
 
void projectField(Field* to, Field* from);
 
void axpy(double a, Field* x, Field* y);
 
void renameField(Field* f, const char* name);
 
void getBF(FieldShape* s, MeshElement* e, Vector3 const& p,
    NewArray<double>& BF);
 
void getGradBF(FieldShape* s, MeshElement* e, Vector3 const& p,
    NewArray<Vector3>& gradBF);
 
}
 
#endif