AdaptiveViews.cpp

//
// Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA.
// 
// Please report all bugs and problems to <gmsh@geuz.org>.

#include <stdio.h>
#include <math.h>
#include <list>
#include <set>
#include "AdaptiveViews.h"
#include "Plugin.h"

// A recursive effective implementation

void computeShapeFunctions(Double_Matrix * coeffs, Double_Matrix * eexps,
                           double u, double v, double w, double *sf);

std::set < adapt_point > adapt_point::all_points;
std::list < adapt_edge * >adapt_edge::all_elems;
std::list < adapt_triangle * >adapt_triangle::all_elems;
std::list < adapt_tet * >adapt_tet::all_elems;
std::list < adapt_quad * >adapt_quad::all_elems;
std::list < adapt_hex * >adapt_hex::all_elems;
#define MAX_NB_NOD 8
int adapt_edge::nbNod = 2;
int adapt_triangle::nbNod = 3;
int adapt_tet::nbNod = 4;
int adapt_quad::nbNod = 4;
int adapt_hex::nbNod = 8;

adapt_point *adapt_point::New(double x, double y, double z,
                              Double_Matrix * coeffs, Double_Matrix * eexps)
{
  adapt_point p;
  p.x = x;
  p.y = y;
  p.z = z;
  std::set < adapt_point >::iterator it = all_points.find(p);
  if(it == all_points.end()) {
    all_points.insert(p);
    it = all_points.find(p);
    double *kkk = (double *)(it->shape_functions);
    computeShapeFunctions(coeffs, eexps, x, y, z, kkk);
    return (adapt_point *) & (*it);
  }
  else
    return (adapt_point *) & (*it);
}

void adapt_edge::Create(int maxlevel, Double_Matrix * coeffs,
                        Double_Matrix * eexps)
{
  int level = 0;
  cleanElement < adapt_edge > ();
  adapt_point *p1 = adapt_point::New(-1, 0, 0, coeffs, eexps);