adaptiveData.cpp

// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.

#include <math.h>
#include <list>
#include <set>
#include "adaptiveData.h"
#include "Plugin.h"
#include "OS.h"
#include "GmshDefines.h"

//#define TIMER

std::set<adaptivePoint> adaptiveLine::allPoints;
std::set<adaptivePoint> adaptiveTriangle::allPoints;
std::set<adaptivePoint> adaptiveQuadrangle::allPoints;
std::set<adaptivePoint> adaptiveTetrahedron::allPoints;
std::set<adaptivePoint> adaptiveHexahedron::allPoints;
std::set<adaptivePoint> adaptivePrism::allPoints;

std::list<adaptiveLine*> adaptiveLine::all;
std::list<adaptiveTriangle*> adaptiveTriangle::all;
std::list<adaptiveQuadrangle*> adaptiveQuadrangle::all;
std::list<adaptiveTetrahedron*> adaptiveTetrahedron::all;
std::list<adaptiveHexahedron*> adaptiveHexahedron::all;
std::list<adaptivePrism*> adaptivePrism::all;

int adaptiveLine::numNodes = 2;
int adaptiveTriangle::numNodes = 3;
int adaptiveQuadrangle::numNodes = 4;
int adaptivePrism::numNodes = 6;
int adaptiveTetrahedron::numNodes = 4;
int adaptiveHexahedron::numNodes = 8;

int adaptiveLine::numEdges = 1;
int adaptiveTriangle::numEdges = 3;
int adaptiveQuadrangle::numEdges = 4;
int adaptivePrism::numEdges = 9;
int adaptiveTetrahedron::numEdges = 6;
int adaptiveHexahedron::numEdges = 12;

template <class T>
static void cleanElement()
{  
  for(typename std::list<T*>::iterator it = T::all.begin(); it != T::all.end(); ++it)
    delete *it;
  T::all.clear();
  T::allPoints.clear();
}

static void computeShapeFunctions(fullMatrix<double> *coeffs, fullMatrix<double> *eexps,
                                  double u, double v, double w, fullVector<double> *sf,
                                  fullVector<double> *tmp)
{
  for(int i = 0; i < eexps->size1(); i++) {
    (*tmp)(i) = pow(u, (*eexps)(i, 0));
    if(eexps->size2() > 1) (*tmp)(i) *= pow(v, (*eexps)(i, 1));
    if(eexps->size2() > 2) (*tmp)(i) *= pow(w, (*eexps)(i, 2));
  }
  coeffs->mult(*tmp, *sf);
}

adaptivePoint *adaptivePoint::add(double x, double y, double z,
                                  std::set<adaptivePoint> &allPoints)
{
  adaptivePoint p;
  p.x = x;
  p.y = y;