TriNedelecBasis.cpp

#include "TriNedelecBasis.h"

TriNedelecBasis::TriNedelecBasis(void){
  // Set Basis Type //
  order   = 1;
  type    = 1;
  size    = 3;
  nodeNbr = 3;
  dim     = 2;

  // Lagrange //
  Polynomial* lagrange    = new Polynomial[3];

  lagrange[0] = 
    Polynomial(1, 0, 0, 0) - Polynomial(1, 1, 0, 0) - Polynomial(1, 0, 1, 0);

  lagrange[1] = 
    Polynomial(1, 1, 0, 0);

  lagrange[2] = 
    Polynomial(1, 0, 1, 0);

  // Basis //
  basis = new std::vector<std::vector<Polynomial> >(size);

  for(int i = 0; i < size; i++)
    (*basis)[i].resize(3);

  for(int i = 0, j = 1; i < 3; i++, j = (j + 1) % 3){
    std::vector<Polynomial> tmp = lagrange[j].gradient();
    tmp[0].mul(lagrange[i]);
    tmp[1].mul(lagrange[i]);
    tmp[2].mul(lagrange[i]);

    (*basis)[i] = lagrange[i].gradient();

    (*basis)[i][0].mul(lagrange[j]);
    (*basis)[i][1].mul(lagrange[j]);
    (*basis)[i][2].mul(lagrange[j]);
   
    (*basis)[i][0].sub(tmp[0]);
    (*basis)[i][1].sub(tmp[1]);
    (*basis)[i][2].sub(tmp[2]);
  }

  // Free Temporary Sapce //
  delete[] lagrange;
}

TriNedelecBasis::~TriNedelecBasis(void){
  delete basis;
}

/*
#include <cstdio>
int main(void){
  const int P = 1;
  const double d = 0.05;
  const char x[2] = {'X', 'Y'};

  TriNedelecBasis b;
  
  const std::vector<std::vector<Polynomial> >& basis = b.getBasis();
  
  printf("\n");
  printf("clear all;\n");
  printf("\n");
  
  printf("\n");
  printf("Order      = %d\n", b.getOrder());