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"

#include <vector>
#include <limits>

namespace apf {

class Field;
class Element;
class Mesh;
class MeshEntity;
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();
  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);


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))

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 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