src/fem/nedelecelement.cpp

Go to the documentation of this file.

/***************************************************************************
                          nedelecelement.cpp  -  description
                             -------------------
    begin                : Fri Dec 19 2003
    copyright            : (C) 2003 by Roman Geus
    email                : roman.geus@psi.ch
***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/

#include <iomanip>
#include "colarray/linalg.h"
#include "nedelecelement.h"

#define CROSS(A,B,k) (A((k+1)%3,0)*B((k+2)%3,0) - A((k+2)%3,0)*B((k+1)%3,0))
#define CROSS2(A,B,c,k) (A((k+1)%3,c)*B((k+2)%3,0) - A((k+2)%3,c)*B((k+1)%3,0))

using namespace colarray;
using namespace mesh;

#ifndef GAUSSIANQUADRATURE
NedelecElement::NedelecElement()
{ }

#else
NedelecElement::NedelecElement() :
  _mu(3,3),
  _eps(3,3),
  _Cur(3,24),
  _Ned(3,24)
{
  _intOrder = 0;

  // CAUTION: We use a fixed order _order. In case of higher order
  // basis functions, this has to be adapted.
  constructGQ(4);
}
#endif


NedelecElement::~NedelecElement()
{ }




//
// The following code is used in case of numerical quadrature ...
//

void NedelecElement::constructGQ(int order)
{
  _intOrder = order;
  Matrix<double> J(order, order);
  Matrix<double> Q(order, order);
  ColumnVector<double> xi1D(order);
  ColumnVector<double> w1D(order);

  // Start with determining the 1D quadrature rule for weight 1 (Legendre)
  J.fill(0.0);
  J(0,0) = 0.5;
  for(int i=1; i<order; i++){
    J(i,i)   = 0.5;
    J(i-1,i) = (double)i/sqrt(16.0*i*i - 4.0);
    J(i,i-1) = (double)i/sqrt(16.0*i*i - 4.0);
  }
  linalg::eigenvalueDecomp(J, Q, xi1D);
  for(int i=0; i<order; i++){
    w1D[i] = Q(0,i)*Q(0,i);
  }


  // Now, with the 1D rule, construct the tensor product rule,
  // i.e. the evaluation points and the weights for a unit tetrahedron
  _w.resize(order*order*order);
  _xi.resize(order*order*order);
  _eta.resize(order*order*order);
  _zeta.resize(order*order*order);
  int l = 0;
  for(int i=0; i<order; i++){
    for(int j=0; j<order; j++){
      for(int k=0; k<order; k++){
        _w[l]    = w1D[i]*w1D[j]*w1D[k]*(1.0 - xi1D[i])*(1.0 - xi1D[i])*(1.0 - xi1D[j]);
        _xi[l]   = xi1D[i];
        _eta[l]  = xi1D[j]*(1.0 - xi1D[i]);
        _zeta[l] = xi1D[k]*(1.0 - xi1D[j])*(1.0 - xi1D[i]);
        l++;
      }
    }    
  }

}


void NedelecElement::evalElementFunctions(Matrix<double>& xi, Matrix<double>& iM, 
                                           Matrix<double>& Cur, Matrix<double>& Ned) 
{

  // Compute the nodal funcions Li and their gradients Gi
  double  L[4];
  Matrix<double> G1(3,1), G2(3,1), G3(3,1), G4(3,1);
  L[0] = 1.0 - xi(0,0) - xi(1,0) - xi(2,0);
  L[1] = xi(0,0);
  L[2] = xi(1,0);
  L[3] = xi(2,0);
  for(int i=0; i<3; i++){
    G1(i,0) = -iM(0,i) - iM(1,i) - iM(2,i);
    G2(i,0) =  iM(0,i);
    G3(i,0) =  iM(1,i);
    G4(i,0) =  iM(2,i);
  }


  // Compute the linear Nedelec functions [1..6]
  for(int i=0; i<3; i++){
    Ned(i,0) = L[0]*G2(i,0) - L[1]*G1(i,0);
    Ned(i,1) = L[0]*G3(i,0) - L[2]*G1(i,0);
    Ned(i,2) = L[0]*G4(i,0) - L[3]*G1(i,0);
    Ned(i,3) = L[1]*G3(i,0) - L[2]*G2(i,0);
    Ned(i,4) = L[1]*G4(i,0) - L[3]*G2(i,0);
    Ned(i,5) = L[2]*G4(i,0) - L[3]*G3(i,0);
  }


  //  Now compute the curls of the Nedelec functions, i.e. curl(N(:,i)
  // 
  //  Note, that curl(f x F) = grad(f) x F  + f * curl(F) and since F
  //  is a gradient field in all the expressions above, i.e. 
  //  curl(grad(phi)) = 0, we can throw away the better part of the 
  //  expressions. We are left with curl(f x F) = grad(f) x F. Fortunately, 
  //  most of the gradients have already been computed ...
  //  
  // 
  for(int i=0; i<3; i++){
    Cur(i,0)  = CROSS(G1,G2,i) - CROSS(G2,G1,i); 
    Cur(i,1)  = CROSS(G1,G3,i) - CROSS(G3,G1,i); 
    Cur(i,2)  = CROSS(G1,G4,i) - CROSS(G4,G1,i);
    Cur(i,3)  = CROSS(G2,G3,i) - CROSS(G3,G2,i);
    Cur(i,4)  = CROSS(G2,G4,i) - CROSS(G4,G2,i);
    Cur(i,5)  = CROSS(G3,G4,i) - CROSS(G4,G3,i);
  }


  //  Compute the quadratic, edge based Nedelec functions [7..12]
  for(int i=0; i<3; i++){
    Ned(i,6)  = L[0]*G2(i,0) + L[1]*G1(i,0);
    Ned(i,7)  = L[0]*G3(i,0) + L[2]*G1(i,0);
    Ned(i,8)  = L[0]*G4(i,0) + L[3]*G1(i,0);
    Ned(i,9)  = L[1]*G3(i,0) + L[2]*G2(i,0);
    Ned(i,10) = L[1]*G4(i,0) + L[3]*G2(i,0);
    Ned(i,11) = L[2]*G4(i,0) + L[3]*G3(i,0);
  }


  // Compute the quadratic, face based and interior Nedelec functions [13..24]
  for(int i=0; i<3; i++){
    Ned(i,12) = L[1]*L[2]*G1(i,0);
    Ned(i,13) = L[2]*L[3]*G2(i,0);
    Ned(i,14) = L[3]*L[0]*G3(i,0);
    Ned(i,15) = L[0]*L[1]*G4(i,0);
    
    Ned(i,16) = L[2]*L[0]*G2(i,0);
    Ned(i,17) = L[3]*L[1]*G3(i,0);
    Ned(i,18) = L[0]*L[2]*G4(i,0);
    Ned(i,19) = L[1]*L[3]*G1(i,0);
    
    Ned(i,20) = L[0]*L[1]*G3(i,0);
    Ned(i,21) = L[1]*L[2]*G4(i,0);
    Ned(i,22) = L[2]*L[3]*G1(i,0);
    Ned(i,23) = L[3]*L[0]*G2(i,0);    
  }


  // Compute the curls of the second order edge based Nedelec element functions
  for(int i=0; i<3; i++){
    Cur(i,6)   = CROSS(G1,G2,i) + CROSS(G2,G1,i); 
    Cur(i,7)   = CROSS(G1,G3,i) + CROSS(G3,G1,i); 
    Cur(i,8)   = CROSS(G1,G4,i) + CROSS(G4,G1,i);
    Cur(i,9)   = CROSS(G2,G3,i) + CROSS(G3,G2,i);
    Cur(i,10)  = CROSS(G2,G4,i) + CROSS(G4,G2,i);
    Cur(i,11)  = CROSS(G3,G4,i) + CROSS(G4,G3,i);
  }  


  // Compute the curls of the second order face and interior Nedelec element functions
  // Note, that grad(Li*Lj) = grad(Li)*Lj + Li*grad(Lj), a quantity which
  // has been computed in the expressions N(:,7) ... N(:,12)
  //
  // pair = [
  //     0  7  8  9
  //     7  0 10 11
  //     8 10  0 12
  //     9 11 12  0
  //        ];
  for(int i=0; i<3; i++){
    Cur(i,12) = CROSS2(Ned,G1, 9,i);
    Cur(i,13) = CROSS2(Ned,G2,11,i);
    Cur(i,14) = CROSS2(Ned,G3, 8,i);
    Cur(i,15) = CROSS2(Ned,G4, 6,i);
    
    Cur(i,16) = CROSS2(Ned,G2, 7,i);
    Cur(i,17) = CROSS2(Ned,G3,10,i);
    Cur(i,18) = CROSS2(Ned,G4, 7,i);
    Cur(i,19) = CROSS2(Ned,G1,10,i);

    Cur(i,20) = CROSS2(Ned,G3, 6,i);
    Cur(i,21) = CROSS2(Ned,G4, 9,i);
    Cur(i,22) = CROSS2(Ned,G1,11,i);
    Cur(i,23) = CROSS2(Ned,G2, 8,i);
  }  
}


void NedelecElement::get_AeMe(Tet* tet,
                               int order,
                               Matrix<double>& A, Matrix<double>& M) {

  assert(_intOrder > 0);
  
  // Retrieve the reference to the mesh this tet is part of
  TetMesh* tetMesh = tet->get_mesh();

  // Size depending on the element order
  int size;
  if (order==1) {
    size = 6;
  }else{
    size = 24;
  }

  // Construct the mapping matrix
  Matrix<double> Map(3,3), iMap(3,3);
  double         detMap;
  Vector3        v0, vi;
  Matrix<double> xi(3,1);

  v0 =  tet->get_corner(0)->get_coord();
  for(int i=1; i<4; i++){
    vi =  tet->get_corner(i)->get_coord();
    Map(0,i-1) = vi.x - v0.x;
    Map(1,i-1) = vi.y - v0.y;
    Map(2,i-1) = vi.z - v0.z;
  }
  detMap = linalg::determinant(Map);
  linalg::inverse(Map, iMap);

  // Initialise the matrix to be returned with zeros
  A.fill(0.0);
  M.fill(0.0);

  // Iterate over all points to be evaluated and compute first the
  // material tensor (mu) and the the element function values
  for(int i=0; i<_intOrder*_intOrder*_intOrder; i++){

    // Construct the evaluation point ...
    xi(0,0) = _xi[i];
    xi(1,0) = _eta[i];
    xi(2,0) = _zeta[i];

    // Compute the Nedelec functions and their curls
    evalElementFunctions(xi, iMap, _Cur, _Ned); 

    // Compute the material tensors _mu and _eps
    tetMesh->get_materials()->getParameters(tet->get_material_ID(), xi, _mu, _eps);

    // Now build the bilinear form (evaluated at position xi[i]
    for(int ii=0; ii<size; ii++){
      for(int jj=0; jj<size; jj++){
        A(ii,jj) += detMap*_w[i]*(
                                  _Cur(0,ii)*(1.0/_mu(0,0))*_Cur(0,jj) + 
                                  _Cur(1,ii)*(1.0/_mu(1,1))*_Cur(1,jj) +
                                  _Cur(2,ii)*(1.0/_mu(2,2))*_Cur(2,jj)
                                  );

        M(ii,jj) += detMap*_w[i]*(
                                  _Ned(0,ii)*_eps(0,0)*_Ned(0,jj) + 
                                  _Ned(1,ii)*_eps(1,1)*_Ned(1,jj) +
                                  _Ned(2,ii)*_eps(2,2)*_Ned(2,jj)
                                  );
      }      
    }
  }


  // Fix the edge orientation if necessary
  for (int k = 0; k < 6; k ++) {
    if (!tet->get_edge_orientation(k)) {
      for (int j = 0; j < A.get_m(); j ++){
        A(k,j) = -A(k,j);
        M(k,j) = -M(k,j);
      }
      for (int i = 0; i < A.get_m(); i ++){
        A(i,k) = -A(i,k);
        M(i,k) = -M(i,k);
      }
    }
  }
}






void NedelecElement::get_Ae(Tet* tet,
                            int order,
                            const Tet::GradientMatrix& J,
                            Matrix<double>& A) {
    // cross products of all combinations of gradients
    static const int ind1[] = {0, 0, 0, 1, 1, 2};
    static const int ind2[] = {1, 2, 3, 2, 3, 3};
    Matrix<double> c(3, 6);
    for (int ii = 0; ii < 6; ii ++) {
        int i = ind1[ii];
        int j = ind2[ii];
        c(0,ii) = J(1,i)*J(2,j) - J(2,i)*J(1,j);
        c(1,ii) = J(2,i)*J(0,j) - J(0,i)*J(2,j);
        c(2,ii) = J(0,i)*J(1,j) - J(1,i)*J(0,j);
    }

    // inner products of all cross products
    Matrix<double> C(6, 6);
    blas::gemm('t', 'n', c, c, C, 1.0, 0.0);

    double volume = tet->get_volume();
  
    // element matrix
    for (int i = 0; i < 6; i ++)
        for (int j = 0; j < 6; j ++)
            A(i,j) = C(i,j) * (4.0*volume);

    if (order == 2) {

        for (int i = 6; i < 12; i ++)
            for (int j = 0; j < 24; j ++)
                A(i,j) = 0.0;

        for (int i = 0; i < 24; i ++)
            for (int j = 6; j < 12; j ++)
                A(i,j) = 0.0;

        for (int i = 0; i < 6; i ++) {
            A(i,12) = (-C(i,0) - C(i,1))*(volume/2);
            A(i,13) = (-C(i,3) - C(i,4))*(volume/2);
            A(i,14) = (-C(i,5) + C(i,1))*(volume/2);
            A(i,15) = ( C(i,2) + C(i,4))*(volume/2);
            A(i,16) = (-C(i,3) + C(i,0))*(volume/2);
            A(i,17) = (-C(i,5) + C(i,3))*(volume/2);
            A(i,18) = ( C(i,2) + C(i,5))*(volume/2);
            A(i,19) = (-C(i,0) - C(i,2))*(volume/2);
            A(i,20) = ( C(i,1) + C(i,3))*(volume/2);
            A(i,21) = ( C(i,4) + C(i,5))*(volume/2);
            A(i,22) = (-C(i,1) - C(i,2))*(volume/2);
            A(i,23) = (-C(i,4) + C(i,0))*(volume/2);
        }

        for (int i = 0; i < 6; i ++)
            for (int j = 12; j < 24; j ++)
                A(j,i) = A(i,j);

        A(12,12) = 2*C(0,0) + 2*C(0,1) + 2*C(1,1);
        A(12,13) = C(0,3) + 2*C(0,4) + C(1,3) + C(1,4);
        A(12,14) = C(0,5) - C(0,1) + C(1,5) - C(1,1);
        A(12,15) = - C(0,2) - C(0,4) - 2*C(1,2) - C(1,4);
        A(12,16) = C(0,3) - 2*C(0,0) + C(1,3) - C(1,0);
        A(12,17) = C(0,5) - C(0,3) + 2*C(1,5) - C(1,3);
        A(12,18) = -2*C(0,2) - C(0,5) - C(1,2) - C(1,5);
        A(12,19) = C(0,0) + C(0,2) + C(1,0) + 2*C(1,2);
        A(12,20) = -C(0,1) - C(0,3) - 2*C(1,1) - C(1,3);
        A(12,21) = -2*C(0,4) - C(0,5) - C(1,4) - 2*C(1,5);
        A(12,22) = C(0,1) + 2*C(0,2) + C(1,1) + C(1,2);
        A(12,23) = C(0,4) - C(0,0) + C(1,4) - C(1,0);
        ;
        A(13,13) = 2*C(3,3) + 2*C(3,4) + 2*C(4,4);
        A(13,14) = C(3,5) - 2*C(3,1) + C(4,5) - C(4,1);
        A(13,15) = - C(3,2) - C(3,4) - C(4,2) - C(4,4);
        A(13,16) = C(3,3) - C(3,0) + C(4,3) - 2*C(4,0);
        A(13,17) = C(3,5) - 2*C(3,3) + C(4,5) - C(4,3);
        A(13,18) = -C(3,2) - C(3,5) - 2*C(4,2) - C(4,5);
        A(13,19) = 2*C(3,0) + C(3,2) + C(4,0) + C(4,2);
        A(13,20) = -C(3,1) - C(3,3) - C(4,1) - C(4,3);
        A(13,21) = -C(3,4) - C(3,5) - 2*C(4,4) - C(4,5);
        A(13,22) = 2*C(3,1) + C(3,2) + C(4,1) + 2*C(4,2);
        A(13,23) = C(3,4) - 2*C(3,0) + C(4,4) - C(4,0);
        ;
        A(14,14) = 2*C(5,5) - 2*C(5,1) + 2*C(1,1);
        A(14,15) = - C(5,2) - 2*C(5,4) + C(1,2) + C(1,4);
        A(14,16) = 2*C(5,3) - C(5,0) - C(1,3) + C(1,0);
        A(14,17) = C(5,5) - C(5,3) - C(1,5) + 2*C(1,3);
        A(14,18) = -C(5,2) - 2*C(5,5) + C(1,2) + C(1,5);
        A(14,19) = C(5,0) + C(5,2) - 2*C(1,0) - C(1,2);
        A(14,20) = -C(5,1) - 2*C(5,3) + C(1,1) + C(1,3);
        A(14,21) = -C(5,4) - C(5,5) + C(1,4) + C(1,5);
        A(14,22) = C(5,1) + C(5,2) - 2*C(1,1) - C(1,2);
        A(14,23) = 2*C(5,4) - C(5,0) - C(1,4) + 2*C(1,0);
        ;
        A(15,15) = 2*C(2,2) + 2*C(2,4) + 2*C(4,4);
        A(15,16) = - C(2,3) + C(2,0) - 2*C(4,3) + C(4,0);
        A(15,17) = - 2*C(2,5) + C(2,3) - C(4,5) + C(4,3);
        A(15,18) = C(2,2) + C(2,5) + C(4,2) + 2*C(4,5);
        A(15,19) = - C(2,0) - 2*C(2,2) - C(4,0) - C(4,2);
        A(15,20) = 2*C(2,1) + C(2,3) + C(4,1) + 2*C(4,3);
        A(15,21) = C(2,4) + 2*C(2,5) + C(4,4) + C(4,5);
        A(15,22) = -C(2,1) - C(2,2) - C(4,1) - C(4,2);
        A(15,23) = -C(2,4) + C(2,0) - 2*C(4,4) + C(4,0);
        ;
        A(16,16) = 2*C(3,3) - 2*C(0,3) + 2*C(0,0);
        A(16,17) = C(3,5) - C(3,3) - C(0,5) + C(0,3);
        A(16,18) = - C(3,2) - 2*C(3,5) + 2*C(0,2) + C(0,5);
        A(16,19) = C(3,0) + C(3,2) - C(0,0) - C(0,2);
        A(16,20) = -C(3,1) - 2*C(3,3) + C(0,1) + C(0,3);
        A(16,21) = -C(3,4) - C(3,5) + 2*C(0,4) + C(0,5);
        A(16,22) = C(3,1) + C(3,2) - C(0,1) - 2*C(0,2);
        A(16,23) = 2*C(3,4) - C(3,0) - C(0,4) + C(0,0);
        ;
        A(17,17) = 2*C(5,5) - 2*C(5,3)  + 2*C(3,3);
        A(17,18) = - C(5,2) - C(5,5) + C(3,2) + C(3,5);
        A(17,19) = C(5,0) + 2*C(5,2) - 2*C(3,0) - C(3,2);
        A(17,20) = -2*C(5,1) - C(5,3) + C(3,1) + C(3,3);
        A(17,21) = -C(5,4) - 2*C(5,5) + C(3,4) + C(3,5);
        A(17,22) = C(5,1) + C(5,2) - 2*C(3,1) - C(3,2);
        A(17,23) = C(5,4) - C(5,0) - C(3,4) + 2*C(3,0);
        ;
        A(18,18) = 2*C(2,2) + 2*C(2,5) + 2*C(5,5);
        A(18,19) = - C(2,0) - C(2,2) - C(5,0) - C(5,2);
        A(18,20) = C(2,1) + C(2,3) + C(5,1) + 2*C(5,3);
        A(18,21) = 2*C(2,4) + C(2,5) + C(5,4) + C(5,5);
        A(18,22) = -C(2,1) - 2*C(2,2) - C(5,1) - C(5,2);
        A(18,23) = -C(2,4) + C(2,0) - 2*C(5,4) + C(5,0);
        ;
        A(19,19) = 2*C(0,0) + 2*C(0,2) + 2*C(2,2);
        A(19,20) = -C(0,1) - C(0,3) - 2*C(2,1) - C(2,3);
        A(19,21) = -C(0,4) - C(0,5) - C(2,4) - 2*C(2,5);
        A(19,22) = 2*C(0,1) + C(0,2) + C(2,1) + C(2,2);
        A(19,23) = C(0,4) - 2*C(0,0) + C(2,4) - C(2,0);
        ;
        A(20,20) = 2*C(1,1) + 2*C(1,3) + 2*C(3,3);
        A(20,21) = C(1,4) + 2*C(1,5) + C(3,4) + C(3,5);
        A(20,22) = -C(1,1) - C(1,2) - C(3,1) - C(3,2);
        A(20,23) = -C(1,4) + C(1,0) - 2*C(3,4) + C(3,0);
        ;
        A(21,21) = 2*C(4,4) + 2*C(4,5) + 2*C(5,5);
        A(21,22) = -C(4,1) - 2*C(4,2) - C(5,1) - C(5,2);
        A(21,23) = -C(4,4) + C(4,0) - C(5,4) + C(5,0);
        ;
        A(22,22) = 2*C(1,1) + 2*C(1,2) + 2*C(2,2);
        A(22,23) = C(1,4) - 2*C(1,0) + C(2,4) - C(2,0);
        ;
        A(23,23) = 2*C(4,4) - 2*C(4,0) + 2*C(0,0);

        for (int i = 13; i < 24; i ++)
            for (int j = 12; j < i; j ++)
                A(i,j) = A(j,i);

        for (int i = 12; i < 24; i ++)
            for (int j = 12; j < 24; j ++)
                A(i,j) *= volume/20.0;
    }

    // edge orientation
    for (int k = 0; k < 6; k ++) {
        if (!tet->get_edge_orientation(k)) {
            for (int j = 0; j < A.get_m(); j ++)
                A(k,j) = -A(k,j);
            for (int i = 0; i < A.get_m(); i ++)
                A(i,k) = -A(i,k);
        }
    }

}

void NedelecElement::get_Me(Tet* tet,
                            int order,
                            const Tet::GradientMatrix& J,
                            Matrix<double>& M) {
    // innner products of all gradients
    Matrix<double> G(4, 4);
    blas::gemm('t', 'n', J, J, G, 1.0, 0.0);

    // element matrix
    M(0,0) = (G(0,0)+G(1,1)-G(0,1))*84;
    M(0,1) = (G(0,0)-G(0,1)-G(0,2)+2*G(1,2))*42;
    M(0,2) = (G(0,0)-G(0,1)-G(0,3)+2*G(1,3))*42;
    M(0,3) = (-G(1,1)+G(0,1)+G(1,2)-2*G(0,2))*42;
    M(0,4) = (-G(1,1)+G(0,1)+G(1,3)-2*G(0,3))*42;
    M(0,5) = (G(0,2)+G(1,3)-G(0,3)-G(1,2))*42;
    ;
    M(1,1) = (G(0,0)+G(2,2)-G(0,2))*84;
    M(1,2) = (G(0,0)-G(0,2)-G(0,3)+2*G(2,3))*42;
    M(1,3) = (G(2,2)-G(0,2)-G(1,2)+2*G(0,1))*42;
    M(1,4) = (G(0,1)+G(2,3)-G(0,3)-G(1,2))*42;
    M(1,5) = (-G(2,2)+G(0,2)+G(2,3)-2*G(0,3))*42;
    ;
    M(2,2) = (G(0,0)+G(3,3)-G(0,3))*84;
    M(2,3) = (G(0,1)+G(2,3)-G(0,2)-G(1,3))*42;
    M(2,4) = (G(3,3)-G(0,3)-G(1,3)+2*G(0,1))*42;
    M(2,5) = (G(3,3)-G(0,3)-G(2,3)+2*G(0,2))*42;
    ;
    M(3,3) = (G(1,1)+G(2,2)-G(1,2))*84;
    M(3,4) = (G(1,1)-G(1,2)-G(1,3)+2*G(2,3))*42;
    M(3,5) = (-G(2,2)+G(1,2)+G(2,3)-2*G(1,3))*42;
    ;
    M(4,4) = (G(1,1)+G(3,3)-G(1,3))*84;
    M(4,5) = (G(3,3)-G(1,3)-G(2,3)+2*G(1,2))*42;
    ;
    M(5,5) = (G(2,2)+G(3,3)-G(2,3))*84;
    ;

    if (order == 2) {
        M(0,6) = (G(1,1)-G(0,0))*84;
        M(0,7) = (-G(0,0)+G(0,1)-G(0,2)+2*G(1,2))*42;
        M(0,8) = (-G(0,0)+G(0,1)-G(0,3)+2*G(1,3))*42;
        M(0,9) = (G(1,1)-G(0,1)+G(1,2)-2*G(0,2))*42;
        M(0,10) = (G(1,1)-G(0,1)+G(1,3)-2*G(0,3))*42;
        M(0,11) = (-G(0,2)-G(0,3)+G(1,2)+G(1,3))*42;
        ;
        M(0,12) = (G(0,1)-2*G(0,0))*7;
        M(0,13) = (G(1,1)-G(0,1))*7;
        M(0,14) = (-G(0,2)+2*G(1,2))*7;
        M(0,15) = (-G(0,3)+G(1,3))*14;
        M(0,16) = (2*G(1,1)-G(0,1))*7;
        M(0,17) = (-2*G(0,2)+G(1,2))*7;
        M(0,18) = (-G(0,3)+2*G(1,3))*7;
        M(0,19) = (-2*G(0,0)+G(0,1))*7;
        ;
        M(0,20) = (2*G(1,2)-2*G(0,2))*7;
        M(0,21) = (G(1,3)-2*G(0,3))*7;
        M(0,22) = (G(0,1)-G(0,0))*7;
        M(0,23) = (2*G(1,1)-G(0,1))*7;
        ;
        M(1,6) = (-G(0,0)-G(0,1)+G(0,2)+2*G(1,2))*42;
        M(1,7) = (G(2,2)-G(0,0))*84;
        M(1,8) = (-G(0,0)+G(0,2)-G(0,3)+2*G(2,3))*42;
        M(1,9) = (G(2,2)-G(0,2)+G(1,2)-2*G(0,1))*42;
        M(1,10) = (-G(0,1)-G(0,3)+G(1,2)+G(2,3))*42;
        M(1,11) = (G(2,2)-G(0,2)+G(2,3)-2*G(0,3))*42;
        ;
        M(1,12) = (-2*G(0,0)+G(0,2))*7;
        M(1,13) = (-2*G(0,1)+G(1,2))*7;
        M(1,14) = (2*G(2,2)-G(0,2))*7;
        M(1,15) = (-G(0,3)+2*G(2,3))*7;
        M(1,16) = (-G(0,1)+G(1,2))*14;
        M(1,17) = (G(2,2)-G(0,2))*7;
        M(1,18) = (-G(0,3)+G(2,3))*14;
        M(1,19) = (-G(0,0)+G(0,2))*7;
        ;
        M(1,20) = (2*G(2,2)-G(0,2))*7;
        M(1,21) = (G(2,3)-2*G(0,3))*7;
        M(1,22) = (G(2,0)-2*G(0,0))*7;
        M(1,23) = (2*G(1,2)-G(0,1))*7;
        ;
        M(2,6) = (-G(0,0)-G(0,1)+G(0,3)+2*G(1,3))*42;
        M(2,7) = (-G(0,0)-G(0,2)+G(0,3)+2*G(2,3))*42;
        M(2,8) = (G(3,3)-G(0,0))*84;
        M(2,9) = (-G(0,2)-G(0,1)+G(2,3)+G(1,3))*42;
        M(2,10) = (G(3,3)-G(0,3)+G(1,3)-2*G(0,1))*42;
        M(2,11) = (G(3,3)-G(0,3)+G(2,3)-2*G(0,2))*42;
        ;
        M(2,12) = (-G(0,0)+G(0,3))*7;
        M(2,13) = (-2*G(0,1)+G(1,3))*7;
        M(2,14) = (-G(0,2)+G(2,3))*14;
        M(2,15) = (2*G(3,3)-G(0,3))*7;
        M(2,16) = (-G(0,1)+2*G(1,3))*7;
        M(2,17) = (-2*G(0,2)+G(2,3))*7;
        M(2,18) = (2*G(3,3)-G(0,3))*7;
        M(2,19) = (-2*G(0,0)+G(0,3))*7;
        ;
        M(2,20) = (2*G(2,3)-G(0,2))*7;
        M(2,21) = (G(3,3)-G(0,3))*7;
        M(2,22) = (G(3,0)-2*G(0,0))*7;
        M(2,23) = (2*G(1,3)-2*G(0,1))*7;
        ;
        M(3,6) = (-G(1,1)-G(0,1)+G(1,2)+2*G(0,2))*42;
        M(3,7) = (G(2,2)+G(0,2)-G(1,2)-2*G(0,1))*42;
        M(3,8) = (-G(0,1)+G(0,2)-G(1,3)+G(2,3))*42;
        M(3,9) = (G(2,2)-G(1,1))*84;
        M(3,10) = (-G(1,1)+G(1,2)-G(1,3)+2*G(2,3))*42;
        M(3,11) = (G(2,2)-G(1,2)+G(2,3)-2*G(1,3))*42;
        ;
        M(3,12) = (-G(0,1)+G(0,2))*14;
        M(3,13) = (-2*G(1,1)+G(1,2))*7;
        M(3,14) = (G(2,2)-G(1,2))*7;
        M(3,15) = (-G(1,3)+2*G(2,3))*7;
        M(3,16) = (-2*G(1,1)+G(1,2))*7;
        M(3,17) = (2*G(2,2)-G(1,2))*7;
        M(3,18) = (-2*G(1,3)+G(2,3))*7;
        M(3,19) = (-G(0,1)+2*G(0,2))*7;
        ;
        M(3,20) = (2*G(2,2)-G(1,2))*7;
        M(3,21) = (2*G(2,3)-2*G(1,3))*7;
        M(3,22) = (G(2,0)-2*G(0,1))*7;
        M(3,23) = (G(1,2)-G(1,1))*7;
        ;
        M(4,6) = (-G(1,1)+G(1,3)-G(0,1)+2*G(0,3))*42;
        M(4,7) = (-G(0,1)+G(0,3)-G(1,2)+G(2,3))*42;
        M(4,8) = (G(3,3)+G(0,3)-G(1,3)-2*G(0,1))*42;
        M(4,9) = (-G(1,1)-G(1,2)+G(1,3)+2*G(2,3))*42;
        M(4,10) = (G(3,3)-G(1,1))*84;
        M(4,11) = (G(3,3)+G(2,3)-G(1,3)-2*G(1,2))*42;
        ;
        M(4,12) = (-G(0,1)+2*G(0,3))*7;
        M(4,13) = (-2*G(1,1)+G(1,3))*7;
        M(4,14) = (-2*G(1,2)+G(2,3))*7;
        M(4,15) = (2*G(3,3)-G(1,3))*7;
        M(4,16) = (-G(1,1)+G(1,3))*7;
        M(4,17) = (-G(1,2)+G(2,3))*14;
        M(4,18) = (G(3,3)-G(1,3))*7;
        M(4,19) = (-G(0,1)+G(0,3))*14;
        ;
        M(4,20) = (2*G(3,2)-G(1,2))*7;
        M(4,21) = (2*G(3,3)-G(1,3))*7;
        M(4,22) = (G(3,0)-2*G(0,1))*7;
        M(4,23) = (G(3,1)-2*G(1,1))*7;
        ;
        M(5,6) = (-G(0,2)+G(0,3)-G(1,2)+G(1,3))*42;
        M(5,7) = (-G(2,2)+G(2,3)-G(0,2)+2*G(0,3))*42;
        M(5,8) = (G(3,3)+G(0,3)-G(2,3)-2*G(0,2))*42;
        M(5,9) = (-G(2,2)+G(2,3)-G(1,2)+2*G(1,3))*42;
        M(5,10) = (G(3,3)+G(1,3)-G(2,3)-2*G(1,2))*42;
        M(5,11) = (G(3,3)-G(2,2))*84;
        ;
        M(5,12) = (-G(0,2)+2*G(0,3))*7;
        M(5,13) = (-G(1,2)+G(1,3))*14;
        M(5,14) = (-2*G(2,2)+G(2,3))*7;
        M(5,15) = (G(3,3)-G(2,3))*7;
        M(5,16) = (2*G(1,3)-G(1,2))*7;
        M(5,17) = (-2*G(2,2)+G(2,3))*7;
        M(5,18) = (2*G(3,3)-G(2,3))*7;
        M(5,19) = (-2*G(0,2)+G(0,3))*7;
        ;
        M(5,20) = (G(3,2)-G(2,2))*7;
        M(5,21) = (2*G(3,3)-G(2,3))*7;
        M(5,22) = (2*G(3,0)-2*G(0,2))*7;
        M(5,23) = (G(3,1)-2*G(1,2))*7;
        ;
        M(6,6) = (G(0,0)+G(1,1)+G(0,1))*84;
        M(6,7) = (G(0,0)+G(0,1)+G(0,2)+2*G(1,2))*42;
        M(6,8) = (G(0,0)+G(0,1)+G(0,3)+2*G(1,3))*42;
        M(6,9) = (G(1,1)+G(0,1)+G(1,2)+2*G(0,2))*42;
        M(6,10)= (G(1,1)+G(0,1)+G(1,3)+2*G(0,3))*42;
        M(6,11) =(G(0,2)+G(0,3)+G(1,2)+G(1,3))*42;
        ;
        M(6,12) = (G(0,1)+2*G(0,0))*7;
        M(6,13) = (G(1,1)+G(0,1))*7;
        M(6,14) = (G(0,2)+2*G(1,2))*7;
        M(6,15) = (G(0,3)+G(1,3))*14;
        M(6,16) = (2*G(1,1)+G(0,1))*7;
        M(6,17) = (2*G(0,2)+G(1,2))*7;
        M(6,18) = (G(0,3)+2*G(1,3))*7;
        M(6,19) = (2*G(0,0)+G(0,1))*7;
        ;
        M(6,20) = (2*G(1,2)+2*G(0,2))*7;
        M(6,21) = (G(1,3)+2*G(0,3))*7;
        M(6,22) = (G(0,1)+G(0,0))*7;
        M(6,23) = (2*G(1,1)+G(0,1))*7;
        ;
        M(7,7) = (G(0,0)+G(2,2)+G(0,2))*84;
        M(7,8) = (G(0,0)+G(0,2)+G(0,3)+2*G(2,3))*42;
        M(7,9) = (G(2,2)+G(0,2)+G(1,2)+2*G(0,1))*42;
        M(7,10) = (G(0,1)+G(0,3)+G(1,2)+G(2,3))*42;
        M(7,11) = (G(2,2)+G(0,2)+G(2,3)+2*G(0,3))*42;
        ;
        M(7,12) = (2*G(0,0)+G(0,2))*7;
        M(7,13) = (2*G(0,1)+G(1,2))*7;
        M(7,14) = (2*G(2,2)+G(0,2))*7;
        M(7,15) = (G(0,3)+2*G(2,3))*7;
        M(7,16) = (G(0,1)+G(1,2))*14;
        M(7,17) = (G(2,2)+G(0,2))*7;
        M(7,18) = (G(0,3)+G(2,3))*14;
        M(7,19) = (G(0,0)+G(0,2))*7;
        ;
        M(7,20) = (2*G(2,2)+G(0,2))*7;
        M(7,21) = (G(2,3)+2*G(0,3))*7;
        M(7,22) = (G(2,0)+2*G(0,0))*7;
        M(7,23) = (2*G(1,2)+G(0,1))*7;
        ;
        M(8,8) = (G(0,0)+G(3,3)+G(0,3))*84;
        M(8,9) = (G(0,1)+G(0,2)+G(1,3)+G(2,3))*42;
        M(8,10) = (G(3,3)+G(0,3)+G(1,3)+2*G(0,1))*42;
        M(8,11) = (G(3,3)+G(0,3)+G(2,3)+2*G(0,2))*42;
        ;
        M(8,12) = (G(0,0)+G(0,3))*7;
        M(8,13) = (2*G(0,1)+G(1,3))*7;
        M(8,14) = (G(0,2)+G(2,3))*14;
        M(8,15) = (2*G(3,3)+G(0,3))*7;
        M(8,16) = (G(0,1)+2*G(1,3))*7;
        M(8,17) = (2*G(0,2)+G(2,3))*7;
        M(8,18) = (2*G(3,3)+G(0,3))*7;
        M(8,19) = (2*G(0,0)+G(0,3))*7;
        ;
        M(8,20) = (2*G(2,3)+G(0,2))*7;
        M(8,21) = (G(3,3)+G(0,3))*7;
        M(8,22) = (G(3,0)+2*G(0,0))*7;
        M(8,23) = (2*G(1,3)+2*G(0,1))*7;
        ;
        M(9,9) = (G(1,1)+G(2,2)+G(1,2))*84;
        M(9,10) = (G(1,1)+G(1,2)+G(1,3)+2*G(2,3))*42;
        M(9,11) = (G(2,2)+G(1,2)+G(2,3)+2*G(1,3))*42;
        ;
        M(9,12) = (G(0,1)+G(0,2))*14;
        M(9,13) = (2*G(1,1)+G(1,2))*7;
        M(9,14) = (G(2,2)+G(1,2))*7;
        M(9,15) = (G(1,3)+2*G(2,3))*7;
        M(9,16) = (2*G(1,1)+G(1,2))*7;
        M(9,17) = (2*G(2,2)+G(1,2))*7;
        M(9,18) = (2*G(1,3)+G(2,3))*7;
        M(9,19) = (G(0,1)+2*G(0,2))*7;
        ;
        M(9,20) = (2*G(2,2)+G(1,2))*7;
        M(9,21) = (2*G(2,3)+2*G(1,3))*7;
        M(9,22) = (G(2,0)+2*G(0,1))*7;
        M(9,23) = (G(1,2)+G(1,1))*7;
        ;
        M(10,10) = (G(1,1)+G(3,3)+G(1,3))*84;
        M(10,11) = (G(3,3)+G(1,3)+G(2,3)+2*G(1,2))*42;
        ;
        M(10,12) = (G(0,1)+2*G(0,3))*7;
        M(10,13) = (2*G(1,1)+G(1,3))*7;
        M(10,14) = (2*G(1,2)+G(2,3))*7;
        M(10,15) = (2*G(3,3)+G(1,3))*7;
        M(10,16) = (G(1,1)+G(1,3))*7;
        M(10,17) = (G(1,2)+G(2,3))*14;
        M(10,18) = (G(3,3)+G(1,3))*7;
        M(10,19) = (G(0,1)+G(0,3))*14;
        ;
        M(10,20) = (2*G(3,2)+G(1,2))*7;
        M(10,21) = (2*G(3,3)+G(1,3))*7;
        M(10,22) = (G(3,0)+2*G(0,1))*7;
        M(10,23) = (G(3,1)+2*G(1,1))*7;
        ;
        M(11,11) = (G(2,2)+G(3,3)+G(2,3))*84;
        ;
        M(11,12) = (G(0,2)+2*G(0,3))*7;
        M(11,13) = (G(1,2)+G(1,3))*14;
        M(11,14) = (2*G(2,2)+G(2,3))*7;
        M(11,15) = (G(3,3)+G(2,3))*7;
        M(11,16) = (2*G(1,3)+G(1,2))*7;
        M(11,17) = (2*G(2,2)+G(2,3))*7;
        M(11,18) = (2*G(3,3)+G(2,3))*7;
        M(11,19) = (2*G(0,2)+G(0,3))*7;
        ;
        M(11,20) = (G(3,2)+G(2,2))*7;
        M(11,21) = (2*G(3,3)+G(2,3))*7;
        M(11,22) = (2*G(3,0)+2*G(0,2))*7;
        M(11,23) = (G(3,1)+2*G(1,2))*7;
        ;
        M(12,12) = G(0,0)*4;
        M(12,13) = G(0,1)*2;
        M(12,14) = G(0,2);
        M(12,15) = G(0,3)*2;
        M(12,16) = G(0,1)*2;
        M(12,17) = G(0,2)*2;
        M(12,18) = G(0,3)*2;
        M(12,19) = G(0,0)*2;
        M(12,20) = G(0,2)*2;
        M(12,21) = G(0,3)*4;
        M(12,22) = G(0,0)*2;
        M(12,23) = G(0,1);
        ;
        ;
        M(13,13) = G(1,1)*4;
        M(13,14) = G(1,2)*2;
        M(13,15) = G(1,3);
        M(13,16) = G(1,1)*2;
        M(13,17) = G(1,2)*2;
        M(13,18) = G(1,3)*2;
        M(13,19) = G(0,1)*2;
        M(13,20) = G(1,2);
        M(13,21) = G(1,3)*2;
        M(13,22) = G(1,0)*4;
        M(13,23) = G(1,1)*2;
        ;
        M(14,14) = G(2,2)*4;
        M(14,15) = G(2,3)*2;
        M(14,16) = G(1,2)*2;
        M(14,17) = G(2,2)*2;
        M(14,18) = G(2,3)*2;
        M(14,19) = G(0,2)*2;
        M(14,20) = G(2,2)*2;
        M(14,21) = G(2,3);
        M(14,22) = G(2,0)*2;
        M(14,23) = G(2,1)*4;
        ;
        M(15,15) = G(3,3)*4;
        M(15,16) = G(1,3)*2;
        M(15,17) = G(2,3)*2;
        M(15,18) = G(3,3)*2;
        M(15,19) = G(0,3)*2;
        M(15,20) = G(3,2)*4;
        M(15,21) = G(3,3)*2;
        M(15,22) = G(3,0);
        M(15,23) = G(3,1)*2;
        ;
        M(16,16) = G(1,1)*4;
        M(16,17) = G(1,2);
        M(16,18) = G(1,3)*4;
        M(16,19) = G(0,1);
        M(16,20) = G(1,2)*2;
        M(16,21) = G(1,3)*2;
        M(16,22) = G(1,0)*2;
        M(16,23) = G(1,1)*2;
        ;
        M(17,17) = G(2,2)*4;
        M(17,18) = G(2,3);
        M(17,19) = G(0,2)*4;
        M(17,20) = G(2,2)*2;
        M(17,21) = G(2,3)*2;
        M(17,22) = G(2,0)*2;
        M(17,23) = G(2,1)*2;
        ;
        M(18,18) = G(3,3)*4;
        M(18,19) = G(0,3);
        M(18,20) = G(3,2)*2;
        M(18,21) = G(3,3)*2;
        M(18,22) = G(3,0)*2;
        M(18,23) = G(3,1)*2;
        ;
        M(19,19) = G(0,0)*4;
        M(19,20) = G(0,2)*2;
        M(19,21) = G(0,3)*2;
        M(19,22) = G(0,0)*2;
        M(19,23) = G(0,1)*2;
        ;
        M(20,20) = G(2,2)*4;
        M(20,21) = G(2,3)*2;
        M(20,22) = G(2,0);
        M(20,23) = G(2,1)*2;
        ;
        M(21,21) = G(3,3)*4;
        M(21,22) = G(3,0)*2;
        M(21,23) = G(3,1);
        ;
        M(22,22) = G(0,0)*4;
        M(22,23) = G(0,1)*2;
        ;
        M(23,23) = G(1,1)*4;
    } // order == 2
  
    double volume = tet->get_volume();

    for (int i = 0; i < M.get_m(); i ++)
        for (int j = i; j < M.get_n(); j ++)
            M(i,j) *= volume/840.0;

    for (int i = 0; i < M.get_m(); i ++)
        for (int j = 0; j < i; j ++)
            M(i,j) = M(j,i);

    // edge orientation
    for (int k = 0; k < 6; k ++) {
        if (!tet->get_edge_orientation(k)) {
            for (int j = 0; j < M.get_m(); j ++)
                M(k,j) = -M(k,j);
            for (int i = 0; i < M.get_m(); i ++)
                M(i,k) = -M(i,k);
        }
    }
}


Vector3 NedelecElement::eval_cartesian(mesh::Tet* tet,
                                       int order,
                                       const Vector<double>& dof,
                                       const mesh::Tet::GradientMatrix& J,
                                       const Vector3& x) {
    Vector4 xi;
    tet->cartesian_to_simplex(x, xi);
    return eval_simplex(tet, order, dof, J, xi);
}
Vector3 NedelecElement::eval_curl_cartesian(mesh::Tet* tet,
                                            int order,
                                            const Vector<double>& dof,
                                            const mesh::Tet::GradientMatrix& J,
                                            const Vector3& x) {
    Vector4 xi;
    tet->cartesian_to_simplex(x, xi);
    return eval_curl_simplex(tet, order, dof, J, xi);
}

Vector3 NedelecElement::eval_simplex(mesh::Tet* tet,
                                     int order,
                                     const Vector<double>& dof,
                                     const mesh::Tet::GradientMatrix& J,
                                     const Vector4& xi) {
    Vector3 f(Vector3::ZERO);
    double coeff;

#define DOF_LIN_1(DOF, IND0, IND1) \
  coeff = (tet->get_edge_orientation(DOF) ? 1.0 : -1.0) * dof(DOF); \
  f.x += coeff * (xi[IND0] * J(0, IND1) - xi[IND1] * J(0, IND0)); \
  f.y += coeff * (xi[IND0] * J(1, IND1) - xi[IND1] * J(1, IND0)); \
  f.z += coeff * (xi[IND0] * J(2, IND1) - xi[IND1] * J(2, IND0));
  
#define DOF_LIN_2(DOF, IND0, IND1) \
  coeff = dof(DOF); \
  f.x += coeff * (xi[IND0] * J(0, IND1) + xi[IND1] * J(0, IND0)); \
  f.y += coeff * (xi[IND0] * J(1, IND1) + xi[IND1] * J(1, IND0)); \
  f.z += coeff * (xi[IND0] * J(2, IND1) + xi[IND1] * J(2, IND0));
  
#define DOF_QUAD(DOF, IND0, IND1, IND2) \
  coeff = dof(DOF); \
  f.x += coeff * J(0, IND0) * xi[IND1] * xi[IND2]; \
  f.y += coeff * J(1, IND0) * xi[IND1] * xi[IND2]; \
  f.z += coeff * J(2, IND0) * xi[IND1] * xi[IND2];

    DOF_LIN_1(0, 0, 1);
    DOF_LIN_1(1, 0, 2);
    DOF_LIN_1(2, 0, 3);
    DOF_LIN_1(3, 1, 2);
    DOF_LIN_1(4, 1, 3);
    DOF_LIN_1(5, 2, 3);

    if (order > 1) {

        DOF_LIN_2(6, 0, 1);
        DOF_LIN_2(7, 0, 2);
        DOF_LIN_2(8, 0, 3);
        DOF_LIN_2(9, 1, 2);
        DOF_LIN_2(10, 1, 3);
        DOF_LIN_2(11, 2, 3);

        DOF_QUAD(12, 0, 1, 2);
        DOF_QUAD(13, 1, 2, 3);
        DOF_QUAD(14, 2, 3, 0);
        DOF_QUAD(15, 3, 0, 1);

        DOF_QUAD(16, 1, 2, 0);
        DOF_QUAD(17, 2, 3, 1);
        DOF_QUAD(18, 3, 0, 2);
        DOF_QUAD(19, 0, 1, 3);

        DOF_QUAD(20, 2, 0, 1);
        DOF_QUAD(21, 3, 1, 2);
        DOF_QUAD(22, 0, 2, 3);
        DOF_QUAD(23, 1, 3, 0);
    
    }

#undef DOF_LIN_1
#undef DOF_LIN_2
#undef DOF_QUAD
  
    return f;
}

Vector3 NedelecElement::eval_curl_simplex(mesh::Tet* tet,
                                          int order,
                                          const Vector<double>& dof,
                                          const mesh::Tet::GradientMatrix& J,
                                          const Vector4& xi) {
    // Compute crossproducts of gradients of basis function
    //   C[i][j][:] = crossprod(grad(L_i), grad(L_j))
    Matrix<Vector3> C(4, 4);
    for (int i = 0; i < 4; i ++)
        for (int j = 0; j < 4; j ++) {
            C(i,j).x = J(1,i)*J(2,j) - J(2,i)*J(1,j);
            C(i,j).y = J(2,i)*J(0,j) - J(0,i)*J(2,j);
            C(i,j).z = J(0,i)*J(1,j) - J(1,i)*J(0,j);
        }

    Vector3 f(Vector3::ZERO);
    double coeff;

#define CURL_DOF_LIN_1(DOF, IND0, IND1) \
  coeff = (tet->get_edge_orientation(DOF) ? 1.0 : -1.0) * dof[DOF]; \
  f.x += coeff * 2 * C(IND0,IND1).x; \
  f.y += coeff * 2 * C(IND0,IND1).y; \
  f.z += coeff * 2 * C(IND0,IND1).z;
#define CURL_DOF_QUAD(DOF, IND0, IND1, IND2) \
  coeff = dof[DOF]; \
  f.x += coeff * (xi[IND1] * C(IND2,IND0).x + xi[IND2] * C(IND1,IND0).x); \
  f.y += coeff * (xi[IND1] * C(IND2,IND0).y + xi[IND2] * C(IND1,IND0).y); \
  f.z += coeff * (xi[IND1] * C(IND2,IND0).z + xi[IND2] * C(IND1,IND0).z);

    CURL_DOF_LIN_1(0, 0, 1);
    CURL_DOF_LIN_1(1, 0, 2);
    CURL_DOF_LIN_1(2, 0, 3);
    CURL_DOF_LIN_1(3, 1, 2);
    CURL_DOF_LIN_1(4, 1, 3);
    CURL_DOF_LIN_1(5, 2, 3);

    if (order > 1) {

        CURL_DOF_QUAD(12, 0, 1, 2);
        CURL_DOF_QUAD(13, 1, 2, 3);
        CURL_DOF_QUAD(14, 2, 3, 0);
        CURL_DOF_QUAD(15, 3, 0, 1);

        CURL_DOF_QUAD(16, 1, 2, 0);
        CURL_DOF_QUAD(17, 2, 3, 1);
        CURL_DOF_QUAD(18, 3, 0, 2);
        CURL_DOF_QUAD(19, 0, 1, 3);

        CURL_DOF_QUAD(20, 2, 0, 1);
        CURL_DOF_QUAD(21, 3, 1, 2);
        CURL_DOF_QUAD(22, 0, 2, 3);
        CURL_DOF_QUAD(23, 1, 3, 0);

    }
  
#undef CURL_DOF_LIN_1
#undef CURL_DOF_QUAD

    return f;
}

double NedelecElement::surface_integral_curl(mesh::Tet* tet,
                                             mesh::Face* face,
                                             const colarray::Vector<double>& dofs) {
    // allocate element matrix
    colarray::Matrix<double> A(get_nof_dofs(), get_nof_dofs());
    A = 0.0;

    // gradients of simplex coordinates
    mesh::Tet::GradientMatrix J(*tet);

    // cross products of all combinations of gradients
    static const int ind1[] = {0, 0, 0, 1, 1, 2};
    static const int ind2[] = {1, 2, 3, 2, 3, 3};
    colarray::Matrix<double> c(3, 6);
    for (int ii = 0; ii < 6; ii ++) {
        int i = ind1[ii];
        int j = ind2[ii];
        c(0,ii) = J(1,i)*J(2,j) - J(2,i)*J(1,j);
        c(1,ii) = J(2,i)*J(0,j) - J(0,i)*J(2,j);
        c(2,ii) = J(0,i)*J(1,j) - J(1,i)*J(0,j);
    }

    // inner products of all cross products
    colarray::Matrix<double> C(6, 6);
    blas::gemm('t', 'n', c, c, C, 1.0, 0.0);

    // compute area integration element
    mesh::Vector3 p0(face->get_corner(0)->get_coord());
    mesh::Vector3 p1(face->get_corner(1)->get_coord());
    mesh::Vector3 p2(face->get_corner(2)->get_coord());
    mesh::Vector3 u = p1 - p0;
    mesh::Vector3 v = p2 - p0;
    double E = u.dot_product(u);
    double F = u.dot_product(v);
    double G = v.dot_product(v);
    double S = sqrt(E*G - F*F);

    // get face index
    int face_id_local;
    for (face_id_local = 0; face_id_local < 4; ++ face_id_local)
        if (tet->get_face_id(face_id_local) == face->get_id())
            break;
    assert(face_id_local < 4);

    // first order part
    colarray::Matrix<double> A1(A, 0, 6, 0, 6);
    A1.inject(C);
    A1 *= 2.0*S;

    // second order part
    if (get_order() == 2) {
        switch (face_id_local) {
        case 0:
            A(0,12) = -S*(C(0,0)+C(0,1))/3.0;
            A(0,13) = -S*C(0,4)/3.0;
            A(0,14) = -S*C(0,5)/3.0;
            A(0,15) = S*(C(0,2)+C(0,4))/3.0;
            A(0,16) = S*(-C(0,3)+C(0,0))/3.0;
            A(0,17) = -S*C(0,5)/3.0;
            A(0,18) = S*(C(0,2)+C(0,5))/3.0;
            A(0,19) = -S*C(0,2)/3.0;
            A(0,20) = S*(C(0,1)+C(0,3))/3.0;
            A(0,21) = S*(C(0,4)+C(0,5))/3.0;
            A(0,22) = -S*C(0,2)/3.0;
            A(0,23) = -S*C(0,4)/3.0;

            A(1,12) = -S*(C(1,0)+C(1,1))/3.0;
            A(1,13) = -S*C(1,4)/3.0;
            A(1,14) = -S*C(1,5)/3.0;
            A(1,15) = S*(C(1,2)+C(1,4))/3.0;
            A(1,16) = S*(-C(1,3)+C(1,0))/3.0;
            A(1,17) = -S*C(1,5)/3.0;
            A(1,18) = S*(C(1,2)+C(1,5))/3.0;
            A(1,19) = -S*C(1,2)/3.0;
            A(1,20) = S*(C(1,1)+C(1,3))/3.0;
            A(1,21) = S*(C(1,4)+C(1,5))/3.0;
            A(1,22) = -S*C(1,2)/3.0;
            A(1,23) = -S*C(1,4)/3.0;

            A(2,12) = -S*(C(2,0)+C(2,1))/3.0;
            A(2,13) = -S*C(2,4)/3.0;
            A(2,14) = -S*C(2,5)/3.0;
            A(2,15) = S*(C(2,2)+C(2,4))/3.0;
            A(2,16) = S*(-C(2,3)+C(2,0))/3.0;
            A(2,17) = -S*C(2,5)/3.0;
            A(2,18) = S*(C(2,2)+C(2,5))/3.0;
            A(2,19) = -S*C(2,2)/3.0;
            A(2,20) = S*(C(2,1)+C(2,3))/3.0;
            A(2,21) = S*(C(2,4)+C(2,5))/3.0;
            A(2,22) = -S*C(2,2)/3.0;
            A(2,23) = -S*C(2,4)/3.0;

            A(3,12) = -S*(C(3,0)+C(3,1))/3.0;
            A(3,13) = -S*C(3,4)/3.0;
            A(3,14) = -S*C(3,5)/3.0;
            A(3,15) = S*(C(3,2)+C(3,4))/3.0;
            A(3,16) = S*(-C(3,3)+C(3,0))/3.0;
            A(3,17) = -S*C(3,5)/3.0;
            A(3,18) = S*(C(3,2)+C(3,5))/3.0;
            A(3,19) = -S*C(3,2)/3.0;
            A(3,20) = S*(C(3,1)+C(3,3))/3.0;
            A(3,21) = S*(C(3,4)+C(3,5))/3.0;
            A(3,22) = -S*C(3,2)/3.0;
            A(3,23) = -S*C(3,4)/3.0;

            A(4,12) = -S*(C(4,0)+C(4,1))/3.0;
            A(4,13) = -S*C(4,4)/3.0;
            A(4,14) = -S*C(4,5)/3.0;
            A(4,15) = S*(C(4,2)+C(4,4))/3.0;
            A(4,16) = S*(-C(4,3)+C(4,0))/3.0;
            A(4,17) = -S*C(4,5)/3.0;
            A(4,18) = S*(C(4,2)+C(4,5))/3.0;
            A(4,19) = -S*C(4,2)/3.0;
            A(4,20) = S*(C(4,1)+C(4,3))/3.0;
            A(4,21) = S*(C(4,4)+C(4,5))/3.0;
            A(4,22) = -S*C(4,2)/3.0;
            A(4,23) = -S*C(4,4)/3.0;

            A(5,12) = -S*(C(5,0)+C(5,1))/3.0;
            A(5,13) = -S*C(5,4)/3.0;
            A(5,14) = -S*C(5,5)/3.0;
            A(5,15) = S*(C(5,2)+C(5,4))/3.0;
            A(5,16) = S*(-C(5,3)+C(5,0))/3.0;
            A(5,17) = -S*C(5,5)/3.0;
            A(5,18) = S*(C(5,2)+C(5,5))/3.0;
            A(5,19) = -S*C(5,2)/3.0;
            A(5,20) = S*(C(5,1)+C(5,3))/3.0;
            A(5,21) = S*(C(5,4)+C(5,5))/3.0;
            A(5,22) = -S*C(5,2)/3.0;
            A(5,23) = -S*C(5,4)/3.0;

            A(12,12) = S*(C(0,0)+C(0,1)+C(1,1))/12.0;

            A(13,12) = S*(2.0*C(0,4)+C(1,4))/24.0;
            A(13,13) = S*C(4,4)/12.0;

            A(14,12) = S*(C(0,5)+C(1,5))/24.0;
            A(14,13) = S*C(4,5)/24.0;
            A(14,14) = S*C(5,5)/12.0;

            A(15,12) = -S*(C(0,2)+2.0*C(1,2)+C(0,4)+C(1,4))/24.0;
            A(15,13) = -S*(C(2,4)+C(4,4))/24.0;
            A(15,14) = -S*(C(2,5)+2.0*C(4,5))/24.0;
            A(15,15) = S*(C(2,2)+C(2,4)+C(4,4))/12.0;

            A(16,12) = -S*(-C(0,3)-C(1,3)+2.0*C(0,0)+C(0,1))/24.0;
            A(16,13) = -S*(-C(3,4)+2.0*C(0,4))/24.0;
            A(16,14) = -S*(-2.0*C(3,5)+C(0,5))/24.0;
            A(16,15) = S*(-C(2,3)-2.0*C(3,4)+C(0,2)+C(0,4))/24.0;
            A(16,16) = S*(C(3,3)-C(0,3)+C(0,0))/12.0;

            A(17,12) = S*(C(0,5)+2.0*C(1,5))/24.0;
            A(17,13) = S*C(4,5)/24.0;
            A(17,14) = S*C(5,5)/24.0;
            A(17,15) = -S*(2.0*C(2,5)+C(4,5))/24.0;
            A(17,16) = -S*(-C(3,5)+C(0,5))/24.0;
            A(17,17) = S*C(5,5)/12.0;

            A(18,12) = -S*(2.0*C(0,2)+C(1,2)+C(0,5)+C(1,5))/24.0;
            A(18,13) = -S*(2.0*C(2,4)+C(4,5))/24.0;
            A(18,14) = -S*(C(2,5)+2.0*C(5,5))/24.0;
            A(18,15) = S*(C(2,2)+C(2,4)+C(2,5)+2.0*C(4,5))/24.0;
            A(18,16) = -S*(C(2,3)-2.0*C(0,2)+2.0*C(3,5)-C(0,5))/24.0;
            A(18,17) = -S*(C(2,5)+C(5,5))/24.0;
            A(18,18) = S*(C(2,2)+C(2,5)+C(5,5))/12.0;

            A(19,12) = S*(C(0,2)+2.0*C(1,2))/24.0;
            A(19,13) = S*C(2,4)/24.0;
            A(19,14) = S*C(2,5)/24.0;
            A(19,15) = -S*(2.0*C(2,2)+C(2,4))/24.0;
            A(19,16) = S*(C(2,3)-C(0,2))/24.0;
            A(19,17) = S*C(2,5)/12.0;
            A(19,18) = -S*(C(2,2)+C(2,5))/24.0;
            A(19,19) = S*C(2,2)/12.0;

            A(20,12) = -S*(C(0,1)+2.0*C(1,1)+C(0,3)+C(1,3))/24.0;
            A(20,13) = -S*(C(1,4)+C(3,4))/24.0;
            A(20,14) = -S*(C(1,5)+2.0*C(3,5))/24.0;
            A(20,15) = S*(2.0*C(1,2)+C(1,4)+C(2,3)+2.0*C(3,4))/24.0;
            A(20,16) = S*(-C(1,3)+C(0,1)-2.0*C(3,3)+C(0,3))/24.0;
            A(20,17) = -S*(2.0*C(1,5)+C(3,5))/24.0;
            A(20,18) = S*(C(1,2)+C(1,5)+C(2,3)+2.0*C(3,5))/24.0;
            A(20,19) = -S*(2.0*C(1,2)+C(2,3))/24.0;
            A(20,20) = S*(C(1,1)+C(1,3)+C(3,3))/12.0;

            A(21,12) = -S*(2.0*C(0,4)+C(1,4)+C(0,5)+2.0*C(1,5))/24.0;
            A(21,13) = -S*(2.0*C(4,4)+C(4,5))/24.0;
            A(21,14) = -S*(C(4,5)+C(5,5))/24.0;
            A(21,15) = S*(C(2,4)+C(4,4)+2.0*C(2,5)+C(4,5))/24.0;
            A(21,16) = S*(-C(3,4)+2.0*C(0,4)-C(3,5)+C(0,5))/24.0;
            A(21,17) = -S*(C(4,5)+2.0*C(5,5))/24.0;
            A(21,18) = S*(2.0*C(2,4)+C(4,5)+C(2,5)+C(5,5))/24.0;
            A(21,19) = -S*(C(2,4)+2.0*C(2,5))/24.0;
            A(21,20) = S*(C(1,4)+C(3,4)+2.0*C(1,5)+C(3,5))/24.0;
            A(21,21) = S*(C(4,4)+C(4,5)+C(5,5))/12.0;

            A(22,12) = S*(2.0*C(0,2)+C(1,2))/24.0;
            A(22,13) = S*C(2,4)/12.0;
            A(22,14) = S*C(2,5)/24.0;
            A(22,15) = -S*(C(2,2)+C(2,4))/24.0;
            A(22,16) = S*(C(2,3)-2.0*C(0,2))/24.0;
            A(22,17) = S*C(2,5)/24.0;
            A(22,18) = -S*(2.0*C(2,2)+C(2,5))/24.0;
            A(22,19) = S*C(2,2)/24.0;
            A(22,20) = -S*(C(1,2)+C(2,3))/24.0;
            A(22,21) = -S*(2.0*C(2,4)+C(2,5))/24.0;
            A(22,22) = S*C(2,2)/12.0;

            A(23,12) = S*(C(0,4)+C(1,4))/24.0;
            A(23,13) = S*C(4,4)/24.0;
            A(23,14) = S*C(4,5)/12.0;
            A(23,15) = -S*(C(2,4)+2.0*C(4,4))/24.0;
            A(23,16) = -S*(-2.0*C(3,4)+C(0,4))/24.0;
            A(23,17) = S*C(4,5)/24.0;
            A(23,18) = -S*(C(2,4)+2.0*C(4,5))/24.0;
            A(23,19) = S*C(2,4)/24.0;
            A(23,20) = -S*(C(1,4)+2.0*C(3,4))/24.0;
            A(23,21) = -S*(C(4,4)+C(4,5))/24.0;
            A(23,22) = S*C(2,4)/24.0;
            A(23,23) = S*C(4,4)/12.0;
            break;
        case 1:
            A(0,12) = -S*(C(0,0)+C(0,1))/3.0;
            A(0,13) = -S*(C(0,3)+C(0,4))/3.0;
            A(0,14) = S*C(0,1)/3.0;
            A(0,15) = S*C(0,2)/3.0;
            A(0,16) = S*C(0,0)/3.0;
            A(0,17) = S*(-C(0,5)+C(0,3))/3.0;
            A(0,18) = S*C(0,2)/3.0;
            A(0,19) = -S*(C(0,0)+C(0,2))/3.0;
            A(0,20) = S*C(0,1)/3.0;
            A(0,21) = S*(C(0,4)+C(0,5))/3.0;
            A(0,22) = -S*(C(0,1)+C(0,2))/3.0;
            A(0,23) = S*C(0,0)/3.0;

            A(1,12) = -S*(C(1,0)+C(1,1))/3.0;
            A(1,13) = -S*(C(1,3)+C(1,4))/3.0;
            A(1,14) = S*C(1,1)/3.0;
            A(1,15) = S*C(1,2)/3.0;
            A(1,16) = S*C(1,0)/3.0;
            A(1,17) = S*(-C(1,5)+C(1,3))/3.0;
            A(1,18) = S*C(1,2)/3.0;
            A(1,19) = -S*(C(1,0)+C(1,2))/3.0;
            A(1,20) = S*C(1,1)/3.0;
            A(1,21) = S*(C(1,4)+C(1,5))/3.0;
            A(1,22) = -S*(C(1,1)+C(1,2))/3.0;
            A(1,23) = S*C(1,0)/3.0;

            A(2,12) = -S*(C(2,0)+C(2,1))/3.0;
            A(2,13) = -S*(C(2,3)+C(2,4))/3.0;
            A(2,14) = S*C(2,1)/3.0;
            A(2,15) = S*C(2,2)/3.0;
            A(2,16) = S*C(2,0)/3.0;
            A(2,17) = S*(-C(2,5)+C(2,3))/3.0;
            A(2,18) = S*C(2,2)/3.0;
            A(2,19) = -S*(C(2,0)+C(2,2))/3.0;
            A(2,20) = S*C(2,1)/3.0;
            A(2,21) = S*(C(2,4)+C(2,5))/3.0;
            A(2,22) = -S*(C(2,1)+C(2,2))/3.0;
            A(2,23) = S*C(2,0)/3.0;

            A(3,12) = -S*(C(3,0)+C(3,1))/3.0;
            A(3,13) = -S*(C(3,3)+C(3,4))/3.0;
            A(3,14) = S*C(3,1)/3.0;
            A(3,15) = S*C(3,2)/3.0;
            A(3,16) = S*C(3,0)/3.0;
            A(3,17) = S*(-C(3,5)+C(3,3))/3.0;
            A(3,18) = S*C(3,2)/3.0;
            A(3,19) = -S*(C(3,0)+C(3,2))/3.0;
            A(3,20) = S*C(3,1)/3.0;
            A(3,21) = S*(C(3,4)+C(3,5))/3.0;
            A(3,22) = -S*(C(3,1)+C(3,2))/3.0;
            A(3,23) = S*C(3,0)/3.0;

            A(4,12) = -S*(C(4,0)+C(4,1))/3.0;
            A(4,13) = -S*(C(4,3)+C(4,4))/3.0;
            A(4,14) = S*C(4,1)/3.0;
            A(4,15) = S*C(4,2)/3.0;
            A(4,16) = S*C(4,0)/3.0;
            A(4,17) = S*(-C(4,5)+C(4,3))/3.0;
            A(4,18) = S*C(4,2)/3.0;
            A(4,19) = -S*(C(4,0)+C(4,2))/3.0;
            A(4,20) = S*C(4,1)/3.0;
            A(4,21) = S*(C(4,4)+C(4,5))/3.0;
            A(4,22) = -S*(C(4,1)+C(4,2))/3.0;
            A(4,23) = S*C(4,0)/3.0;

            A(5,12) = -S*(C(5,0)+C(5,1))/3.0;
            A(5,13) = -S*(C(5,3)+C(5,4))/3.0;
            A(5,14) = S*C(5,1)/3.0;
            A(5,15) = S*C(5,2)/3.0;
            A(5,16) = S*C(5,0)/3.0;
            A(5,17) = S*(-C(5,5)+C(5,3))/3.0;
            A(5,18) = S*C(5,2)/3.0;
            A(5,19) = -S*(C(5,0)+C(5,2))/3.0;
            A(5,20) = S*C(5,1)/3.0;
            A(5,21) = S*(C(5,4)+C(5,5))/3.0;
            A(5,22) = -S*(C(5,1)+C(5,2))/3.0;
            A(5,23) = S*C(5,0)/3.0;

            A(12,12) = S*(C(0,0)+C(0,1)+C(1,1))/12.0;

            A(13,12) = S*(C(0,3)+C(1,3)+2.0*C(0,4)+C(1,4))/24.0;
            A(13,13) = S*(C(3,3)+C(3,4)+C(4,4))/12.0;

            A(14,12) = -S*(C(0,1)+C(1,1))/24.0;
            A(14,13) = -S*(2.0*C(1,3)+C(1,4))/24.0;
            A(14,14) = S*C(1,1)/12.0;

            A(15,12) = -S*(C(0,2)+2.0*C(1,2))/24.0;
            A(15,13) = -S*(C(2,3)+C(2,4))/24.0;
            A(15,14) = S*C(1,2)/24.0;
            A(15,15) = S*C(2,2)/12.0;

            A(16,12) = -S*(2.0*C(0,0)+C(0,1))/24.0;
            A(16,13) = -S*(C(0,3)+2.0*C(0,4))/24.0;
            A(16,14) = S*C(0,1)/24.0;
            A(16,15) = S*C(0,2)/24.0;
            A(16,16) = S*C(0,0)/12.0;

            A(17,12) = -S*(-C(0,5)-2.0*C(1,5)+C(0,3)+C(1,3))/24.0;
            A(17,13) = -S*(-C(3,5)-C(4,5)+2.0*C(3,3)+C(3,4))/24.0;
            A(17,14) = S*(-C(1,5)+2.0*C(1,3))/24.0;
            A(17,15) = S*(-2.0*C(2,5)+C(2,3))/24.0;
            A(17,16) = S*(-C(0,5)+C(0,3))/24.0;
            A(17,17) = S*(C(5,5)-C(3,5)+C(3,3))/12.0;

            A(18,12) = -S*(2.0*C(0,2)+C(1,2))/24.0;
            A(18,13) = -S*(C(2,3)+2.0*C(2,4))/24.0;
            A(18,14) = S*C(1,2)/24.0;
            A(18,15) = S*C(2,2)/24.0;
            A(18,16) = S*C(0,2)/12.0;
            A(18,17) = S*(-C(2,5)+C(2,3))/24.0;
            A(18,18) = S*C(2,2)/12.0;

            A(19,12) = S*(C(0,0)+C(0,1)+C(0,2)+2.0*C(1,2))/24.0;
            A(19,13) = S*(2.0*C(0,3)+C(0,4)+C(2,3)+C(2,4))/24.0;
            A(19,14) = -S*(2.0*C(0,1)+C(1,2))/24.0;
            A(19,15) = -S*(C(0,2)+2.0*C(2,2))/24.0;
            A(19,16) = -S*(C(0,0)+C(0,2))/24.0;
            A(19,17) = -S*(-C(0,5)+2.0*C(0,3)-2.0*C(2,5)+C(2,3))/24.0;
            A(19,18) = -S*(C(0,2)+C(2,2))/24.0;
            A(19,19) = S*(C(0,0)+C(0,2)+C(2,2))/12.0;

            A(20,12) = -S*(C(0,1)+2.0*C(1,1))/24.0;
            A(20,13) = -S*(C(1,3)+C(1,4))/24.0;
            A(20,14) = S*C(1,1)/24.0;
            A(20,15) = S*C(1,2)/12.0;
            A(20,16) = S*C(0,1)/24.0;
            A(20,17) = S*(-2.0*C(1,5)+C(1,3))/24.0;
            A(20,18) = S*C(1,2)/24.0;
            A(20,19) = -S*(C(0,1)+2.0*C(1,2))/24.0;
            A(20,20) = S*C(1,1)/12.0;

            A(21,12) = -S*(2.0*C(0,4)+C(1,4)+C(0,5)+2.0*C(1,5))/24.0;
            A(21,13) = -S*(C(3,4)+2.0*C(4,4)+C(3,5)+C(4,5))/24.0;
            A(21,14) = S*(C(1,4)+C(1,5))/24.0;
            A(21,15) = S*(C(2,4)+2.0*C(2,5))/24.0;
            A(21,16) = S*(2.0*C(0,4)+C(0,5))/24.0;
            A(21,17) = S*(-C(4,5)+C(3,4)-2.0*C(5,5)+C(3,5))/24.0;
            A(21,18) = S*(2.0*C(2,4)+C(2,5))/24.0;
            A(21,19) = -S*(C(0,4)+C(2,4)+C(0,5)+2.0*C(2,5))/24.0;
            A(21,20) = S*(C(1,4)+2.0*C(1,5))/24.0;
            A(21,21) = S*(C(4,4)+C(4,5)+C(5,5))/12.0;

            A(22,12) = S*(C(0,1)+C(1,1)+2.0*C(0,2)+C(1,2))/24.0;
            A(22,13) = S*(2.0*C(1,3)+C(1,4)+C(2,3)+2.0*C(2,4))/24.0;
            A(22,14) = -S*(2.0*C(1,1)+C(1,2))/24.0;
            A(22,15) = -S*(C(1,2)+C(2,2))/24.0;
            A(22,16) = -S*(C(0,1)+2.0*C(0,2))/24.0;
            A(22,17) = -S*(-C(1,5)+2.0*C(1,3)-C(2,5)+C(2,3))/24.0;
            A(22,18) = -S*(C(1,2)+2.0*C(2,2))/24.0;
            A(22,19) = S*(2.0*C(0,1)+C(1,2)+C(0,2)+C(2,2))/24.0;
            A(22,20) = -S*(C(1,1)+C(1,2))/24.0;
            A(22,21) = -S*(C(1,4)+C(1,5)+2.0*C(2,4)+C(2,5))/24.0;
            A(22,22) = S*(C(1,1)+C(1,2)+C(2,2))/12.0;

            A(23,12) = -S*(C(0,0)+C(0,1))/24.0;
            A(23,13) = -S*(2.0*C(0,3)+C(0,4))/24.0;
            A(23,14) = S*C(0,1)/12.0;
            A(23,15) = S*C(0,2)/24.0;
            A(23,16) = S*C(0,0)/24.0;
            A(23,17) = S*(-C(0,5)+2.0*C(0,3))/24.0;
            A(23,18) = S*C(0,2)/24.0;
            A(23,19) = -S*(2.0*C(0,0)+C(0,2))/24.0;
            A(23,20) = S*C(0,1)/24.0;
            A(23,21) = S*(C(0,4)+C(0,5))/24.0;
            A(23,22) = -S*(2.0*C(0,1)+C(0,2))/24.0;
            A(23,23) = S*C(0,0)/12.0;
            break;
        case 2:
            A(0,12) = -S*C(0,0)/3.0;
            A(0,13) = -S*(C(0,3)+C(0,4))/3.0;
            A(0,14) = S*(-C(0,5)+C(0,1))/3.0;
            A(0,15) = S*C(0,4)/3.0;
            A(0,16) = S*(-C(0,3)+C(0,0))/3.0;
            A(0,17) = S*C(0,3)/3.0;
            A(0,18) = S*(C(0,2)+C(0,5))/3.0;
            A(0,19) = -S*C(0,0)/3.0;
            A(0,20) = S*C(0,3)/3.0;
            A(0,21) = S*C(0,4)/3.0;
            A(0,22) = -S*(C(0,1)+C(0,2))/3.0;
            A(0,23) = S*(-C(0,4)+C(0,0))/3.0;

            A(1,12) = -S*C(1,0)/3.0;
            A(1,13) = -S*(C(1,3)+C(1,4))/3.0;
            A(1,14) = S*(-C(1,5)+C(1,1))/3.0;
            A(1,15) = S*C(1,4)/3.0;
            A(1,16) = S*(-C(1,3)+C(1,0))/3.0;
            A(1,17) = S*C(1,3)/3.0;
            A(1,18) = S*(C(1,2)+C(1,5))/3.0;
            A(1,19) = -S*C(1,0)/3.0;
            A(1,20) = S*C(1,3)/3.0;
            A(1,21) = S*C(1,4)/3.0;
            A(1,22) = -S*(C(1,1)+C(1,2))/3.0;
            A(1,23) = S*(-C(1,4)+C(1,0))/3.0;

            A(2,12) = -S*C(2,0)/3.0;
            A(2,13) = -S*(C(2,3)+C(2,4))/3.0;
            A(2,14) = S*(-C(2,5)+C(2,1))/3.0;
            A(2,15) = S*C(2,4)/3.0;
            A(2,16) = S*(-C(2,3)+C(2,0))/3.0;
            A(2,17) = S*C(2,3)/3.0;
            A(2,18) = S*(C(2,2)+C(2,5))/3.0;
            A(2,19) = -S*C(2,0)/3.0;
            A(2,20) = S*C(2,3)/3.0;
            A(2,21) = S*C(2,4)/3.0;
            A(2,22) = -S*(C(2,1)+C(2,2))/3.0;
            A(2,23) = S*(-C(2,4)+C(2,0))/3.0;

            A(3,12) = -S*C(3,0)/3.0;
            A(3,13) = -S*(C(3,3)+C(3,4))/3.0;
            A(3,14) = S*(-C(3,5)+C(3,1))/3.0;
            A(3,15) = S*C(3,4)/3.0;
            A(3,16) = S*(-C(3,3)+C(3,0))/3.0;
            A(3,17) = S*C(3,3)/3.0;
            A(3,18) = S*(C(3,2)+C(3,5))/3.0;
            A(3,19) = -S*C(3,0)/3.0;
            A(3,20) = S*C(3,3)/3.0;
            A(3,21) = S*C(3,4)/3.0;
            A(3,22) = -S*(C(3,1)+C(3,2))/3.0;
            A(3,23) = S*(-C(3,4)+C(3,0))/3.0;

            A(4,12) = -S*C(4,0)/3.0;
            A(4,13) = -S*(C(4,3)+C(4,4))/3.0;
            A(4,14) = S*(-C(4,5)+C(4,1))/3.0;
            A(4,15) = S*C(4,4)/3.0;
            A(4,16) = S*(-C(4,3)+C(4,0))/3.0;
            A(4,17) = S*C(4,3)/3.0;
            A(4,18) = S*(C(4,2)+C(4,5))/3.0;
            A(4,19) = -S*C(4,0)/3.0;
            A(4,20) = S*C(4,3)/3.0;
            A(4,21) = S*C(4,4)/3.0;
            A(4,22) = -S*(C(4,1)+C(4,2))/3.0;
            A(4,23) = S*(-C(4,4)+C(4,0))/3.0;

            A(5,12) = -S*C(5,0)/3.0;
            A(5,13) = -S*(C(5,3)+C(5,4))/3.0;
            A(5,14) = S*(-C(5,5)+C(5,1))/3.0;
            A(5,15) = S*C(5,4)/3.0;
            A(5,16) = S*(-C(5,3)+C(5,0))/3.0;
            A(5,17) = S*C(5,3)/3.0;
            A(5,18) = S*(C(5,2)+C(5,5))/3.0;
            A(5,19) = -S*C(5,0)/3.0;
            A(5,20) = S*C(5,3)/3.0;
            A(5,21) = S*C(5,4)/3.0;
            A(5,22) = -S*(C(5,1)+C(5,2))/3.0;
            A(5,23) = S*(-C(5,4)+C(5,0))/3.0;

            A(12,12) = S*C(0,0)/12.0;

            A(13,12) = S*(C(0,3)+2.0*C(0,4))/24.0;
            A(13,13) = S*(C(3,3)+C(3,4)+C(4,4))/12.0;

            A(14,12) = -S*(-C(0,5)+C(0,1))/24.0;
            A(14,13) = -S*(-C(3,5)-C(4,5)+2.0*C(1,3)+C(1,4))/24.0;
            A(14,14) = S*(C(5,5)-C(1,5)+C(1,1))/12.0;

            A(15,12) = -S*C(0,4)/24.0;
            A(15,13) = -S*(C(3,4)+C(4,4))/24.0;
            A(15,14) = S*(-2.0*C(4,5)+C(1,4))/24.0;
            A(15,15) = S*C(4,4)/12.0;

            A(16,12) = -S*(-C(0,3)+2.0*C(0,0))/24.0;
            A(16,13) = -S*(-C(3,3)-C(3,4)+C(0,3)+2.0*C(0,4))/24.0;
            A(16,14) = S*(2.0*C(3,5)-C(1,3)-C(0,5)+C(0,1))/24.0;
            A(16,15) = S*(-2.0*C(3,4)+C(0,4))/24.0;
            A(16,16) = S*(C(3,3)-C(0,3)+C(0,0))/12.0;

            A(17,12) = -S*C(0,3)/24.0;
            A(17,13) = -S*(2.0*C(3,3)+C(3,4))/24.0;
            A(17,14) = S*(-C(3,5)+2.0*C(1,3))/24.0;
            A(17,15) = S*C(3,4)/24.0;
            A(17,16) = S*(-C(3,3)+C(0,3))/24.0;
            A(17,17) = S*C(3,3)/12.0;

            A(18,12) = -S*(2.0*C(0,2)+C(0,5))/24.0;
            A(18,13) = -S*(C(2,3)+2.0*C(2,4)+C(3,5)+C(4,5))/24.0;
            A(18,14) = S*(-C(2,5)+C(1,2)-2.0*C(5,5)+C(1,5))/24.0;
            A(18,15) = S*(C(2,4)+2.0*C(4,5))/24.0;
            A(18,16) = -S*(C(2,3)-2.0*C(0,2)+2.0*C(3,5)-C(0,5))/24.0;
            A(18,17) = S*(C(2,3)+C(3,5))/24.0;
            A(18,18) = S*(C(2,2)+C(2,5)+C(5,5))/12.0;

            A(19,12) = S*C(0,0)/24.0;
            A(19,13) = S*(2.0*C(0,3)+C(0,4))/24.0;
            A(19,14) = -S*(-C(0,5)+2.0*C(0,1))/24.0;
            A(19,15) = -S*C(0,4)/24.0;
            A(19,16) = -S*(-C(0,3)+C(0,0))/24.0;
            A(19,17) = -S*C(0,3)/12.0;
            A(19,18) = -S*(C(0,2)+C(0,5))/24.0;
            A(19,19) = S*C(0,0)/12.0;

            A(20,12) = -S*C(0,3)/24.0;
            A(20,13) = -S*(C(3,3)+C(3,4))/24.0;
            A(20,14) = S*(-2.0*C(3,5)+C(1,3))/24.0;
            A(20,15) = S*C(3,4)/12.0;
            A(20,16) = S*(-2.0*C(3,3)+C(0,3))/24.0;
            A(20,17) = S*C(3,3)/24.0;
            A(20,18) = S*(C(2,3)+2.0*C(3,5))/24.0;
            A(20,19) = -S*C(0,3)/24.0;
            A(20,20) = S*C(3,3)/12.0;

            A(21,12) = -S*C(0,4)/12.0;
            A(21,13) = -S*(C(3,4)+2.0*C(4,4))/24.0;
            A(21,14) = S*(-C(4,5)+C(1,4))/24.0;
            A(21,15) = S*C(4,4)/24.0;
            A(21,16) = S*(-C(3,4)+2.0*C(0,4))/24.0;
            A(21,17) = S*C(3,4)/24.0;
            A(21,18) = S*(2.0*C(2,4)+C(4,5))/24.0;
            A(21,19) = -S*C(0,4)/24.0;
            A(21,20) = S*C(3,4)/24.0;
            A(21,21) = S*C(4,4)/12.0;

            A(22,12) = S*(C(0,1)+2.0*C(0,2))/24.0;
            A(22,13) = S*(2.0*C(1,3)+C(1,4)+C(2,3)+2.0*C(2,4))/24.0;
            A(22,14) = -S*(-C(1,5)+2.0*C(1,1)-C(2,5)+C(1,2))/24.0;
            A(22,15) = -S*(C(1,4)+C(2,4))/24.0;
            A(22,16) = -S*(-C(1,3)+C(0,1)-C(2,3)+2.0*C(0,2))/24.0;
            A(22,17) = -S*(2.0*C(1,3)+C(2,3))/24.0;
            A(22,18) = -S*(C(1,2)+C(1,5)+2.0*C(2,2)+C(2,5))/24.0;
            A(22,19) = S*(2.0*C(0,1)+C(0,2))/24.0;
            A(22,20) = -S*(C(1,3)+C(2,3))/24.0;
            A(22,21) = -S*(C(1,4)+2.0*C(2,4))/24.0;
            A(22,22) = S*(C(1,1)+C(1,2)+C(2,2))/12.0;

            A(23,12) = -S*(-C(0,4)+C(0,0))/24.0;
            A(23,13) = -S*(-C(3,4)-C(4,4)+2.0*C(0,3)+C(0,4))/24.0;
            A(23,14) = S*(2.0*C(4,5)-C(1,4)-C(0,5)+2.0*C(0,1))/24.0;
            A(23,15) = S*(-2.0*C(4,4)+C(0,4))/24.0;
            A(23,16) = S*(2.0*C(3,4)-C(0,4)-C(0,3)+C(0,0))/24.0;
            A(23,17) = S*(-C(3,4)+2.0*C(0,3))/24.0;
            A(23,18) = -S*(C(2,4)+2.0*C(4,5)-C(0,2)-C(0,5))/24.0;
            A(23,19) = -S*(-C(0,4)+2.0*C(0,0))/24.0;
            A(23,20) = S*(-2.0*C(3,4)+C(0,3))/24.0;
            A(23,21) = S*(-C(4,4)+C(0,4))/24.0;
            A(23,22) = -S*(-C(1,4)-C(2,4)+2.0*C(0,1)+C(0,2))/24.0;
            A(23,23) = S*(C(4,4)-C(0,4)+C(0,0))/12.0;
            break;
        case 3:
            A(0,12) = -S*C(0,1)/3.0;
            A(0,13) = -S*C(0,3)/3.0;
            A(0,14) = S*(-C(0,5)+C(0,1))/3.0;
            A(0,15) = S*(C(0,2)+C(0,4))/3.0;
            A(0,16) = -S*C(0,3)/3.0;
            A(0,17) = S*(-C(0,5)+C(0,3))/3.0;
            A(0,18) = S*C(0,5)/3.0;
            A(0,19) = -S*(C(0,0)+C(0,2))/3.0;
            A(0,20) = S*(C(0,1)+C(0,3))/3.0;
            A(0,21) = S*C(0,5)/3.0;
            A(0,22) = -S*C(0,1)/3.0;
            A(0,23) = S*(-C(0,4)+C(0,0))/3.0;

            A(1,12) = -S*C(1,1)/3.0;
            A(1,13) = -S*C(1,3)/3.0;
            A(1,14) = S*(-C(1,5)+C(1,1))/3.0;
            A(1,15) = S*(C(1,2)+C(1,4))/3.0;
            A(1,16) = -S*C(1,3)/3.0;
            A(1,17) = S*(-C(1,5)+C(1,3))/3.0;
            A(1,18) = S*C(1,5)/3.0;
            A(1,19) = -S*(C(1,0)+C(1,2))/3.0;
            A(1,20) = S*(C(1,1)+C(1,3))/3.0;
            A(1,21) = S*C(1,5)/3.0;
            A(1,22) = -S*C(1,1)/3.0;
            A(1,23) = S*(-C(1,4)+C(1,0))/3.0;

            A(2,12) = -S*C(2,1)/3.0;
            A(2,13) = -S*C(2,3)/3.0;
            A(2,14) = S*(-C(2,5)+C(2,1))/3.0;
            A(2,15) = S*(C(2,2)+C(2,4))/3.0;
            A(2,16) = -S*C(2,3)/3.0;
            A(2,17) = S*(-C(2,5)+C(2,3))/3.0;
            A(2,18) = S*C(2,5)/3.0;
            A(2,19) = -S*(C(2,0)+C(2,2))/3.0;
            A(2,20) = S*(C(2,1)+C(2,3))/3.0;
            A(2,21) = S*C(2,5)/3.0;
            A(2,22) = -S*C(2,1)/3.0;
            A(2,23) = S*(-C(2,4)+C(2,0))/3.0;

            A(3,12) = -S*C(3,1)/3.0;
            A(3,13) = -S*C(3,3)/3.0;
            A(3,14) = S*(-C(3,5)+C(3,1))/3.0;
            A(3,15) = S*(C(3,2)+C(3,4))/3.0;
            A(3,16) = -S*C(3,3)/3.0;
            A(3,17) = S*(-C(3,5)+C(3,3))/3.0;
            A(3,18) = S*C(3,5)/3.0;
            A(3,19) = -S*(C(3,0)+C(3,2))/3.0;
            A(3,20) = S*(C(3,1)+C(3,3))/3.0;
            A(3,21) = S*C(3,5)/3.0;
            A(3,22) = -S*C(3,1)/3.0;
            A(3,23) = S*(-C(3,4)+C(3,0))/3.0;

            A(4,12) = -S*C(4,1)/3.0;
            A(4,13) = -S*C(4,3)/3.0;
            A(4,14) = S*(-C(4,5)+C(4,1))/3.0;
            A(4,15) = S*(C(4,2)+C(4,4))/3.0;
            A(4,16) = -S*C(4,3)/3.0;
            A(4,17) = S*(-C(4,5)+C(4,3))/3.0;
            A(4,18) = S*C(4,5)/3.0;
            A(4,19) = -S*(C(4,0)+C(4,2))/3.0;
            A(4,20) = S*(C(4,1)+C(4,3))/3.0;
            A(4,21) = S*C(4,5)/3.0;
            A(4,22) = -S*C(4,1)/3.0;
            A(4,23) = S*(-C(4,4)+C(4,0))/3.0;

            A(5,12) = -S*C(5,1)/3.0;
            A(5,13) = -S*C(5,3)/3.0;
            A(5,14) = S*(-C(5,5)+C(5,1))/3.0;
            A(5,15) = S*(C(5,2)+C(5,4))/3.0;
            A(5,16) = -S*C(5,3)/3.0;
            A(5,17) = S*(-C(5,5)+C(5,3))/3.0;
            A(5,18) = S*C(5,5)/3.0;
            A(5,19) = -S*(C(5,0)+C(5,2))/3.0;
            A(5,20) = S*(C(5,1)+C(5,3))/3.0;
            A(5,21) = S*C(5,5)/3.0;
            A(5,22) = -S*C(5,1)/3.0;
            A(5,23) = S*(-C(5,4)+C(5,0))/3.0;

            A(12,12) = S*C(1,1)/12.0;

            A(13,12) = S*C(1,3)/24.0;
            A(13,13) = S*C(3,3)/12.0;

            A(14,12) = -S*(-C(1,5)+C(1,1))/24.0;
            A(14,13) = -S*(-C(3,5)+2.0*C(1,3))/24.0;
            A(14,14) = S*(C(5,5)-C(1,5)+C(1,1))/12.0;

            A(15,12) = -S*(2.0*C(1,2)+C(1,4))/24.0;
            A(15,13) = -S*(C(2,3)+C(3,4))/24.0;
            A(15,14) = S*(-C(2,5)-2.0*C(4,5)+C(1,2)+C(1,4))/24.0;
            A(15,15) = S*(C(2,2)+C(2,4)+C(4,4))/12.0;

            A(16,12) = S*C(1,3)/24.0;
            A(16,13) = S*C(3,3)/24.0;
            A(16,14) = -S*(-2.0*C(3,5)+C(1,3))/24.0;
            A(16,15) = -S*(C(2,3)+2.0*C(3,4))/24.0;
            A(16,16) = S*C(3,3)/12.0;

            A(17,12) = -S*(-2.0*C(1,5)+C(1,3))/24.0;
            A(17,13) = -S*(-C(3,5)+2.0*C(3,3))/24.0;
            A(17,14) = S*(C(5,5)-C(1,5)-C(3,5)+2.0*C(1,3))/24.0;
            A(17,15) = S*(-2.0*C(2,5)-C(4,5)+C(2,3)+C(3,4))/24.0;
            A(17,16) = -S*(-C(3,5)+C(3,3))/24.0;
            A(17,17) = S*(C(5,5)-C(3,5)+C(3,3))/12.0;

            A(18,12) = -S*C(1,5)/24.0;
            A(18,13) = -S*C(3,5)/24.0;
            A(18,14) = S*(-2.0*C(5,5)+C(1,5))/24.0;
            A(18,15) = S*(C(2,5)+2.0*C(4,5))/24.0;
            A(18,16) = -S*C(3,5)/12.0;
            A(18,17) = S*(-C(5,5)+C(3,5))/24.0;
            A(18,18) = S*C(5,5)/12.0;

            A(19,12) = S*(C(0,1)+2.0*C(1,2))/24.0;
            A(19,13) = S*(2.0*C(0,3)+C(2,3))/24.0;
            A(19,14) = -S*(-C(0,5)+2.0*C(0,1)-C(2,5)+C(1,2))/24.0;
            A(19,15) = -S*(C(0,2)+C(0,4)+2.0*C(2,2)+C(2,4))/24.0;
            A(19,16) = S*(C(0,3)+C(2,3))/24.0;
            A(19,17) = -S*(-C(0,5)+2.0*C(0,3)-2.0*C(2,5)+C(2,3))/24.0;
            A(19,18) = -S*(C(0,5)+C(2,5))/24.0;
            A(19,19) = S*(C(0,0)+C(0,2)+C(2,2))/12.0;

            A(20,12) = -S*(2.0*C(1,1)+C(1,3))/24.0;
            A(20,13) = -S*(C(1,3)+C(3,3))/24.0;
            A(20,14) = S*(-C(1,5)+C(1,1)-2.0*C(3,5)+C(1,3))/24.0;
            A(20,15) = S*(2.0*C(1,2)+C(1,4)+C(2,3)+2.0*C(3,4))/24.0;
            A(20,16) = -S*(C(1,3)+2.0*C(3,3))/24.0;
            A(20,17) = S*(-2.0*C(1,5)+C(1,3)-C(3,5)+C(3,3))/24.0;
            A(20,18) = S*(C(1,5)+2.0*C(3,5))/24.0;
            A(20,19) = -S*(C(0,1)+2.0*C(1,2)+C(0,3)+C(2,3))/24.0;
            A(20,20) = S*(C(1,1)+C(1,3)+C(3,3))/12.0;

            A(21,12) = -S*C(1,5)/12.0;
            A(21,13) = -S*C(3,5)/24.0;
            A(21,14) = S*(-C(5,5)+C(1,5))/24.0;
            A(21,15) = S*(2.0*C(2,5)+C(4,5))/24.0;
            A(21,16) = -S*C(3,5)/24.0;
            A(21,17) = S*(-2.0*C(5,5)+C(3,5))/24.0;
            A(21,18) = S*C(5,5)/24.0;
            A(21,19) = -S*(C(0,5)+2.0*C(2,5))/24.0;
            A(21,20) = S*(2.0*C(1,5)+C(3,5))/24.0;
            A(21,21) = S*C(5,5)/12.0;

            A(22,12) = S*C(1,1)/24.0;
            A(22,13) = S*C(1,3)/12.0;
            A(22,14) = -S*(-C(1,5)+2.0*C(1,1))/24.0;
            A(22,15) = -S*(C(1,2)+C(1,4))/24.0;
            A(22,16) = S*C(1,3)/24.0;
            A(22,17) = -S*(-C(1,5)+2.0*C(1,3))/24.0;
            A(22,18) = -S*C(1,5)/24.0;
            A(22,19) = S*(2.0*C(0,1)+C(1,2))/24.0;
            A(22,20) = -S*(C(1,1)+C(1,3))/24.0;
            A(22,21) = -S*C(1,5)/24.0;
            A(22,22) = S*C(1,1)/12.0;

            A(23,12) = -S*(-C(1,4)+C(0,1))/24.0;
            A(23,13) = -S*(-C(3,4)+2.0*C(0,3))/24.0;
            A(23,14) = S*(2.0*C(4,5)-C(1,4)-C(0,5)+2.0*C(0,1))/24.0;
            A(23,15) = S*(-C(2,4)-2.0*C(4,4)+C(0,2)+C(0,4))/24.0;
            A(23,16) = -S*(-2.0*C(3,4)+C(0,3))/24.0;
            A(23,17) = S*(C(4,5)-C(3,4)-C(0,5)+2.0*C(0,3))/24.0;
            A(23,18) = S*(-2.0*C(4,5)+C(0,5))/24.0;
            A(23,19) = -S*(-C(0,4)-C(2,4)+2.0*C(0,0)+C(0,2))/24.0;
            A(23,20) = S*(-C(1,4)-2.0*C(3,4)+C(0,1)+C(0,3))/24.0;
            A(23,21) = S*(-C(4,5)+C(0,5))/24.0;
            A(23,22) = -S*(-C(1,4)+2.0*C(0,1))/24.0;
            A(23,23) = S*(C(4,4)-C(0,4)+C(0,0))/12.0;
            break;
        } // switch(face_id_local)
    } // order == 2

    // edge orientation
    Vector<double> ddofs(dofs._n);
    ddofs.inject(dofs);
    for (int k = 0; k < 6; ++ k) {
        if (!tet->get_edge_orientation(k))
            ddofs(k) = -ddofs(k);    
    }
    
    // evaluate integral
    double I = 0.0;
    for (int i = 0; i < 6; ++ i) {
        for (int j = 0; j < i; ++ j) {
            I += 2*ddofs(i)*ddofs(j)*A(i,j);
        }
        I += ddofs(i)*ddofs(i)*A(i,i);
    }

    if (get_order() == 2) {
        for (int i = 0; i < 6; ++ i)
            for (int j = 12; j < 24; ++ j)
                I += 2*ddofs(i)*ddofs(j)*A(i,j);
        for (int i = 12; i < 24; ++ i) {
            for (int j = 12; j < i; ++ j) {
                I += 2*ddofs(i)*ddofs(j)*A(i,j);
            }
            I += ddofs(i)*ddofs(i)*A(i,i);
        }
    }

    return I;
}

---