Skip to content

Instantly share code, notes, and snippets.

@GaZ3ll3

GaZ3ll3/Poisson.h

Last active September 6, 2016 05:05
Show Gist options
  • Save GaZ3ll3/d508c1ba31a6327247ae to your computer and use it in GitHub Desktop.
Save GaZ3ll3/d508c1ba31a6327247ae to your computer and use it in GitHub Desktop.
Download ZIP
mex fenics with matlab
Raw
makefile
CXX = g++
CC = gcc
FF = gfortran
Opt = -Ofast
MATLAB_ROOT = /usr/local/MATLAB/R2014b/
FILE_NAME = test
CXX_FLAGS = -DMATLAB_MEX_FILE -std=c++11 -fopenmp -march=native \
-D_GNU_SOURCE -fPIC -fno-omit-frame-pointer -Wno-write-strings -pthread\
-DMX_COMPAT_32 $(Opt) -DNDEBUG -fopenmp -ffast-math \
-DDOLFIN_LA_INDEX_SIZE=4 -DDOLFIN_SIZE_T=8 \
-DDOLFIN_VERSION=\"1.6.0\" -DHAS_CHOLMOD -DHAS_HDF5 \
-DHAS_MPI -DHAS_OPENMP -DHAS_PETSC -DHAS_QT4 -DHAS_QVTK \
-DHAS_SCOTCH -DHAS_SLEPC -DHAS_UMFPACK -DHAS_VTK -DHAS_ZLIB \
-DNDEBUG -D_BSD_SOURCE -D_FORTIFY_SOURCE=2 -D_LARGEFILE64_SOURCE \
-D_LARGEFILE_SOURCE -g -O2 -fstack-protector --param=ssp-buffer-size=4 \
-Wformat -Werror=format-security -D_FORTIFY_SOURCE=2 \
-Wno-deprecated
CXX_INCLUDE = -I./include/mexplus \
-I./include/pprint \
-I./include/Utility \
-I./include/exprtk \
-I/usr/include/eigen3\
-I$(MATLAB_ROOT)extern/include \
-I$(MATLAB_ROOT)simulink/include \
-isystem /usr/lib/slepcdir/3.4.2 \
-isystem /usr/lib/slepcdir/3.4.2/linux-gnu-c-opt/include\
-isystem /usr/lib/slepcdir/3.4.2/include\
-isystem /usr/include/suitesparse \
-isystem /usr/include/scotch \
-isystem /usr/lib/openmpi/include\
-isystem /usr/lib/openmpi/include/openmpi \
-isystem /usr/include/eigen3 \
-isystem /usr/lib/petscdir/3.4.2/include \
-isystem /usr/lib/petscdir/3.4.2/linux-gnu-c-opt/include \
-isystem /usr/include/qt4 \
-isystem /usr/include/vtk-5.8
MATLAB_LINKS = $(Opt) -pthread -shared\
-Wl,--version-script,$(MATLAB_ROOT)extern/lib/glnxa64/mexFunction.map \
-Wl,--no-undefined
CXX_LIBS = -Wl,-rpath-link,$(MATLAB_ROOT)bin/glnxa64 \
-L$(MATLAB_ROOT)bin/glnxa64 -lmx -lmex -lmat -lm -fopenmp
EXTRA_LINKS = -g -O2 -fstack-protector \
--param=ssp-buffer-size=4 -Wformat -Werror=format-security \
-D_FORTIFY_SOURCE=2 -std=c++11 -Wno-deprecated -fopenmp \
-O2 -g -DNDEBUG -rdynamic /usr/lib/x86_64-linux-gnu/libdolfin.so.1.6.0\
-lboost_filesystem -lboost_program_options -lboost_system -lboost_thread \
-lboost_iostreams -lboost_timer -lpthread -lhdf5 -lpthread -lhdf5 -lz -ldl \
-lm /usr/lib/slepcdir/3.4.2/linux-gnu-c-opt/lib/libslepc.so \
/usr/lib/petscdir/3.4.2/linux-gnu-c-opt/lib/libpetsc.so -lumfpack \
-lamd -lblas -lcholmod -lcamd -lcolamd -lccolamd -Wl,-Bstatic \
-lsuitesparseconfig -Wl,-Bdynamic -lrt -llapack -lgfortran \
-lptscotch -lptesmumps -lptscotcherr -lmpi_cxx -lmpi -lz -ldl\
-lm /usr/lib/slepcdir/3.4.2/linux-gnu-c-opt/lib/libslepc.so \
/usr/lib/petscdir/3.4.2/linux-gnu-c-opt/lib/libpetsc.so -lumfpack \
-lamd -lblas -lcholmod -lcamd -lcolamd -lccolamd -Wl,-Bstatic \
-lsuitesparseconfig -Wl,-Bdynamic -lrt -llapack -lgfortran -lptscotch \
-lptesmumps -lptscotcherr -lmpi_cxx -lmpi -lhwloc -lQtGui -lQtCore \
/usr/lib/libvtkGenericFiltering.so.5.8.0 /usr/lib/libvtkGeovis.so.5.8.0 \
/usr/lib/libvtkCharts.so.5.8.0 /usr/lib/libvtkViews.so.5.8.0 \
/usr/lib/libvtkInfovis.so.5.8.0 /usr/lib/libvtkWidgets.so.5.8.0 \
/usr/lib/libvtkVolumeRendering.so.5.8.0 /usr/lib/libvtkHybrid.so.5.8.0 \
/usr/lib/libvtkParallel.so.5.8.0 /usr/lib/libvtkRendering.so.5.8.0 \
/usr/lib/libvtkImaging.so.5.8.0 /usr/lib/libvtkGraphics.so.5.8.0 \
/usr/lib/libvtkIO.so.5.8.0 /usr/lib/libvtkFiltering.so.5.8.0 \
/usr/lib/libvtkCommon.so.5.8.0 -lm /usr/lib/libvtksys.so.5.8.0 \
-ldl -lQVTK -Wl,-rpath,/usr/lib/slepcdir/3.4.2/linux-gnu-c-opt/lib:/usr/lib/petscdir/3.4.2/linux-gnu-c-opt/lib
all:
$(CXX) $(CXX_FLAGS) $(CXX_INCLUDE) $(FILE_NAME).cpp $(MATLAB_LINKS) -o $(FILE_NAME).mexa64 $(CXX_LIBS) $(EXTRA_LINKS)
Raw
Poisson.h
// This code conforms with the UFC specification version 1.6.0
// and was automatically generated by FFC version 1.6.0.
//
// This code was generated with the option '-l dolfin' and
// contains DOLFIN-specific wrappers that depend on DOLFIN.
//
// This code was generated with the following parameters:
//
// convert_exceptions_to_warnings: False
// cpp_optimize: True
// cpp_optimize_flags: '-O2'
// epsilon: 1e-14
// error_control: False
// form_postfix: True
// format: 'dolfin'
// optimize: True
// precision: 15
// quadrature_degree: 'auto'
// quadrature_rule: 'auto'
// representation: 'auto'
// restrict_keyword: ''
// split: False
#ifndef __POISSON_H
#define __POISSON_H
#include <cmath>
#include <stdexcept>
#include <fstream>
#include <ufc.h>
/// This class defines the interface for a finite element.
class poisson_finite_element_0: public ufc::finite_element
{
public:
/// Constructor
poisson_finite_element_0() : ufc::finite_element()
{
// Do nothing
}
/// Destructor
virtual ~poisson_finite_element_0()
{
// Do nothing
}
/// Return a string identifying the finite element
virtual const char* signature() const
{
return "FiniteElement('Lagrange', Domain(Cell('triangle', 2)), 1, None)";
}
/// Return the cell shape
virtual ufc::shape cell_shape() const
{
return ufc::triangle;
}
/// Return the topological dimension of the cell shape
virtual std::size_t topological_dimension() const
{
return 2;
}
/// Return the geometric dimension of the cell shape
virtual std::size_t geometric_dimension() const
{
return 2;
}
/// Return the dimension of the finite element function space
virtual std::size_t space_dimension() const
{
return 3;
}
/// Return the rank of the value space
virtual std::size_t value_rank() const
{
return 0;
}
/// Return the dimension of the value space for axis i
virtual std::size_t value_dimension(std::size_t i) const
{
return 1;
}
/// Evaluate basis function i at given point x in cell (actual implementation)
static void _evaluate_basis(std::size_t i,
double* values,
const double* x,
const double* vertex_coordinates,
int cell_orientation)
{
// Compute Jacobian
double J[4];
compute_jacobian_triangle_2d(J, vertex_coordinates);
// Compute Jacobian inverse and determinant
double K[4];
double detJ;
compute_jacobian_inverse_triangle_2d(K, detJ, J);
// Compute constants
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
// Get coordinates and map to the reference (FIAT) element
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
// Reset values
*values = 0.0;
switch (i)
{
case 0:
{
// Array of basisvalues
double basisvalues[3] = {0.0, 0.0, 0.0};
// Declare helper variables
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
// Compute basisvalues
basisvalues[0] = 1.0;
basisvalues[1] = tmp0;
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
basisvalues[0] *= std::sqrt(0.5);
basisvalues[2] *= std::sqrt(1.0);
basisvalues[1] *= std::sqrt(3.0);
// Table(s) of coefficients
static const double coefficients0[3] = \
{0.471404520791032, -0.288675134594813, -0.166666666666667};
// Compute value(s)
for (unsigned int r = 0; r < 3; r++)
{
*values += coefficients0[r]*basisvalues[r];
} // end loop over 'r'
break;
}
case 1:
{
// Array of basisvalues
double basisvalues[3] = {0.0, 0.0, 0.0};
// Declare helper variables
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
// Compute basisvalues
basisvalues[0] = 1.0;
basisvalues[1] = tmp0;
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
basisvalues[0] *= std::sqrt(0.5);
basisvalues[2] *= std::sqrt(1.0);
basisvalues[1] *= std::sqrt(3.0);
// Table(s) of coefficients
static const double coefficients0[3] = \
{0.471404520791032, 0.288675134594813, -0.166666666666667};
// Compute value(s)
for (unsigned int r = 0; r < 3; r++)
{
*values += coefficients0[r]*basisvalues[r];
} // end loop over 'r'
break;
}
case 2:
{
// Array of basisvalues
double basisvalues[3] = {0.0, 0.0, 0.0};
// Declare helper variables
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
// Compute basisvalues
basisvalues[0] = 1.0;
basisvalues[1] = tmp0;
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
basisvalues[0] *= std::sqrt(0.5);
basisvalues[2] *= std::sqrt(1.0);
basisvalues[1] *= std::sqrt(3.0);
// Table(s) of coefficients
static const double coefficients0[3] = \
{0.471404520791032, 0.0, 0.333333333333333};
// Compute value(s)
for (unsigned int r = 0; r < 3; r++)
{
*values += coefficients0[r]*basisvalues[r];
} // end loop over 'r'
break;
}
}
}
/// Evaluate basis function i at given point x in cell (non-static member function)
virtual void evaluate_basis(std::size_t i,
double* values,
const double* x,
const double* vertex_coordinates,
int cell_orientation) const
{
_evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
}
/// Evaluate all basis functions at given point x in cell (actual implementation)
static void _evaluate_basis_all(double* values,
const double* x,
const double* vertex_coordinates,
int cell_orientation)
{
// Helper variable to hold values of a single dof.
double dof_values = 0.0;
// Loop dofs and call evaluate_basis
for (unsigned int r = 0; r < 3; r++)
{
_evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
values[r] = dof_values;
} // end loop over 'r'
}
/// Evaluate all basis functions at given point x in cell (non-static member function)
virtual void evaluate_basis_all(double* values,
const double* x,
const double* vertex_coordinates,
int cell_orientation) const
{
_evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
}
/// Evaluate order n derivatives of basis function i at given point x in cell (actual implementation)
static void _evaluate_basis_derivatives(std::size_t i,
std::size_t n,
double* values,
const double* x,
const double* vertex_coordinates,
int cell_orientation)
{
// Compute number of derivatives.
unsigned int num_derivatives = 1;
for (unsigned int r = 0; r < n; r++)
{
num_derivatives *= 2;
} // end loop over 'r'
// Reset values. Assuming that values is always an array.
for (unsigned int r = 0; r < num_derivatives; r++)
{
values[r] = 0.0;
} // end loop over 'r'
// Call evaluate_basis if order of derivatives is equal to zero.
if (n == 0)
{
_evaluate_basis(i, values, x, vertex_coordinates, cell_orientation);
return ;
}
// If order of derivatives is greater than the maximum polynomial degree, return zeros.
if (n > 1)
{
return ;
}
// Compute Jacobian
double J[4];
compute_jacobian_triangle_2d(J, vertex_coordinates);
// Compute Jacobian inverse and determinant
double K[4];
double detJ;
compute_jacobian_inverse_triangle_2d(K, detJ, J);
// Compute constants
const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
// Get coordinates and map to the reference (FIAT) element
double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
// Declare two dimensional array that holds combinations of derivatives and initialise
unsigned int combinations[2][1];
for (unsigned int row = 0; row < 2; row++)
{
for (unsigned int col = 0; col < 1; col++)
combinations[row][col] = 0;
}
// Generate combinations of derivatives
for (unsigned int row = 1; row < num_derivatives; row++)
{
for (unsigned int num = 0; num < row; num++)
{
for (unsigned int col = n-1; col+1 > 0; col--)
{
if (combinations[row][col] + 1 > 1)
combinations[row][col] = 0;
else
{
combinations[row][col] += 1;
break;
}
}
}
}
// Compute inverse of Jacobian
const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
// Declare transformation matrix
// Declare pointer to two dimensional array and initialise
double transform[2][2];
for (unsigned int j = 0; j < num_derivatives; j++)
{
for (unsigned int k = 0; k < num_derivatives; k++)
transform[j][k] = 1;
}
// Construct transformation matrix
for (unsigned int row = 0; row < num_derivatives; row++)
{
for (unsigned int col = 0; col < num_derivatives; col++)
{
for (unsigned int k = 0; k < n; k++)
transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
}
}
switch (i)
{
case 0:
{
// Array of basisvalues
double basisvalues[3] = {0.0, 0.0, 0.0};
// Declare helper variables
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
// Compute basisvalues
basisvalues[0] = 1.0;
basisvalues[1] = tmp0;
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
basisvalues[0] *= std::sqrt(0.5);
basisvalues[2] *= std::sqrt(1.0);
basisvalues[1] *= std::sqrt(3.0);
// Table(s) of coefficients
static const double coefficients0[3] = \
{0.471404520791032, -0.288675134594813, -0.166666666666667};
// Tables of derivatives of the polynomial base (transpose).
static const double dmats0[3][3] = \
{{0.0, 0.0, 0.0},
{4.89897948556635, 0.0, 0.0},
{0.0, 0.0, 0.0}};
static const double dmats1[3][3] = \
{{0.0, 0.0, 0.0},
{2.44948974278318, 0.0, 0.0},
{4.24264068711928, 0.0, 0.0}};
// Compute reference derivatives.
// Declare array of derivatives on FIAT element.
double derivatives[2];
for (unsigned int r = 0; r < 2; r++)
{
derivatives[r] = 0.0;
} // end loop over 'r'
// Declare derivative matrix (of polynomial basis).
double dmats[3][3] = \
{{1.0, 0.0, 0.0},
{0.0, 1.0, 0.0},
{0.0, 0.0, 1.0}};
// Declare (auxiliary) derivative matrix (of polynomial basis).
double dmats_old[3][3] = \
{{1.0, 0.0, 0.0},
{0.0, 1.0, 0.0},
{0.0, 0.0, 1.0}};
// Loop possible derivatives.
for (unsigned int r = 0; r < num_derivatives; r++)
{
// Resetting dmats values to compute next derivative.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats[t][u] = 0.0;
if (t == u)
{
dmats[t][u] = 1.0;
}
} // end loop over 'u'
} // end loop over 't'
// Looping derivative order to generate dmats.
for (unsigned int s = 0; s < n; s++)
{
// Updating dmats_old with new values and resetting dmats.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats_old[t][u] = dmats[t][u];
dmats[t][u] = 0.0;
} // end loop over 'u'
} // end loop over 't'
// Update dmats using an inner product.
if (combinations[r][s] == 0)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
} // end loop over 'tu'
} // end loop over 'u'
} // end loop over 't'
}
if (combinations[r][s] == 1)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
} // end loop over 'tu'
} // end loop over 'u'
} // end loop over 't'
}
} // end loop over 's'
for (unsigned int s = 0; s < 3; s++)
{
for (unsigned int t = 0; t < 3; t++)
{
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
} // end loop over 't'
} // end loop over 's'
} // end loop over 'r'
// Transform derivatives back to physical element
for (unsigned int r = 0; r < num_derivatives; r++)
{
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[r] += transform[r][s]*derivatives[s];
} // end loop over 's'
} // end loop over 'r'
break;
}
case 1:
{
// Array of basisvalues
double basisvalues[3] = {0.0, 0.0, 0.0};
// Declare helper variables
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
// Compute basisvalues
basisvalues[0] = 1.0;
basisvalues[1] = tmp0;
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
basisvalues[0] *= std::sqrt(0.5);
basisvalues[2] *= std::sqrt(1.0);
basisvalues[1] *= std::sqrt(3.0);
// Table(s) of coefficients
static const double coefficients0[3] = \
{0.471404520791032, 0.288675134594813, -0.166666666666667};
// Tables of derivatives of the polynomial base (transpose).
static const double dmats0[3][3] = \
{{0.0, 0.0, 0.0},
{4.89897948556635, 0.0, 0.0},
{0.0, 0.0, 0.0}};
static const double dmats1[3][3] = \
{{0.0, 0.0, 0.0},
{2.44948974278318, 0.0, 0.0},
{4.24264068711928, 0.0, 0.0}};
// Compute reference derivatives.
// Declare array of derivatives on FIAT element.
double derivatives[2];
for (unsigned int r = 0; r < 2; r++)
{
derivatives[r] = 0.0;
} // end loop over 'r'
// Declare derivative matrix (of polynomial basis).
double dmats[3][3] = \
{{1.0, 0.0, 0.0},
{0.0, 1.0, 0.0},
{0.0, 0.0, 1.0}};
// Declare (auxiliary) derivative matrix (of polynomial basis).
double dmats_old[3][3] = \
{{1.0, 0.0, 0.0},
{0.0, 1.0, 0.0},
{0.0, 0.0, 1.0}};
// Loop possible derivatives.
for (unsigned int r = 0; r < num_derivatives; r++)
{
// Resetting dmats values to compute next derivative.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats[t][u] = 0.0;
if (t == u)
{
dmats[t][u] = 1.0;
}
} // end loop over 'u'
} // end loop over 't'
// Looping derivative order to generate dmats.
for (unsigned int s = 0; s < n; s++)
{
// Updating dmats_old with new values and resetting dmats.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats_old[t][u] = dmats[t][u];
dmats[t][u] = 0.0;
} // end loop over 'u'
} // end loop over 't'
// Update dmats using an inner product.
if (combinations[r][s] == 0)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
} // end loop over 'tu'
} // end loop over 'u'
} // end loop over 't'
}
if (combinations[r][s] == 1)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
} // end loop over 'tu'
} // end loop over 'u'
} // end loop over 't'
}
} // end loop over 's'
for (unsigned int s = 0; s < 3; s++)
{
for (unsigned int t = 0; t < 3; t++)
{
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
} // end loop over 't'
} // end loop over 's'
} // end loop over 'r'
// Transform derivatives back to physical element
for (unsigned int r = 0; r < num_derivatives; r++)
{
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[r] += transform[r][s]*derivatives[s];
} // end loop over 's'
} // end loop over 'r'
break;
}
case 2:
{
// Array of basisvalues
double basisvalues[3] = {0.0, 0.0, 0.0};
// Declare helper variables
double tmp0 = (1.0 + Y + 2.0*X)/2.0;
// Compute basisvalues
basisvalues[0] = 1.0;
basisvalues[1] = tmp0;
basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
basisvalues[0] *= std::sqrt(0.5);
basisvalues[2] *= std::sqrt(1.0);
basisvalues[1] *= std::sqrt(3.0);
// Table(s) of coefficients
static const double coefficients0[3] = \
{0.471404520791032, 0.0, 0.333333333333333};
// Tables of derivatives of the polynomial base (transpose).
static const double dmats0[3][3] = \
{{0.0, 0.0, 0.0},
{4.89897948556635, 0.0, 0.0},
{0.0, 0.0, 0.0}};
static const double dmats1[3][3] = \
{{0.0, 0.0, 0.0},
{2.44948974278318, 0.0, 0.0},
{4.24264068711928, 0.0, 0.0}};
// Compute reference derivatives.
// Declare array of derivatives on FIAT element.
double derivatives[2];
for (unsigned int r = 0; r < 2; r++)
{
derivatives[r] = 0.0;
} // end loop over 'r'
// Declare derivative matrix (of polynomial basis).
double dmats[3][3] = \
{{1.0, 0.0, 0.0},
{0.0, 1.0, 0.0},
{0.0, 0.0, 1.0}};
// Declare (auxiliary) derivative matrix (of polynomial basis).
double dmats_old[3][3] = \
{{1.0, 0.0, 0.0},
{0.0, 1.0, 0.0},
{0.0, 0.0, 1.0}};
// Loop possible derivatives.
for (unsigned int r = 0; r < num_derivatives; r++)
{
// Resetting dmats values to compute next derivative.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats[t][u] = 0.0;
if (t == u)
{
dmats[t][u] = 1.0;
}
} // end loop over 'u'
} // end loop over 't'
// Looping derivative order to generate dmats.
for (unsigned int s = 0; s < n; s++)
{
// Updating dmats_old with new values and resetting dmats.
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
dmats_old[t][u] = dmats[t][u];
dmats[t][u] = 0.0;
} // end loop over 'u'
} // end loop over 't'
// Update dmats using an inner product.
if (combinations[r][s] == 0)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
} // end loop over 'tu'
} // end loop over 'u'
} // end loop over 't'
}
if (combinations[r][s] == 1)
{
for (unsigned int t = 0; t < 3; t++)
{
for (unsigned int u = 0; u < 3; u++)
{
for (unsigned int tu = 0; tu < 3; tu++)
{
dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
} // end loop over 'tu'
} // end loop over 'u'
} // end loop over 't'
}
} // end loop over 's'
for (unsigned int s = 0; s < 3; s++)
{
for (unsigned int t = 0; t < 3; t++)
{
derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
} // end loop over 't'
} // end loop over 's'
} // end loop over 'r'
// Transform derivatives back to physical element
for (unsigned int r = 0; r < num_derivatives; r++)
{
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[r] += transform[r][s]*derivatives[s];
} // end loop over 's'
} // end loop over 'r'
break;
}
}
}
/// Evaluate order n derivatives of basis function i at given point x in cell (non-static member function)
virtual void evaluate_basis_derivatives(std::size_t i,
std::size_t n,
double* values,
const double* x,
const double* vertex_coordinates,
int cell_orientation) const
{
_evaluate_basis_derivatives(i, n, values, x, vertex_coordinates, cell_orientation);
}
/// Evaluate order n derivatives of all basis functions at given point x in cell (actual implementation)
static void _evaluate_basis_derivatives_all(std::size_t n,
double* values,
const double* x,
const double* vertex_coordinates,
int cell_orientation)
{
// Call evaluate_basis_all if order of derivatives is equal to zero.
if (n == 0)
{
_evaluate_basis_all(values, x, vertex_coordinates, cell_orientation);
return ;
}
// Compute number of derivatives.
unsigned int num_derivatives = 1;
for (unsigned int r = 0; r < n; r++)
{
num_derivatives *= 2;
} // end loop over 'r'
// Set values equal to zero.
for (unsigned int r = 0; r < 3; r++)
{
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[r*num_derivatives + s] = 0.0;
} // end loop over 's'
} // end loop over 'r'
// If order of derivatives is greater than the maximum polynomial degree, return zeros.
if (n > 1)
{
return ;
}
// Helper variable to hold values of a single dof.
double dof_values[2];
for (unsigned int r = 0; r < 2; r++)
{
dof_values[r] = 0.0;
} // end loop over 'r'
// Loop dofs and call evaluate_basis_derivatives.
for (unsigned int r = 0; r < 3; r++)
{
_evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
for (unsigned int s = 0; s < num_derivatives; s++)
{
values[r*num_derivatives + s] = dof_values[s];
} // end loop over 's'
} // end loop over 'r'
}
/// Evaluate order n derivatives of all basis functions at given point x in cell (non-static member function)
virtual void evaluate_basis_derivatives_all(std::size_t n,
double* values,
const double* x,
const double* vertex_coordinates,
int cell_orientation) const
{
_evaluate_basis_derivatives_all(n, values, x, vertex_coordinates, cell_orientation);
}
/// Evaluate linear functional for dof i on the function f
virtual double evaluate_dof(std::size_t i,
const ufc::function& f,
const double* vertex_coordinates,
int cell_orientation,
const ufc::cell& c) const
{
// Declare variables for result of evaluation
double vals[1];
// Declare variable for physical coordinates
double y[2];
switch (i)
{
case 0:
{
y[0] = vertex_coordinates[0];
y[1] = vertex_coordinates[1];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 1:
{
y[0] = vertex_coordinates[2];
y[1] = vertex_coordinates[3];
f.evaluate(vals, y, c);
return vals[0];
break;
}
case 2:
{
y[0] = vertex_coordinates[4];
y[1] = vertex_coordinates[5];
f.evaluate(vals, y, c);
return vals[0];
break;
}
}
return 0.0;
}
/// Evaluate linear functionals for all dofs on the function f
virtual void evaluate_dofs(double* values,
const ufc::function& f,
const double* vertex_coordinates,
int cell_orientation,
const ufc::cell& c) const
{
// Declare variables for result of evaluation
double vals[1];
// Declare variable for physical coordinates
double y[2];
y[0] = vertex_coordinates[0];
y[1] = vertex_coordinates[1];
f.evaluate(vals, y, c);
values[0] = vals[0];
y[0] = vertex_coordinates[2];
y[1] = vertex_coordinates[3];
f.evaluate(vals, y, c);
values[1] = vals[0];
y[0] = vertex_coordinates[4];
y[1] = vertex_coordinates[5];
f.evaluate(vals, y, c);
values[2] = vals[0];
}
/// Interpolate vertex values from dof values
virtual void interpolate_vertex_values(double* vertex_values,
const double* dof_values,
const double* vertex_coordinates,
int cell_orientation,
const ufc::cell& c) const
{
// Evaluate function and change variables
vertex_values[0] = dof_values[0];
vertex_values[1] = dof_values[1];
vertex_values[2] = dof_values[2];
}
/// Map coordinate xhat from reference cell to coordinate x in cell
virtual void map_from_reference_cell(double* x,
const double* xhat,
const ufc::cell& c) const
{
throw std::runtime_error("map_from_reference_cell not yet implemented.");
}
/// Map from coordinate x in cell to coordinate xhat in reference cell
virtual void map_to_reference_cell(double* xhat,
const double* x,
const ufc::cell& c) const
{
throw std::runtime_error("map_to_reference_cell not yet implemented.");
}
/// Return the number of sub elements (for a mixed element)
virtual std::size_t num_sub_elements() const
{
return 0;
}
/// Create a new finite element for sub element i (for a mixed element)
virtual ufc::finite_element* create_sub_element(std::size_t i) const
{
return 0;
}
/// Create a new class instance
virtual ufc::finite_element* create() const
{
return new poisson_finite_element_0();
}
};
/// This class defines the interface for a local-to-global mapping of
/// degrees of freedom (dofs).
class poisson_dofmap_0: public ufc::dofmap
{
public:
/// Constructor
poisson_dofmap_0() : ufc::dofmap()
{
// Do nothing
}
/// Destructor
virtual ~poisson_dofmap_0()
{
// Do nothing
}
/// Return a string identifying the dofmap
virtual const char* signature() const
{
return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2)), 1, None)";
}
/// Return true iff mesh entities of topological dimension d are needed
virtual bool needs_mesh_entities(std::size_t d) const
{
switch (d)
{
case 0:
{
return true;
break;
}
case 1:
{
return false;
break;
}
case 2:
{
return false;
break;
}
}
return false;
}
/// Return the topological dimension of the associated cell shape
virtual std::size_t topological_dimension() const
{
return 2;
}
/// Return the geometric dimension of the associated cell shape
virtual std::size_t geometric_dimension() const
{
return 2;
}
/// Return the dimension of the global finite element function space
virtual std::size_t global_dimension(const std::vector<std::size_t>&
num_global_entities) const
{
return num_global_entities[0];
}
/// Return the dimension of the local finite element function space for a cell
virtual std::size_t num_element_dofs() const
{
return 3;
}
/// Return the number of dofs on each cell facet
virtual std::size_t num_facet_dofs() const
{
return 2;
}
/// Return the number of dofs associated with each cell entity of dimension d
virtual std::size_t num_entity_dofs(std::size_t d) const
{
switch (d)
{
case 0:
{
return 1;
break;
}
case 1:
{
return 0;
break;
}
case 2:
{
return 0;
break;
}
}
return 0;
}
/// Tabulate the local-to-global mapping of dofs on a cell
virtual void tabulate_dofs(std::size_t* dofs,
const std::vector<std::size_t>& num_global_entities,
const ufc::cell& c) const
{
dofs[0] = c.entity_indices[0][0];
dofs[1] = c.entity_indices[0][1];
dofs[2] = c.entity_indices[0][2];
}
/// Tabulate the local-to-local mapping from facet dofs to cell dofs
virtual void tabulate_facet_dofs(std::size_t* dofs,
std::size_t facet) const
{
switch (facet)
{
case 0:
{
dofs[0] = 1;
dofs[1] = 2;
break;
}
case 1:
{
dofs[0] = 0;
dofs[1] = 2;
break;
}
case 2:
{
dofs[0] = 0;
dofs[1] = 1;
break;
}
}
}
/// Tabulate the local-to-local mapping of dofs on entity (d, i)
virtual void tabulate_entity_dofs(std::size_t* dofs,
std::size_t d, std::size_t i) const
{
if (d > 2)
{
throw std::runtime_error("d is larger than dimension (2)");
}
switch (d)
{
case 0:
{
if (i > 2)
{
throw std::runtime_error("i is larger than number of entities (2)");
}
switch (i)
{
case 0:
{
dofs[0] = 0;
break;
}
case 1:
{
dofs[0] = 1;
break;
}
case 2:
{
dofs[0] = 2;
break;
}
}
break;
}
case 1:
{
break;
}
case 2:
{
break;
}
}
}
/// Tabulate the coordinates of all dofs on a cell
virtual void tabulate_coordinates(double* dof_coordinates,
const double* vertex_coordinates) const
{
dof_coordinates[0] = vertex_coordinates[0];
dof_coordinates[1] = vertex_coordinates[1];
dof_coordinates[2] = vertex_coordinates[2];
dof_coordinates[3] = vertex_coordinates[3];
dof_coordinates[4] = vertex_coordinates[4];
dof_coordinates[5] = vertex_coordinates[5];
}
/// Return the number of sub dofmaps (for a mixed element)
virtual std::size_t num_sub_dofmaps() const
{
return 0;
}
/// Create a new dofmap for sub dofmap i (for a mixed element)
virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
{
return 0;
}
/// Create a new class instance
virtual ufc::dofmap* create() const
{
return new poisson_dofmap_0();
}
};
/// This class defines the interface for the tabulation of the cell
/// tensor corresponding to the local contribution to a form from
/// the integral over a cell.
class poisson_cell_integral_0_otherwise: public ufc::cell_integral
{
public:
/// Constructor
poisson_cell_integral_0_otherwise() : ufc::cell_integral()
{
// Do nothing
}
/// Destructor
virtual ~poisson_cell_integral_0_otherwise()
{
// Do nothing
}
/// Tabulate which form coefficients are used by this integral
virtual const std::vector<bool> & enabled_coefficients() const
{
static const std::vector<bool> enabled({});
return enabled;
}
/// Tabulate the tensor for the contribution from a local cell
virtual void tabulate_tensor(double* A,
const double * const * w,
const double* vertex_coordinates,
int cell_orientation) const
{
// Number of operations (multiply-add pairs) for Jacobian data: 3
// Number of operations (multiply-add pairs) for geometry tensor: 8
// Number of operations (multiply-add pairs) for tensor contraction: 11
// Total number of operations (multiply-add pairs): 22
// Compute Jacobian
double J[4];
compute_jacobian_triangle_2d(J, vertex_coordinates);
// Compute Jacobian inverse and determinant
double K[4];
double detJ;
compute_jacobian_inverse_triangle_2d(K, detJ, J);
// Set scale factor
const double det = std::abs(detJ);
// Compute geometry tensor
const double G0_0_0 = det*(K[0]*K[0] + K[1]*K[1]);
const double G0_0_1 = det*(K[0]*K[2] + K[1]*K[3]);
const double G0_1_0 = det*(K[2]*K[0] + K[3]*K[1]);
const double G0_1_1 = det*(K[2]*K[2] + K[3]*K[3]);
// Compute element tensor
A[0] = 0.499999999999999*G0_0_0 + 0.5*G0_0_1 + 0.5*G0_1_0 + 0.5*G0_1_1;
A[1] = -0.499999999999999*G0_0_0 - 0.5*G0_1_0;
A[2] = -0.5*G0_0_1 - 0.5*G0_1_1;
A[3] = -0.499999999999999*G0_0_0 - 0.5*G0_0_1;
A[4] = 0.499999999999999*G0_0_0;
A[5] = 0.5*G0_0_1;
A[6] = -0.5*G0_1_0 - 0.5*G0_1_1;
A[7] = 0.5*G0_1_0;
A[8] = 0.5*G0_1_1;
}
};
/// This class defines the interface for the tabulation of the cell
/// tensor corresponding to the local contribution to a form from
/// the integral over a cell.
class poisson_cell_integral_1_otherwise: public ufc::cell_integral
{
public:
/// Constructor
poisson_cell_integral_1_otherwise() : ufc::cell_integral()
{
// Do nothing
}
/// Destructor
virtual ~poisson_cell_integral_1_otherwise()
{
// Do nothing
}
/// Tabulate which form coefficients are used by this integral
virtual const std::vector<bool> & enabled_coefficients() const
{
static const std::vector<bool> enabled({true, false});
return enabled;
}
/// Tabulate the tensor for the contribution from a local cell
virtual void tabulate_tensor(double* A,
const double * const * w,
const double* vertex_coordinates,
int cell_orientation) const
{
// Number of operations (multiply-add pairs) for Jacobian data: 3
// Number of operations (multiply-add pairs) for geometry tensor: 3
// Number of operations (multiply-add pairs) for tensor contraction: 7
// Total number of operations (multiply-add pairs): 13
// Compute Jacobian
double J[4];
compute_jacobian_triangle_2d(J, vertex_coordinates);
// Compute Jacobian inverse and determinant
double K[4];
double detJ;
compute_jacobian_inverse_triangle_2d(K, detJ, J);
// Set scale factor
const double det = std::abs(detJ);
// Compute geometry tensor
const double G0_0 = det*w[0][0]*(1.0);
const double G0_1 = det*w[0][1]*(1.0);
const double G0_2 = det*w[0][2]*(1.0);
// Compute element tensor
A[0] = 0.0833333333333334*G0_0 + 0.0416666666666667*G0_1 + 0.0416666666666667*G0_2;
A[1] = 0.0416666666666667*G0_0 + 0.0833333333333333*G0_1 + 0.0416666666666666*G0_2;
A[2] = 0.0416666666666667*G0_0 + 0.0416666666666666*G0_1 + 0.0833333333333333*G0_2;
}
};
/// This class defines the interface for the tabulation of the
/// exterior facet tensor corresponding to the local contribution to
/// a form from the integral over an exterior facet.
class poisson_exterior_facet_integral_1_otherwise: public ufc::exterior_facet_integral
{
public:
/// Constructor
poisson_exterior_facet_integral_1_otherwise() : ufc::exterior_facet_integral()
{
// Do nothing
}
/// Destructor
virtual ~poisson_exterior_facet_integral_1_otherwise()
{
// Do nothing
}
/// Tabulate which form coefficients are used by this integral
virtual const std::vector<bool> & enabled_coefficients() const
{
static const std::vector<bool> enabled({false, true});
return enabled;
}
/// Tabulate the tensor for the contribution from a local exterior facet
virtual void tabulate_tensor(double* A,
const double * const * w,
const double* vertex_coordinates,
std::size_t facet,
int cell_orientation) const
{
// Number of operations (multiply-add pairs) for Jacobian data: 10
// Number of operations (multiply-add pairs) for geometry tensor: 3
// Number of operations (multiply-add pairs) for tensor contraction: 9
// Total number of operations (multiply-add pairs): 22
// Compute Jacobian
double J[4];
compute_jacobian_triangle_2d(J, vertex_coordinates);
// Compute Jacobian inverse and determinant
double K[4];
double detJ;
compute_jacobian_inverse_triangle_2d(K, detJ, J);
// Get vertices on edge
static unsigned int edge_vertices[3][2] = {{1, 2}, {0, 2}, {0, 1}};
const unsigned int v0 = edge_vertices[facet][0];
const unsigned int v1 = edge_vertices[facet][1];
// Compute scale factor (length of edge scaled by length of reference interval)
const double dx0 = vertex_coordinates[2*v1 + 0] - vertex_coordinates[2*v0 + 0];
const double dx1 = vertex_coordinates[2*v1 + 1] - vertex_coordinates[2*v0 + 1];
const double det = std::sqrt(dx0*dx0 + dx1*dx1);
// Compute geometry tensor
const double G0_0 = det*w[1][0]*(1.0);
const double G0_1 = det*w[1][1]*(1.0);
const double G0_2 = det*w[1][2]*(1.0);
// Compute element tensor
switch (facet)
{
case 0:
{
A[0] = 0.0;
A[1] = 0.333333333333333*G0_1 + 0.166666666666667*G0_2;
A[2] = 0.166666666666667*G0_1 + 0.333333333333333*G0_2;
break;
}
case 1:
{
A[0] = 0.333333333333333*G0_0 + 0.166666666666667*G0_2;
A[1] = 0.0;
A[2] = 0.166666666666667*G0_0 + 0.333333333333333*G0_2;
break;
}
case 2:
{
A[0] = 0.333333333333333*G0_0 + 0.166666666666667*G0_1;
A[1] = 0.166666666666667*G0_0 + 0.333333333333333*G0_1;
A[2] = 0.0;
break;
}
}
}
};
/// This class defines the interface for the assembly of the global
/// tensor corresponding to a form with r + n arguments, that is, a
/// mapping
///
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
///
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
/// global tensor A is defined by
///
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
///
/// where each argument Vj represents the application to the
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
/// fixed functions (coefficients).
class poisson_form_0: public ufc::form
{
public:
/// Constructor
poisson_form_0() : ufc::form()
{
// Do nothing
}
/// Destructor
virtual ~poisson_form_0()
{
// Do nothing
}
/// Return a string identifying the form
virtual const char* signature() const
{
return "7f3815c8e31ab18b282532dab8742e6961eab53ad4e81fa7cb09d956cec41e4af64946224d6780069be4537652b627b48dcbe12cb4f41895baf526eb4e5dcc56";
}
/// Return the rank of the global tensor (r)
virtual std::size_t rank() const
{
return 2;
}
/// Return the number of coefficients (n)
virtual std::size_t num_coefficients() const
{
return 0;
}
/// Return original coefficient position for each coefficient (0 <= i < n)
virtual std::size_t original_coefficient_position(std::size_t i) const
{
static const std::vector<std::size_t> position({});
return position[i];
}
/// Create a new finite element for argument function i
virtual ufc::finite_element* create_finite_element(std::size_t i) const
{
switch (i)
{
case 0:
{
return new poisson_finite_element_0();
break;
}
case 1:
{
return new poisson_finite_element_0();
break;
}
}
return 0;
}
/// Create a new dofmap for argument function i
virtual ufc::dofmap* create_dofmap(std::size_t i) const
{
switch (i)
{
case 0:
{
return new poisson_dofmap_0();
break;
}
case 1:
{
return new poisson_dofmap_0();
break;
}
}
return 0;
}
/// Return the number of cell domains
virtual std::size_t max_cell_subdomain_id() const
{
return 0;
}
/// Return the number of exterior facet domains
virtual std::size_t max_exterior_facet_subdomain_id() const
{
return 0;
}
/// Return the number of interior facet domains
virtual std::size_t max_interior_facet_subdomain_id() const
{
return 0;
}
/// Return the number of vertex domains
virtual std::size_t max_vertex_subdomain_id() const
{
return 0;
}
/// Return the number of custom domains
virtual std::size_t max_custom_subdomain_id() const
{
return 0;
}
/// Return whether the form has any cell integrals
virtual bool has_cell_integrals() const
{
return true;
}
/// Return whether the form has any exterior facet integrals
virtual bool has_exterior_facet_integrals() const
{
return false;
}
/// Return whether the form has any interior facet integrals
virtual bool has_interior_facet_integrals() const
{
return false;
}
/// Return whether the form has any vertex integrals
virtual bool has_vertex_integrals() const
{
return false;
}
/// Return whether the form has any custom integrals
virtual bool has_custom_integrals() const
{
return false;
}
/// Create a new cell integral on sub domain subdomain_id
virtual ufc::cell_integral* create_cell_integral(std::size_t subdomain_id) const
{
return 0;
}
/// Create a new exterior facet integral on sub domain subdomain_id
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t subdomain_id) const
{
return 0;
}
/// Create a new interior facet integral on sub domain subdomain_id
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t subdomain_id) const
{
return 0;
}
/// Create a new vertex integral on sub domain subdomain_id
virtual ufc::vertex_integral* create_vertex_integral(std::size_t subdomain_id) const
{
return 0;
}
/// Create a new custom integral on sub domain subdomain_id
virtual ufc::custom_integral* create_custom_integral(std::size_t subdomain_id) const
{
return 0;
}
/// Create a new cell integral on everywhere else
virtual ufc::cell_integral* create_default_cell_integral() const
{
return new poisson_cell_integral_0_otherwise();
}
/// Create a new exterior facet integral on everywhere else
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
{
return 0;
}
/// Create a new interior facet integral on everywhere else
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
{
return 0;
}
/// Create a new vertex integral on everywhere else
virtual ufc::vertex_integral* create_default_vertex_integral() const
{
return 0;
}
/// Create a new custom integral on everywhere else
virtual ufc::custom_integral* create_default_custom_integral() const
{
return 0;
}
};
/// This class defines the interface for the assembly of the global
/// tensor corresponding to a form with r + n arguments, that is, a
/// mapping
///
/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
///
/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
/// global tensor A is defined by
///
/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
///
/// where each argument Vj represents the application to the
/// sequence of basis functions of Vj and w1, w2, ..., wn are given
/// fixed functions (coefficients).
class poisson_form_1: public ufc::form
{
public:
/// Constructor
poisson_form_1() : ufc::form()
{
// Do nothing
}
/// Destructor
virtual ~poisson_form_1()
{
// Do nothing
}
/// Return a string identifying the form
virtual const char* signature() const
{
return "e58507906d4ed11ee413a96a0fc6c639025a4127e4d90f9e99c7bee85f398d70855d6f5e3ae70dcc968b6b0f7967bb5fa54f876c26904fcb75b8950e4bc622c4";
}
/// Return the rank of the global tensor (r)
virtual std::size_t rank() const
{
return 1;
}
/// Return the number of coefficients (n)
virtual std::size_t num_coefficients() const
{
return 2;
}
/// Return original coefficient position for each coefficient (0 <= i < n)
virtual std::size_t original_coefficient_position(std::size_t i) const
{
static const std::vector<std::size_t> position({0, 1});
return position[i];
}
/// Create a new finite element for argument function i
virtual ufc::finite_element* create_finite_element(std::size_t i) const
{
switch (i)
{
case 0:
{
return new poisson_finite_element_0();
break;
}
case 1:
{
return new poisson_finite_element_0();
break;
}
case 2:
{
return new poisson_finite_element_0();
break;
}
}
return 0;
}
/// Create a new dofmap for argument function i
virtual ufc::dofmap* create_dofmap(std::size_t i) const
{
switch (i)
{
case 0:
{
return new poisson_dofmap_0();
break;
}
case 1:
{
return new poisson_dofmap_0();
break;
}
case 2:
{
return new poisson_dofmap_0();
break;
}
}
return 0;
}
/// Return the number of cell domains
virtual std::size_t max_cell_subdomain_id() const
{
return 0;
}
/// Return the number of exterior facet domains
virtual std::size_t max_exterior_facet_subdomain_id() const
{
return 0;
}
/// Return the number of interior facet domains
virtual std::size_t max_interior_facet_subdomain_id() const
{
return 0;
}
/// Return the number of vertex domains
virtual std::size_t max_vertex_subdomain_id() const
{
return 0;
}
/// Return the number of custom domains
virtual std::size_t max_custom_subdomain_id() const
{
return 0;
}
/// Return whether the form has any cell integrals
virtual bool has_cell_integrals() const
{
return true;
}
/// Return whether the form has any exterior facet integrals
virtual bool has_exterior_facet_integrals() const
{
return true;
}
/// Return whether the form has any interior facet integrals
virtual bool has_interior_facet_integrals() const
{
return false;
}
/// Return whether the form has any vertex integrals
virtual bool has_vertex_integrals() const
{
return false;
}
/// Return whether the form has any custom integrals
virtual bool has_custom_integrals() const
{
return false;
}
/// Create a new cell integral on sub domain subdomain_id
virtual ufc::cell_integral* create_cell_integral(std::size_t subdomain_id) const
{
return 0;
}
/// Create a new exterior facet integral on sub domain subdomain_id
virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t subdomain_id) const
{
return 0;
}
/// Create a new interior facet integral on sub domain subdomain_id
virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t subdomain_id) const
{
return 0;
}
/// Create a new vertex integral on sub domain subdomain_id
virtual ufc::vertex_integral* create_vertex_integral(std::size_t subdomain_id) const
{
return 0;
}
/// Create a new custom integral on sub domain subdomain_id
virtual ufc::custom_integral* create_custom_integral(std::size_t subdomain_id) const
{
return 0;
}
/// Create a new cell integral on everywhere else
virtual ufc::cell_integral* create_default_cell_integral() const
{
return new poisson_cell_integral_1_otherwise();
}
/// Create a new exterior facet integral on everywhere else
virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
{
return new poisson_exterior_facet_integral_1_otherwise();
}
/// Create a new interior facet integral on everywhere else
virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
{
return 0;
}
/// Create a new vertex integral on everywhere else
virtual ufc::vertex_integral* create_default_vertex_integral() const
{
return 0;
}
/// Create a new custom integral on everywhere else
virtual ufc::custom_integral* create_default_custom_integral() const
{
return 0;
}
};
// DOLFIN wrappers
// Standard library includes
#include <string>
// DOLFIN includes
#include <dolfin/common/NoDeleter.h>
#include <dolfin/fem/FiniteElement.h>
#include <dolfin/fem/DofMap.h>
#include <dolfin/fem/Form.h>
#include <dolfin/function/FunctionSpace.h>
#include <dolfin/function/GenericFunction.h>
#include <dolfin/function/CoefficientAssigner.h>
#include <dolfin/adaptivity/ErrorControl.h>
#include <dolfin/adaptivity/GoalFunctional.h>
namespace Poisson
{
class CoefficientSpace_f: public dolfin::FunctionSpace
{
public:
//--- Constructors for standard function space, 2 different versions ---
// Create standard function space (reference version)
CoefficientSpace_f(const dolfin::Mesh& mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), mesh)))
{
// Do nothing
}
// Create standard function space (shared pointer version)
CoefficientSpace_f(std::shared_ptr<const dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), *mesh)))
{
// Do nothing
}
//--- Constructors for constrained function space, 2 different versions ---
// Create standard function space (reference version)
CoefficientSpace_f(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), mesh,
dolfin::reference_to_no_delete_pointer(constrained_domain))))
{
// Do nothing
}
// Create standard function space (shared pointer version)
CoefficientSpace_f(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
dolfin::FunctionSpace(mesh,
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), *mesh, constrained_domain)))
{
// Do nothing
}
};
class CoefficientSpace_g: public dolfin::FunctionSpace
{
public:
//--- Constructors for standard function space, 2 different versions ---
// Create standard function space (reference version)
CoefficientSpace_g(const dolfin::Mesh& mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), mesh)))
{
// Do nothing
}
// Create standard function space (shared pointer version)
CoefficientSpace_g(std::shared_ptr<const dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), *mesh)))
{
// Do nothing
}
//--- Constructors for constrained function space, 2 different versions ---
// Create standard function space (reference version)
CoefficientSpace_g(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), mesh,
dolfin::reference_to_no_delete_pointer(constrained_domain))))
{
// Do nothing
}
// Create standard function space (shared pointer version)
CoefficientSpace_g(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
dolfin::FunctionSpace(mesh,
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), *mesh, constrained_domain)))
{
// Do nothing
}
};
class Form_a_FunctionSpace_0: public dolfin::FunctionSpace
{
public:
//--- Constructors for standard function space, 2 different versions ---
// Create standard function space (reference version)
Form_a_FunctionSpace_0(const dolfin::Mesh& mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), mesh)))
{
// Do nothing
}
// Create standard function space (shared pointer version)
Form_a_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), *mesh)))
{
// Do nothing
}
//--- Constructors for constrained function space, 2 different versions ---
// Create standard function space (reference version)
Form_a_FunctionSpace_0(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), mesh,
dolfin::reference_to_no_delete_pointer(constrained_domain))))
{
// Do nothing
}
// Create standard function space (shared pointer version)
Form_a_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
dolfin::FunctionSpace(mesh,
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), *mesh, constrained_domain)))
{
// Do nothing
}
};
class Form_a_FunctionSpace_1: public dolfin::FunctionSpace
{
public:
//--- Constructors for standard function space, 2 different versions ---
// Create standard function space (reference version)
Form_a_FunctionSpace_1(const dolfin::Mesh& mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), mesh)))
{
// Do nothing
}
// Create standard function space (shared pointer version)
Form_a_FunctionSpace_1(std::shared_ptr<const dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), *mesh)))
{
// Do nothing
}
//--- Constructors for constrained function space, 2 different versions ---
// Create standard function space (reference version)
Form_a_FunctionSpace_1(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), mesh,
dolfin::reference_to_no_delete_pointer(constrained_domain))))
{
// Do nothing
}
// Create standard function space (shared pointer version)
Form_a_FunctionSpace_1(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
dolfin::FunctionSpace(mesh,
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), *mesh, constrained_domain)))
{
// Do nothing
}
};
class Form_a: public dolfin::Form
{
public:
// Constructor
Form_a(const dolfin::FunctionSpace& V1, const dolfin::FunctionSpace& V0):
dolfin::Form(2, 0)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
_function_spaces[1] = reference_to_no_delete_pointer(V1);
_ufc_form = std::shared_ptr<const ufc::form>(new poisson_form_0());
}
// Constructor
Form_a(std::shared_ptr<const dolfin::FunctionSpace> V1, std::shared_ptr<const dolfin::FunctionSpace> V0):
dolfin::Form(2, 0)
{
_function_spaces[0] = V0;
_function_spaces[1] = V1;
_ufc_form = std::shared_ptr<const ufc::form>(new poisson_form_0());
}
// Destructor
~Form_a()
{}
/// Return the number of the coefficient with this name
virtual std::size_t coefficient_number(const std::string& name) const
{
dolfin::dolfin_error("generated code for class Form",
"access coefficient data",
"There are no coefficients");
return 0;
}
/// Return the name of the coefficient with this number
virtual std::string coefficient_name(std::size_t i) const
{
dolfin::dolfin_error("generated code for class Form",
"access coefficient data",
"There are no coefficients");
return "unnamed";
}
// Typedefs
typedef Form_a_FunctionSpace_0 TestSpace;
typedef Form_a_FunctionSpace_1 TrialSpace;
// Coefficients
};
class Form_L_FunctionSpace_0: public dolfin::FunctionSpace
{
public:
//--- Constructors for standard function space, 2 different versions ---
// Create standard function space (reference version)
Form_L_FunctionSpace_0(const dolfin::Mesh& mesh):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), mesh)))
{
// Do nothing
}
// Create standard function space (shared pointer version)
Form_L_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh):
dolfin::FunctionSpace(mesh,
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), *mesh)))
{
// Do nothing
}
//--- Constructors for constrained function space, 2 different versions ---
// Create standard function space (reference version)
Form_L_FunctionSpace_0(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), mesh,
dolfin::reference_to_no_delete_pointer(constrained_domain))))
{
// Do nothing
}
// Create standard function space (shared pointer version)
Form_L_FunctionSpace_0(std::shared_ptr<const dolfin::Mesh> mesh, std::shared_ptr<const dolfin::SubDomain> constrained_domain):
dolfin::FunctionSpace(mesh,
std::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(std::shared_ptr<ufc::finite_element>(new poisson_finite_element_0()))),
std::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(std::shared_ptr<ufc::dofmap>(new poisson_dofmap_0()), *mesh, constrained_domain)))
{
// Do nothing
}
};
typedef CoefficientSpace_f Form_L_FunctionSpace_1;
typedef CoefficientSpace_g Form_L_FunctionSpace_2;
class Form_L: public dolfin::Form
{
public:
// Constructor
Form_L(const dolfin::FunctionSpace& V0):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
_ufc_form = std::shared_ptr<const ufc::form>(new poisson_form_1());
}
// Constructor
Form_L(const dolfin::FunctionSpace& V0, const dolfin::GenericFunction& f, const dolfin::GenericFunction& g):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
this->f = f;
this->g = g;
_ufc_form = std::shared_ptr<const ufc::form>(new poisson_form_1());
}
// Constructor
Form_L(const dolfin::FunctionSpace& V0, std::shared_ptr<const dolfin::GenericFunction> f, std::shared_ptr<const dolfin::GenericFunction> g):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = reference_to_no_delete_pointer(V0);
this->f = *f;
this->g = *g;
_ufc_form = std::shared_ptr<const ufc::form>(new poisson_form_1());
}
// Constructor
Form_L(std::shared_ptr<const dolfin::FunctionSpace> V0):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = V0;
_ufc_form = std::shared_ptr<const ufc::form>(new poisson_form_1());
}
// Constructor
Form_L(std::shared_ptr<const dolfin::FunctionSpace> V0, const dolfin::GenericFunction& f, const dolfin::GenericFunction& g):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = V0;
this->f = f;
this->g = g;
_ufc_form = std::shared_ptr<const ufc::form>(new poisson_form_1());
}
// Constructor
Form_L(std::shared_ptr<const dolfin::FunctionSpace> V0, std::shared_ptr<const dolfin::GenericFunction> f, std::shared_ptr<const dolfin::GenericFunction> g):
dolfin::Form(1, 2), f(*this, 0), g(*this, 1)
{
_function_spaces[0] = V0;
this->f = *f;
this->g = *g;
_ufc_form = std::shared_ptr<const ufc::form>(new poisson_form_1());
}
// Destructor
~Form_L()
{}
/// Return the number of the coefficient with this name
virtual std::size_t coefficient_number(const std::string& name) const
{
if (name == "f")
return 0;
else if (name == "g")
return 1;
dolfin::dolfin_error("generated code for class Form",
"access coefficient data",
"Invalid coefficient");
return 0;
}
/// Return the name of the coefficient with this number
virtual std::string coefficient_name(std::size_t i) const
{
switch (i)
{
case 0:
return "f";
case 1:
return "g";
}
dolfin::dolfin_error("generated code for class Form",
"access coefficient data",
"Invalid coefficient");
return "unnamed";
}
// Typedefs
typedef Form_L_FunctionSpace_0 TestSpace;
typedef Form_L_FunctionSpace_1 CoefficientSpace_f;
typedef Form_L_FunctionSpace_2 CoefficientSpace_g;
// Coefficients
dolfin::CoefficientAssigner f;
dolfin::CoefficientAssigner g;
};
// Class typedefs
typedef Form_a BilinearForm;
typedef Form_a JacobianForm;
typedef Form_L LinearForm;
typedef Form_L ResidualForm;
typedef Form_a::TestSpace FunctionSpace;
}
#endif
Raw
run.sh
export LD_PRELOAD=/usr/lib/openmpi/lib/libmpi.so:/usr/lib/openmpi/lib/libmpi_cxx.so:/usr/lib/liblapack.so:/usr/lib/libblas.so:/lib/x86_64-linux-gnu/libcrypto.so.1.0.0:/usr/lib/x86_64-linux-gnu/libcurl.so.4:/lib/x86_64-linux-gnu/libexpat.so.1:/usr/lib/x86_64-linux-gnu/libfreetype.so.6:/lib/x86_64-linux-gnu/libpng12.so.0:/usr/lib/x86_64-linux-gnu/libQtCore.so.4:/usr/lib/x86_64-linux-gnu/libQtGui.so.4:/usr/lib/x86_64-linux-gnu/libQtGui.so.4:/usr/lib/x86_64-linux-gnu/libQtSql.so.4:/lib/x86_64-linux-gnu/libssl.so.1.0.0:/usr/lib/x86_64-linux-gnu/libtiff.so.5 && matlab -nodesktop -nosplash
Raw
test.cpp
#include "mex.h"
#include <dolfin.h>
#include "Poisson.h"
using namespace dolfin;
// Source term (right-hand side)
class Source : public Expression
{
void eval(Array<double>& values, const Array<double>& x) const
{
double dx = x[0] - 0.5;
double dy = x[1] - 0.5;
values[0] = 10*exp(-(dx*dx + dy*dy) / 0.02);
}
};
// Normal derivative (Neumann boundary condition)
class dUdN : public Expression
{
void eval(Array<double>& values, const Array<double>& x) const
{
values[0] = sin(5*x[0]);
}
};
// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
bool inside(const Array<double>& x, bool on_boundary) const
{
return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS;
}
};
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
UnitSquareMesh mesh(320, 320);
Poisson::FunctionSpace V(mesh);
// Define boundary condition
Constant u0(0.0);
DirichletBoundary boundary;
DirichletBC bc(V, u0, boundary);
// Define variational forms
Poisson::BilinearForm a(V, V);
Poisson::LinearForm L(V);
Source f;
dUdN g;
L.f = f;
L.g = g;
// Compute solution
Function u(V);
solve(a == L, u, bc);
//Save solution in VTK format
// File file("poisson.pvd");
// file << u;
//
// plot(u);
// interactive();
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment