/*
  Finite Element demos

  compile with:  'make fem'
 */

#include <iostream>
#include <sstream>
#include <iomanip>

#include <time.h>
#include <math.h>
#include <complex>



// ng-soft header files
#include <fem.hpp>


using namespace std;
using namespace ngstd;
using namespace ngbla;
using namespace ngfem;


namespace netgen {
  ostream * testout = &cout;
}


int main ()
{
  // ***************** test integration rule

  // integrate  x*y  over triangle [(0,0), (1,0), (0,1)]
  // with integrationrule of order 100 

  const IntegrationRule & ir = 
    GetIntegrationRules().SelectIntegrationRule (ET_TRIG, 100);
  cout << "number of ipts = " << ir.GetNIP() << endl;

  double sum = 0;
  for (int i = 0; i < ir.GetNIP(); i++)
    {
      const IntegrationPoint & ip = ir[i];
      sum += ip.Weight() * ip(0) * ip(1);
    }

  cout << "Integral = " << sum << " =?= " << 1.0/24.0 << endl << endl;


  // **************** Compute shape functions

  FE_Quad1 quad;
  Vector<> shapes(quad.GetNDof());
  for (double x = 0; x <= 1.0001; x += 0.25)
    for (double y = 0; y <= 1.0001; y += 0.25)
      {
	IntegrationPoint ip (x, y);
	quad.CalcShape (ip, shapes);
	cout << x << ", " << y << ": " << shapes << endl;
      }
  cout << endl << endl;


  // ***************** element matrix assembling

  // own memory management
  LocalHeap lh(10000);

  // reference finite element
  FE_Trig1 trig;

  // integrators for (\nabla u, \nabla v) and (u,v)
  ConstantCoefficientFunction coef(1);
  LaplaceIntegrator<2> laplace(&coef);
  MassIntegrator<2> mass(&coef);


  // element geometry:
  ElementTransformation eltrans;

  eltrans.SetElement (&trig, 0, 0);

  // vertex coordinates  
  double pts[3][2] = { { 0, 0 }, { 1, 0 }, { 0, 1 } };
  eltrans.AllocPointMatrix (2,3);
  eltrans.PointMatrix() = Trans (FlatMatrix<double> (3, 2, &pts[0][0]));

  cout << "PointMatrix = " << endl << eltrans.PointMatrix() << endl;


  FlatMatrix<double> elmat_lap;
  laplace.AssembleElementMatrix (trig, eltrans, elmat_lap, lh);

  cout << "elmat laplace = " << endl << elmat_lap << endl;



  FlatMatrix<double> elmat_mass;
  mass.AssembleElementMatrix (trig, eltrans, elmat_mass, lh);
  
  cout << "elmat mass = " << endl << elmat_mass << endl;



  return 0;
}