DOLFIN

Merge lp:~corrado-maurini/dolfin/tao into lp:~fenics-core/dolfin/trunk

  1. tao
  2. Merge into trunk
Proposed by corrado maurini
Status: Merged
Merge reported by: Garth Wells
Merged at revision: not available
Proposed branch: lp:~corrado-maurini/dolfin/tao
Merge into: lp:~fenics-core/dolfin/trunk
Diff against target: 6853 lines (+6648/-4)
23 files modified
CMakeLists.txt (+6/-1)
bench/logs/bench.log (+4/-0)
cmake/modules/FindSLEPc.cmake (+1/-1)
cmake/modules/FindTAO.cmake (+195/-0)
demo/undocumented/contact-vi-tao/cpp/CMakeLists.txt (+40/-0)
demo/undocumented/contact-vi-tao/cpp/Elasticity.h (+4785/-0)
demo/undocumented/contact-vi-tao/cpp/Elasticity.ufl (+45/-0)
demo/undocumented/contact-vi-tao/cpp/main.cpp (+125/-0)
demo/undocumented/contact-vi-tao/python/demo_contact.py (+93/-0)
dolfin/CMakeLists.txt (+7/-0)
dolfin/common/defines.cpp (+9/-0)
dolfin/common/defines.h (+3/-0)
dolfin/nls/TAOLinearBoundSolver.cpp (+440/-0)
dolfin/nls/TAOLinearBoundSolver.h (+210/-0)
dolfin/nls/dolfin_nls.h (+1/-0)
dolfin/swig/modules/la/dependencies.txt (+1/-1)
dolfin/swig/modules/la/dependencies.txt.moved (+1/-0)
dolfin/swig/modules/la/module.i (+2/-0)
dolfin/swig/nls/docstrings.i.THIS (+272/-0)
dolfin/swig/shared_ptr_classes.i (+3/-0)
site-packages/dolfin/compilemodules/swigimportinfo.py.THIS (+295/-0)
test/unit/nls/python/TAOLinearBoundSolver.py (+109/-0)
test/unit/test.py (+1/-1)
To merge this branch: bzr merge lp:~corrado-maurini/dolfin/tao
Reviewer Review Type Date Requested Status
Garth Wells very light Approve
Anders Logg Needs Fixing
Review via email: mp+123195@code.launchpad.net

Description of the change

Propose an interface to TAO for solving variational inequalities with bound constraints.

May you please have look and say if your are still interested in including it in dolfin and suggest possible improvements? The interface is minimal, but it allows to use most solvers and select underlying petsc KSP and PC, when possible. It considers only the minimization of quadratic bound-constrained functionals. If/when you think that it is almost suitable for merging I will do some cleanup and include an unit test.

The interface is implemented in TAOBoundConstrainedSolver.cpp. I added this file in /la.
The choice is questionable. I implemented it as a linear solver (get inspiration by PETScKrylovSolver), even if it is not and may have its place in nls.

To test:
----------------
see demo/la/bound-constraints. There are 3 examples in python and 1 in cpp

Dependences: requires to install TAO.
-------------------------------------
I tested it with petsc-3.2 and tao-2.0. Most probably everything will work also with petsc 3.3 and tao 2.1.

To download tao
 - http://www.mcs.anl.gov/research/projects/tao/download/tao-2.0-p6.tar.gz
 - to install just set TAO_DIR to the current tao directory and "make all" once the PETSC_DIR and PETSC_ARCH variables are set to your current petsc-3.2 configuration. Then I added a FindTAO.cmake module to check for it in dolfin. Please check this, because I am not an expert of cmake. I have just transposed the module for SLEPC.

To post a comment you must log in.
Revision history for this message
Johannes Ring (johannr) wrote :

Hi Corrado,

Have you signed the copyright consent form [1]? If not, please do that.

[1] http://fenicsproject.org/contributing/contributing_code.html#copyright-and-licensing-consent

Revision history for this message
corrado maurini (corrado-maurini) wrote :

Hi,

I will send you the author letter today.

Do you really need the one from my institution? This may take longer ...

Corrado

Le 13 sept. 2012 à 15:20, Johannes Ring a écrit :

> Hi Corrado,
>
> Have you signed the copyright consent form [1]? If not, please do that.
>
> [1] http://fenicsproject.org/contributing/contributing_code.html#copyright-and-licensing-consent
> --
> https://code.launchpad.net/~corrado-maurini/dolfin/tao/+merge/123195
> You are the owner of lp:~corrado-maurini/dolfin/tao.

Revision history for this message
Johannes Ring (johannr) wrote :

> I will send you the author letter today.

Thanks!

> Do you really need the one from my institution? This may take longer ...

Yes, we would like that and it is okay if it takes some time.

Revision history for this message
corrado maurini (corrado-maurini) wrote :

Hi,

After a long process, please find attached the agreement from my "organization".

The signature is not at the right place, but I hope that you may still accept it. It could require 2 further months to have a new version ...

If you confirm it is fine for you, I will send the paper version by mail.

Yours,

Corrado

Corrado Maurini
Institut Jean Le Rond d'Alembert (UMR 7190)
Université Pierre et Marie Curie
4 Place Jussieu, case 162
Tour 55-65, 4ème étage, Bureau 424
75252 Paris Cedex 05
Tél: 01 44 27 87 19, Fax: 01 44 27 52 59
E-mail: <email address hidden>
Web: http://www.lmm.jussieu.fr/~corrado/

Le 13 sept. 2012 à 18:13, Johannes Ring a écrit :

>> I will send you the author letter today.
>
> Thanks!
>
>> Do you really need the one from my institution? This may take longer ...
>
> Yes, we would like that and it is okay if it takes some time.
>
> --
> https://code.launchpad.net/~corrado-maurini/dolfin/tao/+merge/123195
> You are the owner of lp:~corrado-maurini/dolfin/tao.

Revision history for this message
Johannes Ring (johannr) wrote :

Hi Corrado,

Thanks, but I couldn't find the attachment. I think that when replying to a Launchpad merge request, any attachment will be lost. Please send it to me directly instead.

Revision history for this message
Anders Logg (logg) wrote :

Hi, what's the status of this patch? Is it still up for merge? If so, please sync it with trunk.

Other issues:

- Demos solving PDE should not go in the la directory
- Please add a unit test
- Consider adding matching C++ and Python versions for all demos

Revision history for this message
Anders Logg (logg) :
review: Needs Fixing (preliminary check)
Revision history for this message
corrado maurini (corrado-maurini) wrote :

> Hi, what's the status of this patch? Is it still up for merge? If so, please
> sync it with trunk.

done

>
> Other issues:
>
> - Demos solving PDE should not go in the la directory

moved to undocumented/conctact-vi-tao

> - Please add a unit test

done in nls/TAOLinearBoundSolver.py

> - Consider adding matching C++ and Python versions for all demos

done.

I have also move the TAOLinearBoundSolver class from la to nls, because probably more appropriate

I have also corrected here a bug: the name for the unit test of PETScSNESSolver in unit/test.py was wrong and the PETScSNESSolver was skipped.

Revision history for this message
Garth Wells (garth-wells) wrote :

I've merged this without testing (I don't have TAO installed) to get in before the git transition - it will get looked more closely later.

review: Approve (very light)

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'CMakeLists.txt'
2--- CMakeLists.txt 2013-03-25 09:21:56 +0000
3+++ CMakeLists.txt 2013-03-25 21:11:24 +0000
4@@ -66,6 +66,7 @@
5 list(APPEND OPTIONAL_PACKAGES "MPI")
6 list(APPEND OPTIONAL_PACKAGES "PETSc")
7 list(APPEND OPTIONAL_PACKAGES "SLEPc")
8+list(APPEND OPTIONAL_PACKAGES "TAO")
9 list(APPEND OPTIONAL_PACKAGES "Trilinos")
10 list(APPEND OPTIONAL_PACKAGES "UMFPACK")
11 list(APPEND OPTIONAL_PACKAGES "CHOLMOD")
12@@ -284,12 +285,15 @@
13
14 endif()
15
16-# Check for PETSc and SLEPc
17+# Check for PETSc, SLEPc and TAO
18 if (DOLFIN_ENABLE_PETSC)
19 find_package(PETSc 3.2)
20 if (PETSC_FOUND AND DOLFIN_ENABLE_SLEPC)
21 find_package(SLEPc 3.2)
22 endif()
23+ if (PETSC_FOUND AND DOLFIN_ENABLE_TAO)
24+ find_package(TAO)
25+ endif()
26 endif()
27
28 # Check for ParMETIS and SCOTCH
29@@ -406,6 +410,7 @@
30 endif()
31 endif()
32
33+
34 # Check for CGAL
35 if (DOLFIN_ENABLE_CGAL)
36 find_package(CGAL 4.1)
37
38=== modified file 'bench/logs/bench.log'
39--- bench/logs/bench.log 2010-06-21 08:56:48 +0000
40+++ bench/logs/bench.log 2013-03-25 21:11:24 +0000
41@@ -149,3 +149,7 @@
42 (2010, 6, 21, 1, 14, 44) mesh-iteration-cpp 35.36 "Iteration over entities of unit cube of size 128 x 128 x 128 (100 repetitions)"
43 (2010, 6, 21, 1, 15, 6) function-evaluation-cpp 21.8412 "Evaluations of functions at arbitrary points."
44 (2010, 6, 21, 1, 15, 6) function-evaluation-cpp 21.77 "Evaluations of functions at arbitrary points."
45+(2012, 9, 3, 0, 47, 55) fem-jit-python 0.000333905 "JIT compilation (in memory cache)"
46+(2012, 9, 3, 0, 48, 29) function-extrapolation-python 34.026 "Calling FFC just-in-time (JIT) compiler, this may take some time."
47+(2012, 9, 3, 0, 48, 29) function-extrapolation-python 1.19862 "Calling FFC just-in-time (JIT) compiler, this may take some time."
48+(2012, 9, 3, 0, 48, 30) la-cusp-python 0 "Solving the linear system for Poisson's equation using PETSc CUSP"
49
50=== modified file 'cmake/modules/FindSLEPc.cmake'
51--- cmake/modules/FindSLEPc.cmake 2013-02-21 12:11:19 +0000
52+++ cmake/modules/FindSLEPc.cmake 2013-03-25 21:11:24 +0000
53@@ -85,7 +85,7 @@
54
55 find_library(SLEPC_LIBRARY
56 NAMES slepc
57- HINTS ${SLEPC_DIR}/lib $ENV{SLEPC_DIR}/lib
58+ HINTS ${SLEPC_DIR}/lib $ENV{SLEPC_DIR}/lib ${SLEPC_DIR}/${PETSC_ARCH}/lib $ENV{SLEPC_DIR}/$ENV{PETSC_ARCH}/lib
59 DOC "The SLEPc library"
60 )
61 mark_as_advanced(SLEPC_LIBRARY)
62
63=== added file 'cmake/modules/FindTAO.cmake'
64--- cmake/modules/FindTAO.cmake 1970-01-01 00:00:00 +0000
65+++ cmake/modules/FindTAO.cmake 2013-03-25 21:11:24 +0000
66@@ -0,0 +1,195 @@
67+# - Try to find TAO
68+# Once done this will define
69+#
70+# TAO_FOUND - system has TAO
71+# TAO_INCLUDE_DIR - include directories for TAO
72+# TAO_LIBARIES - libraries for TAO
73+# TAO_DIR - directory where TAO is built
74+#
75+# Assumes that PETSC_DIR and PETSC_ARCH has been set by
76+# alredy calling find_package(PETSc)
77+
78+#=============================================================================
79+# Copyright (C) 2010-2012 Garth N. Wells, Anders Logg and Johannes Ring
80+# All rights reserved.
81+#
82+# Redistribution and use in source and binary forms, with or without
83+# modification, are permitted provided that the following conditions
84+# are met:
85+#
86+# 1. Redistributions of source code must retain the above copyright
87+# notice, this list of conditions and the following disclaimer.
88+# 2. Redistributions in binary form must reproduce the above copyright
89+# notice, this list of conditions and the following disclaimer in
90+# the documentation and/or other materials provided with the
91+# distribution.
92+#
93+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
94+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
95+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
96+# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
97+# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
98+# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
99+# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
100+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
101+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
102+# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
103+# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
104+# POSSIBILITY OF SUCH DAMAGE.
105+#=============================================================================
106+
107+message(STATUS "Checking for package 'TAO'")
108+
109+# Set debian_arches (PETSC_ARCH for Debian-style installations)
110+foreach (debian_arches linux kfreebsd)
111+ if ("${CMAKE_BUILD_TYPE}" STREQUAL "Debug")
112+ set(DEBIAN_FLAVORS ${debian_arches}-gnu-c-debug ${debian_arches}-gnu-c-opt ${DEBIAN_FLAVORS})
113+ else()
114+ set(DEBIAN_FLAVORS ${debian_arches}-gnu-c-opt ${debian_arches}-gnu-c-debug ${DEBIAN_FLAVORS})
115+ endif()
116+endforeach()
117+
118+# List of possible locations for TAO_DIR
119+set(tao_dir_locations "")
120+list(APPEND tao_dir_locations "/usr/local/lib/tao")
121+list(APPEND tao_dir_locations "$ENV{HOME}/tao")
122+list(APPEND tao_dir_locations "$ENV{TAO_DIR}")
123+list(APPEND tao_dir_locations "$ENV{TAO_DIR}")
124+list(APPEND tao_dir_locations "$ENV{TAO_DIR}/$ENV{PETSC_ARCH}")
125+
126+# Try to figure out TAO_DIR by finding tao.h
127+find_path(TAO_DIR /include/tao.h
128+ HINTS ${TAO_DIR} $ENV{TAO_DIR}
129+ PATHS ${tao_dir_locations}
130+ DOC "tao directory")
131+
132+#Report result of search for TAO_DIR
133+if (DEFINED TAO_DIR)
134+ message(STATUS "TAO_DIR is ${TAO_DIR}")
135+else()
136+ message(STATUS "TAO_DIR is empty")
137+endif()
138+
139+# Get variables from TAO configuration
140+if (TAO_DIR)
141+
142+ find_library(TAO_LIBRARY
143+ NAMES tao
144+ HINTS ${TAO_DIR} $ENV{TAO_DIR} ${TAO_DIR}/${PETSC_ARCH}/lib $ENV{TAO_DIR}/$ENV{PETSC_ARCH}/lib
145+ DOC "The TAO library"
146+ )
147+ mark_as_advanced(TAO_LIBRARY)
148+
149+ # Create a temporary Makefile to probe the TAOcc configuration
150+ set(tao_config_makefile ${PROJECT_BINARY_DIR}/Makefile.TAO)
151+ file(WRITE ${tao_config_makefile}
152+"# This file was autogenerated by FindTAO.cmake
153+TAO_DIR = ${TAO_DIR}
154+PETSC_ARCH = ${PETSC_ARCH}
155+PETSC_DIR = ${PETSC_DIR}
156+include ${TAO_DIR}/conf/tao_base
157+show :
158+ -@echo -n \${\${VARIABLE}}
159+")
160+
161+ # Define macro for getting TAO variables from Makefile
162+ macro(TAO_GET_VARIABLE var name)
163+ set(${var} "NOTFOUND" CACHE INTERNAL "Cleared" FORCE)
164+ execute_process(COMMAND ${CMAKE_MAKE_PROGRAM} --no-print-directory -f ${tao_config_makefile} show VARIABLE=${name}
165+ OUTPUT_VARIABLE ${var}
166+ RESULT_VARIABLE tao_return)
167+ endmacro()
168+
169+ # Call macro to get the TAO variables
170+ tao_get_variable(TAO_INCLUDE TAO_INCLUDE)
171+ tao_get_variable(TAO_LIB TAO_LIB)
172+
173+ # Remove temporary Makefile
174+ file(REMOVE ${tao_config_makefile})
175+
176+##############################################
177+ # Extract include paths and libraries from compile command line
178+ include(ResolveCompilerPaths)
179+ resolve_includes(TAO_INCLUDE_DIRS "${TAO_INCLUDE}")
180+ resolve_libraries(TAO_EXTERNAL_LIB "${TAO_EXTERNAL_LIB}")
181+
182+ # Add variables to CMake cache and mark as advanced
183+ set(TAO_INCLUDE_DIRS ${TAO_INCLUDE_DIRS} CACHE STRING "TAO include paths." FORCE)
184+ set(TAO_LIBRARIES ${TAO_LIBRARY} ${TAO_EXTERNAL_LIB} CACHE STRING "TAO libraries." FORCE)
185+ mark_as_advanced(TAO_INCLUDE_DIRS TAO_LIBRARIES)
186+endif()
187+
188+if (DOLFIN_SKIP_BUILD_TESTS)
189+ set(TAO_TEST_RUNS TRUE)
190+ set(TAO_VERSION "UNKNOWN")
191+ set(TAO_VERSION_OK TRUE)
192+elseif (TAO_LIBRARIES AND TAO_INCLUDE_DIRS)
193+
194+ # Set flags for building test program
195+ set(CMAKE_REQUIRED_INCLUDES ${TAO_INCLUDE_DIRS} ${PETSC_INCLUDE_DIRS})
196+ set(CMAKE_REQUIRED_LIBRARIES ${TAO_LIBRARIES} ${PETSC_LIBRARIES})
197+
198+ # Add MPI variables if MPI has been found
199+ if (MPI_C_FOUND)
200+ set(CMAKE_REQUIRED_INCLUDES ${CMAKE_REQUIRED_INCLUDES} ${MPI_C_INCLUDE_PATH})
201+
202+ set(CMAKE_REQUIRED_LIBRARIES ${CMAKE_REQUIRED_LIBRARIES} ${MPI_C_LIBRARIES})
203+ set(CMAKE_REQUIRED_FLAGS "${CMAKE_REQUIRED_FLAGS} ${MPI_C_COMPILE_FLAGS}")
204+ endif()
205+
206+##############################################
207+
208+ # Extract include paths and libraries from compile command line
209+ # include(ResolveCompilerPaths)
210+ # resolve_includes(TAO_INCLUDE_DIRS "${TAO_INCLUDE}")
211+ #resolve_libraries(TAO_LIB "${TAO_LIB}")
212+
213+ # Add variables to CMake cache and mark as advanced
214+ #set(TAO_INCLUDE_DIRS ${TAO_INCLUDE_DIRS} CACHE STRING "TAO include paths." FORCE)
215+ #set(TAO_LIBRARIES ${TAO_LIB} CACHE STRING "TAO libraries." FORCE)
216+ #mark_as_advanced(TAO INCLUDE_DIRS TAO_LIBRARIES)
217+
218+ # Set flags for building test program
219+ #set(CMAKE_REQUIRED_INCLUDES ${TAO_INCLUDE_DIRS} ${PETSC_INCLUDE_DIRS})
220+ #set(CMAKE_REQUIRED_LIBRARIES ${TAO_LIBRARIES} ${PETSC_LIBRARIES})
221+
222+ # Add MPI variables if MPI has been found
223+ #if (MPI_FOUND)
224+ # set(CMAKE_REQUIRED_INCLUDES ${CMAKE_REQUIRED_INCLUDES} ${MPI_INCLUDE_PATH})
225+ # set(CMAKE_REQUIRED_LIBRARIES ${CMAKE_REQUIRED_LIBRARIES} ${MPI_LIBRARIES})
226+ # set(CMAKE_REQUIRED_FLAGS "${CMAKE_REQUIRED_FLAGS} ${MPI_COMPILE_FLAGS}")
227+ #endif()
228+
229+
230+ # Run TAO test program
231+ include(CheckCXXSourceRuns)
232+ check_cxx_source_runs("
233+#include \"petsc.h\"
234+#include \"tao.h\"
235+int main()
236+{
237+ PetscErrorCode ierr;
238+ int argc = 0;
239+ char** argv = NULL;
240+ ierr = TaoInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
241+ TaoSolver tao;
242+ ierr = TaoCreate(PETSC_COMM_SELF, &tao); CHKERRQ(ierr);
243+ //ierr = TaoSetFromOptions(tao); CHKERRQ(ierr);
244+ ierr = TaoDestroy(&tao); CHKERRQ(ierr);
245+ ierr = TaoFinalize(); CHKERRQ(ierr);
246+ return 0;
247+}
248+" TAO_TEST_RUNS)
249+ if (TAO_TEST_RUNS)
250+ message(STATUS "TAO test runs")
251+ else()
252+ message(STATUS "TAO test failed")
253+ endif()
254+
255+endif()
256+
257+# Standard package handling
258+include(FindPackageHandleStandardArgs)
259+find_package_handle_standard_args(TAO
260+ "TAO could not be found. Be sure to set TAO_DIR, PETSC_DIR, and PETSC_ARCH."
261+ TAO_LIBRARIES TAO_DIR TAO_INCLUDE_DIRS TAO_TEST_RUNS)
262
263=== added directory 'demo/undocumented/contact-vi-tao'
264=== added directory 'demo/undocumented/contact-vi-tao/cpp'
265=== added file 'demo/undocumented/contact-vi-tao/cpp/CMakeLists.txt'
266--- demo/undocumented/contact-vi-tao/cpp/CMakeLists.txt 1970-01-01 00:00:00 +0000
267+++ demo/undocumented/contact-vi-tao/cpp/CMakeLists.txt 2013-03-25 21:11:24 +0000
268@@ -0,0 +1,40 @@
269+# This file is automatically generated by running
270+#
271+# cmake/scripts/generate-cmakefiles
272+#
273+# Require CMake 2.8
274+cmake_minimum_required(VERSION 2.8)
275+
276+set(PROJECT_NAME demo_contact)
277+project(${PROJECT_NAME})
278+
279+# Set verbose output while testing CMake
280+#set(CMAKE_VERBOSE_MAKEFILE 1)
281+
282+# Set CMake behavior
283+cmake_policy(SET CMP0004 OLD)
284+
285+# Get DOLFIN configuration data (DOLFINConfig.cmake must be in DOLFIN_CMAKE_CONFIG_PATH)
286+find_package(DOLFIN)
287+
288+# Default build type (can be overridden by user)
289+if (NOT CMAKE_BUILD_TYPE)
290+ set(CMAKE_BUILD_TYPE "RelWithDebInfo" CACHE STRING
291+ "Choose the type of build, options are: Debug MinSizeRel Release RelWithDebInfo." FORCE)
292+endif()
293+
294+# Compiler definitions
295+add_definitions(${DOLFIN_CXX_DEFINITIONS})
296+
297+# Add special DOLFIN compiler flags
298+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${DOLFIN_CXX_FLAGS}")
299+
300+# Include directories
301+include_directories(${DOLFIN_INCLUDE_DIRS})
302+include_directories(SYSTEM ${DOLFIN_3RD_PARTY_INCLUDE_DIRS})
303+
304+# Executable
305+add_executable(${PROJECT_NAME} main.cpp)
306+
307+# Target libraries
308+target_link_libraries(${PROJECT_NAME} ${DOLFIN_LIBRARIES} ${DOLFIN_3RD_PARTY_LIBRARIES})
309
310=== added file 'demo/undocumented/contact-vi-tao/cpp/Elasticity.h'
311--- demo/undocumented/contact-vi-tao/cpp/Elasticity.h 1970-01-01 00:00:00 +0000
312+++ demo/undocumented/contact-vi-tao/cpp/Elasticity.h 2013-03-25 21:11:24 +0000
313@@ -0,0 +1,4785 @@
314+// This code conforms with the UFC specification version 2.1.0+
315+// and was automatically generated by FFC version 1.1.0+.
316+//
317+// This code was generated with the option '-l dolfin' and
318+// contains DOLFIN-specific wrappers that depend on DOLFIN.
319+//
320+// This code was generated with the following parameters:
321+//
322+// cache_dir: ''
323+// convert_exceptions_to_warnings: False
324+// cpp_optimize: False
325+// cpp_optimize_flags: '-O2'
326+// epsilon: 1e-14
327+// error_control: False
328+// form_postfix: True
329+// format: 'dolfin'
330+// log_level: 20
331+// log_prefix: ''
332+// optimize: False
333+// output_dir: '.'
334+// precision: 15
335+// quadrature_degree: 'auto'
336+// quadrature_rule: 'auto'
337+// representation: 'auto'
338+// split: False
339+// swig_binary: 'swig'
340+// swig_path: ''
341+
342+#ifndef __ELASTICITY_H
343+#define __ELASTICITY_H
344+
345+#include <cmath>
346+#include <stdexcept>
347+#include <fstream>
348+#include <ufc.h>
349+
350+/// This class defines the interface for a finite element.
351+
352+class elasticity_finite_element_0: public ufc::finite_element
353+{
354+public:
355+
356+ /// Constructor
357+ elasticity_finite_element_0() : ufc::finite_element()
358+ {
359+ // Do nothing
360+ }
361+
362+ /// Destructor
363+ virtual ~elasticity_finite_element_0()
364+ {
365+ // Do nothing
366+ }
367+
368+ /// Return a string identifying the finite element
369+ virtual const char* signature() const
370+ {
371+ return "FiniteElement('Real', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 0, None)";
372+ }
373+
374+ /// Return the cell shape
375+ virtual ufc::shape cell_shape() const
376+ {
377+ return ufc::triangle;
378+ }
379+
380+ /// Return the topological dimension of the cell shape
381+ virtual std::size_t topological_dimension() const
382+ {
383+ return 2;
384+ }
385+
386+ /// Return the geometric dimension of the cell shape
387+ virtual std::size_t geometric_dimension() const
388+ {
389+ return 2;
390+ }
391+
392+ /// Return the dimension of the finite element function space
393+ virtual std::size_t space_dimension() const
394+ {
395+ return 1;
396+ }
397+
398+ /// Return the rank of the value space
399+ virtual std::size_t value_rank() const
400+ {
401+ return 0;
402+ }
403+
404+ /// Return the dimension of the value space for axis i
405+ virtual std::size_t value_dimension(std::size_t i) const
406+ {
407+ return 1;
408+ }
409+
410+ /// Evaluate basis function i at given point x in cell
411+ virtual void evaluate_basis(std::size_t i,
412+ double* values,
413+ const double* x,
414+ const double* vertex_coordinates,
415+ int cell_orientation) const
416+ {
417+ // Compute Jacobian
418+ double J[4];
419+ compute_jacobian_triangle_2d(J, vertex_coordinates);
420+
421+ // Compute Jacobian inverse and determinant
422+ double K[4];
423+ double detJ;
424+ compute_jacobian_inverse_triangle_2d(K, detJ, J);
425+
426+
427+ // Compute constants
428+
429+ // Get coordinates and map to the reference (FIAT) element
430+
431+ // Reset values
432+ *values = 0.0;
433+
434+ // Array of basisvalues
435+ double basisvalues[1] = {0.0};
436+
437+ // Declare helper variables
438+
439+ // Compute basisvalues
440+ basisvalues[0] = 1.0;
441+
442+ // Table(s) of coefficients
443+ static const double coefficients0[1] = \
444+ {1.0};
445+
446+ // Compute value(s)
447+ for (unsigned int r = 0; r < 1; r++)
448+ {
449+ *values += coefficients0[r]*basisvalues[r];
450+ }// end loop over 'r'
451+ }
452+
453+ /// Evaluate all basis functions at given point x in cell
454+ virtual void evaluate_basis_all(double* values,
455+ const double* x,
456+ const double* vertex_coordinates,
457+ int cell_orientation) const
458+ {
459+ // Element is constant, calling evaluate_basis.
460+ evaluate_basis(0, values, x, vertex_coordinates, cell_orientation);
461+ }
462+
463+ /// Evaluate order n derivatives of basis function i at given point x in cell
464+ virtual void evaluate_basis_derivatives(std::size_t i,
465+ std::size_t n,
466+ double* values,
467+ const double* x,
468+ const double* vertex_coordinates,
469+ int cell_orientation) const
470+ {
471+ // Compute Jacobian
472+ double J[4];
473+ compute_jacobian_triangle_2d(J, vertex_coordinates);
474+
475+ // Compute Jacobian inverse and determinant
476+ double K[4];
477+ double detJ;
478+ compute_jacobian_inverse_triangle_2d(K, detJ, J);
479+
480+
481+ // Compute constants
482+
483+ // Get coordinates and map to the reference (FIAT) element
484+
485+ // Compute number of derivatives.
486+ unsigned int num_derivatives = 1;
487+ for (unsigned int r = 0; r < n; r++)
488+ {
489+ num_derivatives *= 2;
490+ }// end loop over 'r'
491+
492+ // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
493+ unsigned int **combinations = new unsigned int *[num_derivatives];
494+ for (unsigned int row = 0; row < num_derivatives; row++)
495+ {
496+ combinations[row] = new unsigned int [n];
497+ for (unsigned int col = 0; col < n; col++)
498+ combinations[row][col] = 0;
499+ }
500+
501+ // Generate combinations of derivatives
502+ for (unsigned int row = 1; row < num_derivatives; row++)
503+ {
504+ for (unsigned int num = 0; num < row; num++)
505+ {
506+ for (unsigned int col = n-1; col+1 > 0; col--)
507+ {
508+ if (combinations[row][col] + 1 > 1)
509+ combinations[row][col] = 0;
510+ else
511+ {
512+ combinations[row][col] += 1;
513+ break;
514+ }
515+ }
516+ }
517+ }
518+
519+ // Compute inverse of Jacobian
520+ const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
521+
522+ // Declare transformation matrix
523+ // Declare pointer to two dimensional array and initialise
524+ double **transform = new double *[num_derivatives];
525+
526+ for (unsigned int j = 0; j < num_derivatives; j++)
527+ {
528+ transform[j] = new double [num_derivatives];
529+ for (unsigned int k = 0; k < num_derivatives; k++)
530+ transform[j][k] = 1;
531+ }
532+
533+ // Construct transformation matrix
534+ for (unsigned int row = 0; row < num_derivatives; row++)
535+ {
536+ for (unsigned int col = 0; col < num_derivatives; col++)
537+ {
538+ for (unsigned int k = 0; k < n; k++)
539+ transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
540+ }
541+ }
542+
543+ // Reset values. Assuming that values is always an array.
544+ for (unsigned int r = 0; r < num_derivatives; r++)
545+ {
546+ values[r] = 0.0;
547+ }// end loop over 'r'
548+
549+
550+ // Array of basisvalues
551+ double basisvalues[1] = {0.0};
552+
553+ // Declare helper variables
554+
555+ // Compute basisvalues
556+ basisvalues[0] = 1.0;
557+
558+ // Table(s) of coefficients
559+ static const double coefficients0[1] = \
560+ {1.0};
561+
562+ // Tables of derivatives of the polynomial base (transpose).
563+ static const double dmats0[1][1] = \
564+ {{0.0}};
565+
566+ static const double dmats1[1][1] = \
567+ {{0.0}};
568+
569+ // Compute reference derivatives.
570+ // Declare pointer to array of derivatives on FIAT element.
571+ double *derivatives = new double[num_derivatives];
572+ for (unsigned int r = 0; r < num_derivatives; r++)
573+ {
574+ derivatives[r] = 0.0;
575+ }// end loop over 'r'
576+
577+ // Declare derivative matrix (of polynomial basis).
578+ double dmats[1][1] = \
579+ {{1.0}};
580+
581+ // Declare (auxiliary) derivative matrix (of polynomial basis).
582+ double dmats_old[1][1] = \
583+ {{1.0}};
584+
585+ // Loop possible derivatives.
586+ for (unsigned int r = 0; r < num_derivatives; r++)
587+ {
588+ // Resetting dmats values to compute next derivative.
589+ for (unsigned int t = 0; t < 1; t++)
590+ {
591+ for (unsigned int u = 0; u < 1; u++)
592+ {
593+ dmats[t][u] = 0.0;
594+ if (t == u)
595+ {
596+ dmats[t][u] = 1.0;
597+ }
598+
599+ }// end loop over 'u'
600+ }// end loop over 't'
601+
602+ // Looping derivative order to generate dmats.
603+ for (unsigned int s = 0; s < n; s++)
604+ {
605+ // Updating dmats_old with new values and resetting dmats.
606+ for (unsigned int t = 0; t < 1; t++)
607+ {
608+ for (unsigned int u = 0; u < 1; u++)
609+ {
610+ dmats_old[t][u] = dmats[t][u];
611+ dmats[t][u] = 0.0;
612+ }// end loop over 'u'
613+ }// end loop over 't'
614+
615+ // Update dmats using an inner product.
616+ if (combinations[r][s] == 0)
617+ {
618+ for (unsigned int t = 0; t < 1; t++)
619+ {
620+ for (unsigned int u = 0; u < 1; u++)
621+ {
622+ for (unsigned int tu = 0; tu < 1; tu++)
623+ {
624+ dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
625+ }// end loop over 'tu'
626+ }// end loop over 'u'
627+ }// end loop over 't'
628+ }
629+
630+ if (combinations[r][s] == 1)
631+ {
632+ for (unsigned int t = 0; t < 1; t++)
633+ {
634+ for (unsigned int u = 0; u < 1; u++)
635+ {
636+ for (unsigned int tu = 0; tu < 1; tu++)
637+ {
638+ dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
639+ }// end loop over 'tu'
640+ }// end loop over 'u'
641+ }// end loop over 't'
642+ }
643+
644+ }// end loop over 's'
645+ for (unsigned int s = 0; s < 1; s++)
646+ {
647+ for (unsigned int t = 0; t < 1; t++)
648+ {
649+ derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
650+ }// end loop over 't'
651+ }// end loop over 's'
652+ }// end loop over 'r'
653+
654+ // Transform derivatives back to physical element
655+ for (unsigned int r = 0; r < num_derivatives; r++)
656+ {
657+ for (unsigned int s = 0; s < num_derivatives; s++)
658+ {
659+ values[r] += transform[r][s]*derivatives[s];
660+ }// end loop over 's'
661+ }// end loop over 'r'
662+
663+ // Delete pointer to array of derivatives on FIAT element
664+ delete [] derivatives;
665+
666+ // Delete pointer to array of combinations of derivatives and transform
667+ for (unsigned int r = 0; r < num_derivatives; r++)
668+ {
669+ delete [] combinations[r];
670+ }// end loop over 'r'
671+ delete [] combinations;
672+ for (unsigned int r = 0; r < num_derivatives; r++)
673+ {
674+ delete [] transform[r];
675+ }// end loop over 'r'
676+ delete [] transform;
677+ }
678+
679+ /// Evaluate order n derivatives of all basis functions at given point x in cell
680+ virtual void evaluate_basis_derivatives_all(std::size_t n,
681+ double* values,
682+ const double* x,
683+ const double* vertex_coordinates,
684+ int cell_orientation) const
685+ {
686+ // Element is constant, calling evaluate_basis_derivatives.
687+ evaluate_basis_derivatives(0, n, values, x, vertex_coordinates, cell_orientation);
688+ }
689+
690+ /// Evaluate linear functional for dof i on the function f
691+ virtual double evaluate_dof(std::size_t i,
692+ const ufc::function& f,
693+ const double* vertex_coordinates,
694+ int cell_orientation,
695+ const ufc::cell& c) const
696+ {
697+ // Declare variables for result of evaluation
698+ double vals[1];
699+
700+ // Declare variable for physical coordinates
701+ double y[2];
702+ switch (i)
703+ {
704+ case 0:
705+ {
706+ y[0] = 0.333333333333333*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[2] + 0.333333333333333*vertex_coordinates[4];
707+ y[1] = 0.333333333333333*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[3] + 0.333333333333333*vertex_coordinates[5];
708+ f.evaluate(vals, y, c);
709+ return vals[0];
710+ break;
711+ }
712+ }
713+
714+ return 0.0;
715+ }
716+
717+ /// Evaluate linear functionals for all dofs on the function f
718+ virtual void evaluate_dofs(double* values,
719+ const ufc::function& f,
720+ const double* vertex_coordinates,
721+ int cell_orientation,
722+ const ufc::cell& c) const
723+ {
724+ // Declare variables for result of evaluation
725+ double vals[1];
726+
727+ // Declare variable for physical coordinates
728+ double y[2];
729+ y[0] = 0.333333333333333*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[2] + 0.333333333333333*vertex_coordinates[4];
730+ y[1] = 0.333333333333333*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[3] + 0.333333333333333*vertex_coordinates[5];
731+ f.evaluate(vals, y, c);
732+ values[0] = vals[0];
733+ }
734+
735+ /// Interpolate vertex values from dof values
736+ virtual void interpolate_vertex_values(double* vertex_values,
737+ const double* dof_values,
738+ const double* vertex_coordinates,
739+ int cell_orientation,
740+ const ufc::cell& c) const
741+ {
742+ // Evaluate function and change variables
743+ vertex_values[0] = dof_values[0];
744+ vertex_values[1] = dof_values[0];
745+ vertex_values[2] = dof_values[0];
746+ }
747+
748+ /// Map coordinate xhat from reference cell to coordinate x in cell
749+ virtual void map_from_reference_cell(double* x,
750+ const double* xhat,
751+ const ufc::cell& c) const
752+ {
753+ throw std::runtime_error("map_from_reference_cell not yet implemented.");
754+ }
755+
756+ /// Map from coordinate x in cell to coordinate xhat in reference cell
757+ virtual void map_to_reference_cell(double* xhat,
758+ const double* x,
759+ const ufc::cell& c) const
760+ {
761+ throw std::runtime_error("map_to_reference_cell not yet implemented.");
762+ }
763+
764+ /// Return the number of sub elements (for a mixed element)
765+ virtual std::size_t num_sub_elements() const
766+ {
767+ return 0;
768+ }
769+
770+ /// Create a new finite element for sub element i (for a mixed element)
771+ virtual ufc::finite_element* create_sub_element(std::size_t i) const
772+ {
773+ return 0;
774+ }
775+
776+ /// Create a new class instance
777+ virtual ufc::finite_element* create() const
778+ {
779+ return new elasticity_finite_element_0();
780+ }
781+
782+};
783+
784+/// This class defines the interface for a finite element.
785+
786+class elasticity_finite_element_1: public ufc::finite_element
787+{
788+public:
789+
790+ /// Constructor
791+ elasticity_finite_element_1() : ufc::finite_element()
792+ {
793+ // Do nothing
794+ }
795+
796+ /// Destructor
797+ virtual ~elasticity_finite_element_1()
798+ {
799+ // Do nothing
800+ }
801+
802+ /// Return a string identifying the finite element
803+ virtual const char* signature() const
804+ {
805+ return "FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
806+ }
807+
808+ /// Return the cell shape
809+ virtual ufc::shape cell_shape() const
810+ {
811+ return ufc::triangle;
812+ }
813+
814+ /// Return the topological dimension of the cell shape
815+ virtual std::size_t topological_dimension() const
816+ {
817+ return 2;
818+ }
819+
820+ /// Return the geometric dimension of the cell shape
821+ virtual std::size_t geometric_dimension() const
822+ {
823+ return 2;
824+ }
825+
826+ /// Return the dimension of the finite element function space
827+ virtual std::size_t space_dimension() const
828+ {
829+ return 3;
830+ }
831+
832+ /// Return the rank of the value space
833+ virtual std::size_t value_rank() const
834+ {
835+ return 0;
836+ }
837+
838+ /// Return the dimension of the value space for axis i
839+ virtual std::size_t value_dimension(std::size_t i) const
840+ {
841+ return 1;
842+ }
843+
844+ /// Evaluate basis function i at given point x in cell
845+ virtual void evaluate_basis(std::size_t i,
846+ double* values,
847+ const double* x,
848+ const double* vertex_coordinates,
849+ int cell_orientation) const
850+ {
851+ // Compute Jacobian
852+ double J[4];
853+ compute_jacobian_triangle_2d(J, vertex_coordinates);
854+
855+ // Compute Jacobian inverse and determinant
856+ double K[4];
857+ double detJ;
858+ compute_jacobian_inverse_triangle_2d(K, detJ, J);
859+
860+
861+ // Compute constants
862+ const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
863+ const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
864+
865+ // Get coordinates and map to the reference (FIAT) element
866+ double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
867+ double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
868+
869+ // Reset values
870+ *values = 0.0;
871+ switch (i)
872+ {
873+ case 0:
874+ {
875+
876+ // Array of basisvalues
877+ double basisvalues[3] = {0.0, 0.0, 0.0};
878+
879+ // Declare helper variables
880+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
881+
882+ // Compute basisvalues
883+ basisvalues[0] = 1.0;
884+ basisvalues[1] = tmp0;
885+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
886+ basisvalues[0] *= std::sqrt(0.5);
887+ basisvalues[2] *= std::sqrt(1.0);
888+ basisvalues[1] *= std::sqrt(3.0);
889+
890+ // Table(s) of coefficients
891+ static const double coefficients0[3] = \
892+ {0.471404520791032, -0.288675134594813, -0.166666666666667};
893+
894+ // Compute value(s)
895+ for (unsigned int r = 0; r < 3; r++)
896+ {
897+ *values += coefficients0[r]*basisvalues[r];
898+ }// end loop over 'r'
899+ break;
900+ }
901+ case 1:
902+ {
903+
904+ // Array of basisvalues
905+ double basisvalues[3] = {0.0, 0.0, 0.0};
906+
907+ // Declare helper variables
908+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
909+
910+ // Compute basisvalues
911+ basisvalues[0] = 1.0;
912+ basisvalues[1] = tmp0;
913+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
914+ basisvalues[0] *= std::sqrt(0.5);
915+ basisvalues[2] *= std::sqrt(1.0);
916+ basisvalues[1] *= std::sqrt(3.0);
917+
918+ // Table(s) of coefficients
919+ static const double coefficients0[3] = \
920+ {0.471404520791032, 0.288675134594813, -0.166666666666667};
921+
922+ // Compute value(s)
923+ for (unsigned int r = 0; r < 3; r++)
924+ {
925+ *values += coefficients0[r]*basisvalues[r];
926+ }// end loop over 'r'
927+ break;
928+ }
929+ case 2:
930+ {
931+
932+ // Array of basisvalues
933+ double basisvalues[3] = {0.0, 0.0, 0.0};
934+
935+ // Declare helper variables
936+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
937+
938+ // Compute basisvalues
939+ basisvalues[0] = 1.0;
940+ basisvalues[1] = tmp0;
941+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
942+ basisvalues[0] *= std::sqrt(0.5);
943+ basisvalues[2] *= std::sqrt(1.0);
944+ basisvalues[1] *= std::sqrt(3.0);
945+
946+ // Table(s) of coefficients
947+ static const double coefficients0[3] = \
948+ {0.471404520791032, 0.0, 0.333333333333333};
949+
950+ // Compute value(s)
951+ for (unsigned int r = 0; r < 3; r++)
952+ {
953+ *values += coefficients0[r]*basisvalues[r];
954+ }// end loop over 'r'
955+ break;
956+ }
957+ }
958+
959+ }
960+
961+ /// Evaluate all basis functions at given point x in cell
962+ virtual void evaluate_basis_all(double* values,
963+ const double* x,
964+ const double* vertex_coordinates,
965+ int cell_orientation) const
966+ {
967+ // Helper variable to hold values of a single dof.
968+ double dof_values = 0.0;
969+
970+ // Loop dofs and call evaluate_basis
971+ for (unsigned int r = 0; r < 3; r++)
972+ {
973+ evaluate_basis(r, &dof_values, x, vertex_coordinates, cell_orientation);
974+ values[r] = dof_values;
975+ }// end loop over 'r'
976+ }
977+
978+ /// Evaluate order n derivatives of basis function i at given point x in cell
979+ virtual void evaluate_basis_derivatives(std::size_t i,
980+ std::size_t n,
981+ double* values,
982+ const double* x,
983+ const double* vertex_coordinates,
984+ int cell_orientation) const
985+ {
986+ // Compute Jacobian
987+ double J[4];
988+ compute_jacobian_triangle_2d(J, vertex_coordinates);
989+
990+ // Compute Jacobian inverse and determinant
991+ double K[4];
992+ double detJ;
993+ compute_jacobian_inverse_triangle_2d(K, detJ, J);
994+
995+
996+ // Compute constants
997+ const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
998+ const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
999+
1000+ // Get coordinates and map to the reference (FIAT) element
1001+ double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
1002+ double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1003+
1004+ // Compute number of derivatives.
1005+ unsigned int num_derivatives = 1;
1006+ for (unsigned int r = 0; r < n; r++)
1007+ {
1008+ num_derivatives *= 2;
1009+ }// end loop over 'r'
1010+
1011+ // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1012+ unsigned int **combinations = new unsigned int *[num_derivatives];
1013+ for (unsigned int row = 0; row < num_derivatives; row++)
1014+ {
1015+ combinations[row] = new unsigned int [n];
1016+ for (unsigned int col = 0; col < n; col++)
1017+ combinations[row][col] = 0;
1018+ }
1019+
1020+ // Generate combinations of derivatives
1021+ for (unsigned int row = 1; row < num_derivatives; row++)
1022+ {
1023+ for (unsigned int num = 0; num < row; num++)
1024+ {
1025+ for (unsigned int col = n-1; col+1 > 0; col--)
1026+ {
1027+ if (combinations[row][col] + 1 > 1)
1028+ combinations[row][col] = 0;
1029+ else
1030+ {
1031+ combinations[row][col] += 1;
1032+ break;
1033+ }
1034+ }
1035+ }
1036+ }
1037+
1038+ // Compute inverse of Jacobian
1039+ const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
1040+
1041+ // Declare transformation matrix
1042+ // Declare pointer to two dimensional array and initialise
1043+ double **transform = new double *[num_derivatives];
1044+
1045+ for (unsigned int j = 0; j < num_derivatives; j++)
1046+ {
1047+ transform[j] = new double [num_derivatives];
1048+ for (unsigned int k = 0; k < num_derivatives; k++)
1049+ transform[j][k] = 1;
1050+ }
1051+
1052+ // Construct transformation matrix
1053+ for (unsigned int row = 0; row < num_derivatives; row++)
1054+ {
1055+ for (unsigned int col = 0; col < num_derivatives; col++)
1056+ {
1057+ for (unsigned int k = 0; k < n; k++)
1058+ transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
1059+ }
1060+ }
1061+
1062+ // Reset values. Assuming that values is always an array.
1063+ for (unsigned int r = 0; r < num_derivatives; r++)
1064+ {
1065+ values[r] = 0.0;
1066+ }// end loop over 'r'
1067+
1068+ switch (i)
1069+ {
1070+ case 0:
1071+ {
1072+
1073+ // Array of basisvalues
1074+ double basisvalues[3] = {0.0, 0.0, 0.0};
1075+
1076+ // Declare helper variables
1077+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1078+
1079+ // Compute basisvalues
1080+ basisvalues[0] = 1.0;
1081+ basisvalues[1] = tmp0;
1082+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1083+ basisvalues[0] *= std::sqrt(0.5);
1084+ basisvalues[2] *= std::sqrt(1.0);
1085+ basisvalues[1] *= std::sqrt(3.0);
1086+
1087+ // Table(s) of coefficients
1088+ static const double coefficients0[3] = \
1089+ {0.471404520791032, -0.288675134594813, -0.166666666666667};
1090+
1091+ // Tables of derivatives of the polynomial base (transpose).
1092+ static const double dmats0[3][3] = \
1093+ {{0.0, 0.0, 0.0},
1094+ {4.89897948556636, 0.0, 0.0},
1095+ {0.0, 0.0, 0.0}};
1096+
1097+ static const double dmats1[3][3] = \
1098+ {{0.0, 0.0, 0.0},
1099+ {2.44948974278318, 0.0, 0.0},
1100+ {4.24264068711928, 0.0, 0.0}};
1101+
1102+ // Compute reference derivatives.
1103+ // Declare pointer to array of derivatives on FIAT element.
1104+ double *derivatives = new double[num_derivatives];
1105+ for (unsigned int r = 0; r < num_derivatives; r++)
1106+ {
1107+ derivatives[r] = 0.0;
1108+ }// end loop over 'r'
1109+
1110+ // Declare derivative matrix (of polynomial basis).
1111+ double dmats[3][3] = \
1112+ {{1.0, 0.0, 0.0},
1113+ {0.0, 1.0, 0.0},
1114+ {0.0, 0.0, 1.0}};
1115+
1116+ // Declare (auxiliary) derivative matrix (of polynomial basis).
1117+ double dmats_old[3][3] = \
1118+ {{1.0, 0.0, 0.0},
1119+ {0.0, 1.0, 0.0},
1120+ {0.0, 0.0, 1.0}};
1121+
1122+ // Loop possible derivatives.
1123+ for (unsigned int r = 0; r < num_derivatives; r++)
1124+ {
1125+ // Resetting dmats values to compute next derivative.
1126+ for (unsigned int t = 0; t < 3; t++)
1127+ {
1128+ for (unsigned int u = 0; u < 3; u++)
1129+ {
1130+ dmats[t][u] = 0.0;
1131+ if (t == u)
1132+ {
1133+ dmats[t][u] = 1.0;
1134+ }
1135+
1136+ }// end loop over 'u'
1137+ }// end loop over 't'
1138+
1139+ // Looping derivative order to generate dmats.
1140+ for (unsigned int s = 0; s < n; s++)
1141+ {
1142+ // Updating dmats_old with new values and resetting dmats.
1143+ for (unsigned int t = 0; t < 3; t++)
1144+ {
1145+ for (unsigned int u = 0; u < 3; u++)
1146+ {
1147+ dmats_old[t][u] = dmats[t][u];
1148+ dmats[t][u] = 0.0;
1149+ }// end loop over 'u'
1150+ }// end loop over 't'
1151+
1152+ // Update dmats using an inner product.
1153+ if (combinations[r][s] == 0)
1154+ {
1155+ for (unsigned int t = 0; t < 3; t++)
1156+ {
1157+ for (unsigned int u = 0; u < 3; u++)
1158+ {
1159+ for (unsigned int tu = 0; tu < 3; tu++)
1160+ {
1161+ dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1162+ }// end loop over 'tu'
1163+ }// end loop over 'u'
1164+ }// end loop over 't'
1165+ }
1166+
1167+ if (combinations[r][s] == 1)
1168+ {
1169+ for (unsigned int t = 0; t < 3; t++)
1170+ {
1171+ for (unsigned int u = 0; u < 3; u++)
1172+ {
1173+ for (unsigned int tu = 0; tu < 3; tu++)
1174+ {
1175+ dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1176+ }// end loop over 'tu'
1177+ }// end loop over 'u'
1178+ }// end loop over 't'
1179+ }
1180+
1181+ }// end loop over 's'
1182+ for (unsigned int s = 0; s < 3; s++)
1183+ {
1184+ for (unsigned int t = 0; t < 3; t++)
1185+ {
1186+ derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1187+ }// end loop over 't'
1188+ }// end loop over 's'
1189+ }// end loop over 'r'
1190+
1191+ // Transform derivatives back to physical element
1192+ for (unsigned int r = 0; r < num_derivatives; r++)
1193+ {
1194+ for (unsigned int s = 0; s < num_derivatives; s++)
1195+ {
1196+ values[r] += transform[r][s]*derivatives[s];
1197+ }// end loop over 's'
1198+ }// end loop over 'r'
1199+
1200+ // Delete pointer to array of derivatives on FIAT element
1201+ delete [] derivatives;
1202+
1203+ // Delete pointer to array of combinations of derivatives and transform
1204+ for (unsigned int r = 0; r < num_derivatives; r++)
1205+ {
1206+ delete [] combinations[r];
1207+ }// end loop over 'r'
1208+ delete [] combinations;
1209+ for (unsigned int r = 0; r < num_derivatives; r++)
1210+ {
1211+ delete [] transform[r];
1212+ }// end loop over 'r'
1213+ delete [] transform;
1214+ break;
1215+ }
1216+ case 1:
1217+ {
1218+
1219+ // Array of basisvalues
1220+ double basisvalues[3] = {0.0, 0.0, 0.0};
1221+
1222+ // Declare helper variables
1223+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1224+
1225+ // Compute basisvalues
1226+ basisvalues[0] = 1.0;
1227+ basisvalues[1] = tmp0;
1228+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1229+ basisvalues[0] *= std::sqrt(0.5);
1230+ basisvalues[2] *= std::sqrt(1.0);
1231+ basisvalues[1] *= std::sqrt(3.0);
1232+
1233+ // Table(s) of coefficients
1234+ static const double coefficients0[3] = \
1235+ {0.471404520791032, 0.288675134594813, -0.166666666666667};
1236+
1237+ // Tables of derivatives of the polynomial base (transpose).
1238+ static const double dmats0[3][3] = \
1239+ {{0.0, 0.0, 0.0},
1240+ {4.89897948556636, 0.0, 0.0},
1241+ {0.0, 0.0, 0.0}};
1242+
1243+ static const double dmats1[3][3] = \
1244+ {{0.0, 0.0, 0.0},
1245+ {2.44948974278318, 0.0, 0.0},
1246+ {4.24264068711928, 0.0, 0.0}};
1247+
1248+ // Compute reference derivatives.
1249+ // Declare pointer to array of derivatives on FIAT element.
1250+ double *derivatives = new double[num_derivatives];
1251+ for (unsigned int r = 0; r < num_derivatives; r++)
1252+ {
1253+ derivatives[r] = 0.0;
1254+ }// end loop over 'r'
1255+
1256+ // Declare derivative matrix (of polynomial basis).
1257+ double dmats[3][3] = \
1258+ {{1.0, 0.0, 0.0},
1259+ {0.0, 1.0, 0.0},
1260+ {0.0, 0.0, 1.0}};
1261+
1262+ // Declare (auxiliary) derivative matrix (of polynomial basis).
1263+ double dmats_old[3][3] = \
1264+ {{1.0, 0.0, 0.0},
1265+ {0.0, 1.0, 0.0},
1266+ {0.0, 0.0, 1.0}};
1267+
1268+ // Loop possible derivatives.
1269+ for (unsigned int r = 0; r < num_derivatives; r++)
1270+ {
1271+ // Resetting dmats values to compute next derivative.
1272+ for (unsigned int t = 0; t < 3; t++)
1273+ {
1274+ for (unsigned int u = 0; u < 3; u++)
1275+ {
1276+ dmats[t][u] = 0.0;
1277+ if (t == u)
1278+ {
1279+ dmats[t][u] = 1.0;
1280+ }
1281+
1282+ }// end loop over 'u'
1283+ }// end loop over 't'
1284+
1285+ // Looping derivative order to generate dmats.
1286+ for (unsigned int s = 0; s < n; s++)
1287+ {
1288+ // Updating dmats_old with new values and resetting dmats.
1289+ for (unsigned int t = 0; t < 3; t++)
1290+ {
1291+ for (unsigned int u = 0; u < 3; u++)
1292+ {
1293+ dmats_old[t][u] = dmats[t][u];
1294+ dmats[t][u] = 0.0;
1295+ }// end loop over 'u'
1296+ }// end loop over 't'
1297+
1298+ // Update dmats using an inner product.
1299+ if (combinations[r][s] == 0)
1300+ {
1301+ for (unsigned int t = 0; t < 3; t++)
1302+ {
1303+ for (unsigned int u = 0; u < 3; u++)
1304+ {
1305+ for (unsigned int tu = 0; tu < 3; tu++)
1306+ {
1307+ dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1308+ }// end loop over 'tu'
1309+ }// end loop over 'u'
1310+ }// end loop over 't'
1311+ }
1312+
1313+ if (combinations[r][s] == 1)
1314+ {
1315+ for (unsigned int t = 0; t < 3; t++)
1316+ {
1317+ for (unsigned int u = 0; u < 3; u++)
1318+ {
1319+ for (unsigned int tu = 0; tu < 3; tu++)
1320+ {
1321+ dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1322+ }// end loop over 'tu'
1323+ }// end loop over 'u'
1324+ }// end loop over 't'
1325+ }
1326+
1327+ }// end loop over 's'
1328+ for (unsigned int s = 0; s < 3; s++)
1329+ {
1330+ for (unsigned int t = 0; t < 3; t++)
1331+ {
1332+ derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1333+ }// end loop over 't'
1334+ }// end loop over 's'
1335+ }// end loop over 'r'
1336+
1337+ // Transform derivatives back to physical element
1338+ for (unsigned int r = 0; r < num_derivatives; r++)
1339+ {
1340+ for (unsigned int s = 0; s < num_derivatives; s++)
1341+ {
1342+ values[r] += transform[r][s]*derivatives[s];
1343+ }// end loop over 's'
1344+ }// end loop over 'r'
1345+
1346+ // Delete pointer to array of derivatives on FIAT element
1347+ delete [] derivatives;
1348+
1349+ // Delete pointer to array of combinations of derivatives and transform
1350+ for (unsigned int r = 0; r < num_derivatives; r++)
1351+ {
1352+ delete [] combinations[r];
1353+ }// end loop over 'r'
1354+ delete [] combinations;
1355+ for (unsigned int r = 0; r < num_derivatives; r++)
1356+ {
1357+ delete [] transform[r];
1358+ }// end loop over 'r'
1359+ delete [] transform;
1360+ break;
1361+ }
1362+ case 2:
1363+ {
1364+
1365+ // Array of basisvalues
1366+ double basisvalues[3] = {0.0, 0.0, 0.0};
1367+
1368+ // Declare helper variables
1369+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1370+
1371+ // Compute basisvalues
1372+ basisvalues[0] = 1.0;
1373+ basisvalues[1] = tmp0;
1374+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1375+ basisvalues[0] *= std::sqrt(0.5);
1376+ basisvalues[2] *= std::sqrt(1.0);
1377+ basisvalues[1] *= std::sqrt(3.0);
1378+
1379+ // Table(s) of coefficients
1380+ static const double coefficients0[3] = \
1381+ {0.471404520791032, 0.0, 0.333333333333333};
1382+
1383+ // Tables of derivatives of the polynomial base (transpose).
1384+ static const double dmats0[3][3] = \
1385+ {{0.0, 0.0, 0.0},
1386+ {4.89897948556636, 0.0, 0.0},
1387+ {0.0, 0.0, 0.0}};
1388+
1389+ static const double dmats1[3][3] = \
1390+ {{0.0, 0.0, 0.0},
1391+ {2.44948974278318, 0.0, 0.0},
1392+ {4.24264068711928, 0.0, 0.0}};
1393+
1394+ // Compute reference derivatives.
1395+ // Declare pointer to array of derivatives on FIAT element.
1396+ double *derivatives = new double[num_derivatives];
1397+ for (unsigned int r = 0; r < num_derivatives; r++)
1398+ {
1399+ derivatives[r] = 0.0;
1400+ }// end loop over 'r'
1401+
1402+ // Declare derivative matrix (of polynomial basis).
1403+ double dmats[3][3] = \
1404+ {{1.0, 0.0, 0.0},
1405+ {0.0, 1.0, 0.0},
1406+ {0.0, 0.0, 1.0}};
1407+
1408+ // Declare (auxiliary) derivative matrix (of polynomial basis).
1409+ double dmats_old[3][3] = \
1410+ {{1.0, 0.0, 0.0},
1411+ {0.0, 1.0, 0.0},
1412+ {0.0, 0.0, 1.0}};
1413+
1414+ // Loop possible derivatives.
1415+ for (unsigned int r = 0; r < num_derivatives; r++)
1416+ {
1417+ // Resetting dmats values to compute next derivative.
1418+ for (unsigned int t = 0; t < 3; t++)
1419+ {
1420+ for (unsigned int u = 0; u < 3; u++)
1421+ {
1422+ dmats[t][u] = 0.0;
1423+ if (t == u)
1424+ {
1425+ dmats[t][u] = 1.0;
1426+ }
1427+
1428+ }// end loop over 'u'
1429+ }// end loop over 't'
1430+
1431+ // Looping derivative order to generate dmats.
1432+ for (unsigned int s = 0; s < n; s++)
1433+ {
1434+ // Updating dmats_old with new values and resetting dmats.
1435+ for (unsigned int t = 0; t < 3; t++)
1436+ {
1437+ for (unsigned int u = 0; u < 3; u++)
1438+ {
1439+ dmats_old[t][u] = dmats[t][u];
1440+ dmats[t][u] = 0.0;
1441+ }// end loop over 'u'
1442+ }// end loop over 't'
1443+
1444+ // Update dmats using an inner product.
1445+ if (combinations[r][s] == 0)
1446+ {
1447+ for (unsigned int t = 0; t < 3; t++)
1448+ {
1449+ for (unsigned int u = 0; u < 3; u++)
1450+ {
1451+ for (unsigned int tu = 0; tu < 3; tu++)
1452+ {
1453+ dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
1454+ }// end loop over 'tu'
1455+ }// end loop over 'u'
1456+ }// end loop over 't'
1457+ }
1458+
1459+ if (combinations[r][s] == 1)
1460+ {
1461+ for (unsigned int t = 0; t < 3; t++)
1462+ {
1463+ for (unsigned int u = 0; u < 3; u++)
1464+ {
1465+ for (unsigned int tu = 0; tu < 3; tu++)
1466+ {
1467+ dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
1468+ }// end loop over 'tu'
1469+ }// end loop over 'u'
1470+ }// end loop over 't'
1471+ }
1472+
1473+ }// end loop over 's'
1474+ for (unsigned int s = 0; s < 3; s++)
1475+ {
1476+ for (unsigned int t = 0; t < 3; t++)
1477+ {
1478+ derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
1479+ }// end loop over 't'
1480+ }// end loop over 's'
1481+ }// end loop over 'r'
1482+
1483+ // Transform derivatives back to physical element
1484+ for (unsigned int r = 0; r < num_derivatives; r++)
1485+ {
1486+ for (unsigned int s = 0; s < num_derivatives; s++)
1487+ {
1488+ values[r] += transform[r][s]*derivatives[s];
1489+ }// end loop over 's'
1490+ }// end loop over 'r'
1491+
1492+ // Delete pointer to array of derivatives on FIAT element
1493+ delete [] derivatives;
1494+
1495+ // Delete pointer to array of combinations of derivatives and transform
1496+ for (unsigned int r = 0; r < num_derivatives; r++)
1497+ {
1498+ delete [] combinations[r];
1499+ }// end loop over 'r'
1500+ delete [] combinations;
1501+ for (unsigned int r = 0; r < num_derivatives; r++)
1502+ {
1503+ delete [] transform[r];
1504+ }// end loop over 'r'
1505+ delete [] transform;
1506+ break;
1507+ }
1508+ }
1509+
1510+ }
1511+
1512+ /// Evaluate order n derivatives of all basis functions at given point x in cell
1513+ virtual void evaluate_basis_derivatives_all(std::size_t n,
1514+ double* values,
1515+ const double* x,
1516+ const double* vertex_coordinates,
1517+ int cell_orientation) const
1518+ {
1519+ // Compute number of derivatives.
1520+ unsigned int num_derivatives = 1;
1521+ for (unsigned int r = 0; r < n; r++)
1522+ {
1523+ num_derivatives *= 2;
1524+ }// end loop over 'r'
1525+
1526+ // Helper variable to hold values of a single dof.
1527+ double *dof_values = new double[num_derivatives];
1528+ for (unsigned int r = 0; r < num_derivatives; r++)
1529+ {
1530+ dof_values[r] = 0.0;
1531+ }// end loop over 'r'
1532+
1533+ // Loop dofs and call evaluate_basis_derivatives.
1534+ for (unsigned int r = 0; r < 3; r++)
1535+ {
1536+ evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
1537+ for (unsigned int s = 0; s < num_derivatives; s++)
1538+ {
1539+ values[r*num_derivatives + s] = dof_values[s];
1540+ }// end loop over 's'
1541+ }// end loop over 'r'
1542+
1543+ // Delete pointer.
1544+ delete [] dof_values;
1545+ }
1546+
1547+ /// Evaluate linear functional for dof i on the function f
1548+ virtual double evaluate_dof(std::size_t i,
1549+ const ufc::function& f,
1550+ const double* vertex_coordinates,
1551+ int cell_orientation,
1552+ const ufc::cell& c) const
1553+ {
1554+ // Declare variables for result of evaluation
1555+ double vals[1];
1556+
1557+ // Declare variable for physical coordinates
1558+ double y[2];
1559+ switch (i)
1560+ {
1561+ case 0:
1562+ {
1563+ y[0] = vertex_coordinates[0];
1564+ y[1] = vertex_coordinates[1];
1565+ f.evaluate(vals, y, c);
1566+ return vals[0];
1567+ break;
1568+ }
1569+ case 1:
1570+ {
1571+ y[0] = vertex_coordinates[2];
1572+ y[1] = vertex_coordinates[3];
1573+ f.evaluate(vals, y, c);
1574+ return vals[0];
1575+ break;
1576+ }
1577+ case 2:
1578+ {
1579+ y[0] = vertex_coordinates[4];
1580+ y[1] = vertex_coordinates[5];
1581+ f.evaluate(vals, y, c);
1582+ return vals[0];
1583+ break;
1584+ }
1585+ }
1586+
1587+ return 0.0;
1588+ }
1589+
1590+ /// Evaluate linear functionals for all dofs on the function f
1591+ virtual void evaluate_dofs(double* values,
1592+ const ufc::function& f,
1593+ const double* vertex_coordinates,
1594+ int cell_orientation,
1595+ const ufc::cell& c) const
1596+ {
1597+ // Declare variables for result of evaluation
1598+ double vals[1];
1599+
1600+ // Declare variable for physical coordinates
1601+ double y[2];
1602+ y[0] = vertex_coordinates[0];
1603+ y[1] = vertex_coordinates[1];
1604+ f.evaluate(vals, y, c);
1605+ values[0] = vals[0];
1606+ y[0] = vertex_coordinates[2];
1607+ y[1] = vertex_coordinates[3];
1608+ f.evaluate(vals, y, c);
1609+ values[1] = vals[0];
1610+ y[0] = vertex_coordinates[4];
1611+ y[1] = vertex_coordinates[5];
1612+ f.evaluate(vals, y, c);
1613+ values[2] = vals[0];
1614+ }
1615+
1616+ /// Interpolate vertex values from dof values
1617+ virtual void interpolate_vertex_values(double* vertex_values,
1618+ const double* dof_values,
1619+ const double* vertex_coordinates,
1620+ int cell_orientation,
1621+ const ufc::cell& c) const
1622+ {
1623+ // Evaluate function and change variables
1624+ vertex_values[0] = dof_values[0];
1625+ vertex_values[1] = dof_values[1];
1626+ vertex_values[2] = dof_values[2];
1627+ }
1628+
1629+ /// Map coordinate xhat from reference cell to coordinate x in cell
1630+ virtual void map_from_reference_cell(double* x,
1631+ const double* xhat,
1632+ const ufc::cell& c) const
1633+ {
1634+ throw std::runtime_error("map_from_reference_cell not yet implemented.");
1635+ }
1636+
1637+ /// Map from coordinate x in cell to coordinate xhat in reference cell
1638+ virtual void map_to_reference_cell(double* xhat,
1639+ const double* x,
1640+ const ufc::cell& c) const
1641+ {
1642+ throw std::runtime_error("map_to_reference_cell not yet implemented.");
1643+ }
1644+
1645+ /// Return the number of sub elements (for a mixed element)
1646+ virtual std::size_t num_sub_elements() const
1647+ {
1648+ return 0;
1649+ }
1650+
1651+ /// Create a new finite element for sub element i (for a mixed element)
1652+ virtual ufc::finite_element* create_sub_element(std::size_t i) const
1653+ {
1654+ return 0;
1655+ }
1656+
1657+ /// Create a new class instance
1658+ virtual ufc::finite_element* create() const
1659+ {
1660+ return new elasticity_finite_element_1();
1661+ }
1662+
1663+};
1664+
1665+/// This class defines the interface for a finite element.
1666+
1667+class elasticity_finite_element_2: public ufc::finite_element
1668+{
1669+public:
1670+
1671+ /// Constructor
1672+ elasticity_finite_element_2() : ufc::finite_element()
1673+ {
1674+ // Do nothing
1675+ }
1676+
1677+ /// Destructor
1678+ virtual ~elasticity_finite_element_2()
1679+ {
1680+ // Do nothing
1681+ }
1682+
1683+ /// Return a string identifying the finite element
1684+ virtual const char* signature() const
1685+ {
1686+ return "VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, 2, None)";
1687+ }
1688+
1689+ /// Return the cell shape
1690+ virtual ufc::shape cell_shape() const
1691+ {
1692+ return ufc::triangle;
1693+ }
1694+
1695+ /// Return the topological dimension of the cell shape
1696+ virtual std::size_t topological_dimension() const
1697+ {
1698+ return 2;
1699+ }
1700+
1701+ /// Return the geometric dimension of the cell shape
1702+ virtual std::size_t geometric_dimension() const
1703+ {
1704+ return 2;
1705+ }
1706+
1707+ /// Return the dimension of the finite element function space
1708+ virtual std::size_t space_dimension() const
1709+ {
1710+ return 6;
1711+ }
1712+
1713+ /// Return the rank of the value space
1714+ virtual std::size_t value_rank() const
1715+ {
1716+ return 1;
1717+ }
1718+
1719+ /// Return the dimension of the value space for axis i
1720+ virtual std::size_t value_dimension(std::size_t i) const
1721+ {
1722+ switch (i)
1723+ {
1724+ case 0:
1725+ {
1726+ return 2;
1727+ break;
1728+ }
1729+ }
1730+
1731+ return 0;
1732+ }
1733+
1734+ /// Evaluate basis function i at given point x in cell
1735+ virtual void evaluate_basis(std::size_t i,
1736+ double* values,
1737+ const double* x,
1738+ const double* vertex_coordinates,
1739+ int cell_orientation) const
1740+ {
1741+ // Compute Jacobian
1742+ double J[4];
1743+ compute_jacobian_triangle_2d(J, vertex_coordinates);
1744+
1745+ // Compute Jacobian inverse and determinant
1746+ double K[4];
1747+ double detJ;
1748+ compute_jacobian_inverse_triangle_2d(K, detJ, J);
1749+
1750+
1751+ // Compute constants
1752+ const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
1753+ const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
1754+
1755+ // Get coordinates and map to the reference (FIAT) element
1756+ double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
1757+ double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1758+
1759+ // Reset values
1760+ values[0] = 0.0;
1761+ values[1] = 0.0;
1762+ switch (i)
1763+ {
1764+ case 0:
1765+ {
1766+
1767+ // Array of basisvalues
1768+ double basisvalues[3] = {0.0, 0.0, 0.0};
1769+
1770+ // Declare helper variables
1771+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1772+
1773+ // Compute basisvalues
1774+ basisvalues[0] = 1.0;
1775+ basisvalues[1] = tmp0;
1776+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1777+ basisvalues[0] *= std::sqrt(0.5);
1778+ basisvalues[2] *= std::sqrt(1.0);
1779+ basisvalues[1] *= std::sqrt(3.0);
1780+
1781+ // Table(s) of coefficients
1782+ static const double coefficients0[3] = \
1783+ {0.471404520791032, -0.288675134594813, -0.166666666666667};
1784+
1785+ // Compute value(s)
1786+ for (unsigned int r = 0; r < 3; r++)
1787+ {
1788+ values[0] += coefficients0[r]*basisvalues[r];
1789+ }// end loop over 'r'
1790+ break;
1791+ }
1792+ case 1:
1793+ {
1794+
1795+ // Array of basisvalues
1796+ double basisvalues[3] = {0.0, 0.0, 0.0};
1797+
1798+ // Declare helper variables
1799+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1800+
1801+ // Compute basisvalues
1802+ basisvalues[0] = 1.0;
1803+ basisvalues[1] = tmp0;
1804+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1805+ basisvalues[0] *= std::sqrt(0.5);
1806+ basisvalues[2] *= std::sqrt(1.0);
1807+ basisvalues[1] *= std::sqrt(3.0);
1808+
1809+ // Table(s) of coefficients
1810+ static const double coefficients0[3] = \
1811+ {0.471404520791032, 0.288675134594813, -0.166666666666667};
1812+
1813+ // Compute value(s)
1814+ for (unsigned int r = 0; r < 3; r++)
1815+ {
1816+ values[0] += coefficients0[r]*basisvalues[r];
1817+ }// end loop over 'r'
1818+ break;
1819+ }
1820+ case 2:
1821+ {
1822+
1823+ // Array of basisvalues
1824+ double basisvalues[3] = {0.0, 0.0, 0.0};
1825+
1826+ // Declare helper variables
1827+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1828+
1829+ // Compute basisvalues
1830+ basisvalues[0] = 1.0;
1831+ basisvalues[1] = tmp0;
1832+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1833+ basisvalues[0] *= std::sqrt(0.5);
1834+ basisvalues[2] *= std::sqrt(1.0);
1835+ basisvalues[1] *= std::sqrt(3.0);
1836+
1837+ // Table(s) of coefficients
1838+ static const double coefficients0[3] = \
1839+ {0.471404520791032, 0.0, 0.333333333333333};
1840+
1841+ // Compute value(s)
1842+ for (unsigned int r = 0; r < 3; r++)
1843+ {
1844+ values[0] += coefficients0[r]*basisvalues[r];
1845+ }// end loop over 'r'
1846+ break;
1847+ }
1848+ case 3:
1849+ {
1850+
1851+ // Array of basisvalues
1852+ double basisvalues[3] = {0.0, 0.0, 0.0};
1853+
1854+ // Declare helper variables
1855+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1856+
1857+ // Compute basisvalues
1858+ basisvalues[0] = 1.0;
1859+ basisvalues[1] = tmp0;
1860+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1861+ basisvalues[0] *= std::sqrt(0.5);
1862+ basisvalues[2] *= std::sqrt(1.0);
1863+ basisvalues[1] *= std::sqrt(3.0);
1864+
1865+ // Table(s) of coefficients
1866+ static const double coefficients0[3] = \
1867+ {0.471404520791032, -0.288675134594813, -0.166666666666667};
1868+
1869+ // Compute value(s)
1870+ for (unsigned int r = 0; r < 3; r++)
1871+ {
1872+ values[1] += coefficients0[r]*basisvalues[r];
1873+ }// end loop over 'r'
1874+ break;
1875+ }
1876+ case 4:
1877+ {
1878+
1879+ // Array of basisvalues
1880+ double basisvalues[3] = {0.0, 0.0, 0.0};
1881+
1882+ // Declare helper variables
1883+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1884+
1885+ // Compute basisvalues
1886+ basisvalues[0] = 1.0;
1887+ basisvalues[1] = tmp0;
1888+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1889+ basisvalues[0] *= std::sqrt(0.5);
1890+ basisvalues[2] *= std::sqrt(1.0);
1891+ basisvalues[1] *= std::sqrt(3.0);
1892+
1893+ // Table(s) of coefficients
1894+ static const double coefficients0[3] = \
1895+ {0.471404520791032, 0.288675134594813, -0.166666666666667};
1896+
1897+ // Compute value(s)
1898+ for (unsigned int r = 0; r < 3; r++)
1899+ {
1900+ values[1] += coefficients0[r]*basisvalues[r];
1901+ }// end loop over 'r'
1902+ break;
1903+ }
1904+ case 5:
1905+ {
1906+
1907+ // Array of basisvalues
1908+ double basisvalues[3] = {0.0, 0.0, 0.0};
1909+
1910+ // Declare helper variables
1911+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
1912+
1913+ // Compute basisvalues
1914+ basisvalues[0] = 1.0;
1915+ basisvalues[1] = tmp0;
1916+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
1917+ basisvalues[0] *= std::sqrt(0.5);
1918+ basisvalues[2] *= std::sqrt(1.0);
1919+ basisvalues[1] *= std::sqrt(3.0);
1920+
1921+ // Table(s) of coefficients
1922+ static const double coefficients0[3] = \
1923+ {0.471404520791032, 0.0, 0.333333333333333};
1924+
1925+ // Compute value(s)
1926+ for (unsigned int r = 0; r < 3; r++)
1927+ {
1928+ values[1] += coefficients0[r]*basisvalues[r];
1929+ }// end loop over 'r'
1930+ break;
1931+ }
1932+ }
1933+
1934+ }
1935+
1936+ /// Evaluate all basis functions at given point x in cell
1937+ virtual void evaluate_basis_all(double* values,
1938+ const double* x,
1939+ const double* vertex_coordinates,
1940+ int cell_orientation) const
1941+ {
1942+ // Helper variable to hold values of a single dof.
1943+ double dof_values[2] = {0.0, 0.0};
1944+
1945+ // Loop dofs and call evaluate_basis
1946+ for (unsigned int r = 0; r < 6; r++)
1947+ {
1948+ evaluate_basis(r, dof_values, x, vertex_coordinates, cell_orientation);
1949+ for (unsigned int s = 0; s < 2; s++)
1950+ {
1951+ values[r*2 + s] = dof_values[s];
1952+ }// end loop over 's'
1953+ }// end loop over 'r'
1954+ }
1955+
1956+ /// Evaluate order n derivatives of basis function i at given point x in cell
1957+ virtual void evaluate_basis_derivatives(std::size_t i,
1958+ std::size_t n,
1959+ double* values,
1960+ const double* x,
1961+ const double* vertex_coordinates,
1962+ int cell_orientation) const
1963+ {
1964+ // Compute Jacobian
1965+ double J[4];
1966+ compute_jacobian_triangle_2d(J, vertex_coordinates);
1967+
1968+ // Compute Jacobian inverse and determinant
1969+ double K[4];
1970+ double detJ;
1971+ compute_jacobian_inverse_triangle_2d(K, detJ, J);
1972+
1973+
1974+ // Compute constants
1975+ const double C0 = vertex_coordinates[2] + vertex_coordinates[4];
1976+ const double C1 = vertex_coordinates[3] + vertex_coordinates[5];
1977+
1978+ // Get coordinates and map to the reference (FIAT) element
1979+ double X = (J[1]*(C1 - 2.0*x[1]) + J[3]*(2.0*x[0] - C0)) / detJ;
1980+ double Y = (J[0]*(2.0*x[1] - C1) + J[2]*(C0 - 2.0*x[0])) / detJ;
1981+
1982+ // Compute number of derivatives.
1983+ unsigned int num_derivatives = 1;
1984+ for (unsigned int r = 0; r < n; r++)
1985+ {
1986+ num_derivatives *= 2;
1987+ }// end loop over 'r'
1988+
1989+ // Declare pointer to two dimensional array that holds combinations of derivatives and initialise
1990+ unsigned int **combinations = new unsigned int *[num_derivatives];
1991+ for (unsigned int row = 0; row < num_derivatives; row++)
1992+ {
1993+ combinations[row] = new unsigned int [n];
1994+ for (unsigned int col = 0; col < n; col++)
1995+ combinations[row][col] = 0;
1996+ }
1997+
1998+ // Generate combinations of derivatives
1999+ for (unsigned int row = 1; row < num_derivatives; row++)
2000+ {
2001+ for (unsigned int num = 0; num < row; num++)
2002+ {
2003+ for (unsigned int col = n-1; col+1 > 0; col--)
2004+ {
2005+ if (combinations[row][col] + 1 > 1)
2006+ combinations[row][col] = 0;
2007+ else
2008+ {
2009+ combinations[row][col] += 1;
2010+ break;
2011+ }
2012+ }
2013+ }
2014+ }
2015+
2016+ // Compute inverse of Jacobian
2017+ const double Jinv[2][2] = {{K[0], K[1]}, {K[2], K[3]}};
2018+
2019+ // Declare transformation matrix
2020+ // Declare pointer to two dimensional array and initialise
2021+ double **transform = new double *[num_derivatives];
2022+
2023+ for (unsigned int j = 0; j < num_derivatives; j++)
2024+ {
2025+ transform[j] = new double [num_derivatives];
2026+ for (unsigned int k = 0; k < num_derivatives; k++)
2027+ transform[j][k] = 1;
2028+ }
2029+
2030+ // Construct transformation matrix
2031+ for (unsigned int row = 0; row < num_derivatives; row++)
2032+ {
2033+ for (unsigned int col = 0; col < num_derivatives; col++)
2034+ {
2035+ for (unsigned int k = 0; k < n; k++)
2036+ transform[row][col] *= Jinv[combinations[col][k]][combinations[row][k]];
2037+ }
2038+ }
2039+
2040+ // Reset values. Assuming that values is always an array.
2041+ for (unsigned int r = 0; r < 2*num_derivatives; r++)
2042+ {
2043+ values[r] = 0.0;
2044+ }// end loop over 'r'
2045+
2046+ switch (i)
2047+ {
2048+ case 0:
2049+ {
2050+
2051+ // Array of basisvalues
2052+ double basisvalues[3] = {0.0, 0.0, 0.0};
2053+
2054+ // Declare helper variables
2055+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2056+
2057+ // Compute basisvalues
2058+ basisvalues[0] = 1.0;
2059+ basisvalues[1] = tmp0;
2060+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2061+ basisvalues[0] *= std::sqrt(0.5);
2062+ basisvalues[2] *= std::sqrt(1.0);
2063+ basisvalues[1] *= std::sqrt(3.0);
2064+
2065+ // Table(s) of coefficients
2066+ static const double coefficients0[3] = \
2067+ {0.471404520791032, -0.288675134594813, -0.166666666666667};
2068+
2069+ // Tables of derivatives of the polynomial base (transpose).
2070+ static const double dmats0[3][3] = \
2071+ {{0.0, 0.0, 0.0},
2072+ {4.89897948556636, 0.0, 0.0},
2073+ {0.0, 0.0, 0.0}};
2074+
2075+ static const double dmats1[3][3] = \
2076+ {{0.0, 0.0, 0.0},
2077+ {2.44948974278318, 0.0, 0.0},
2078+ {4.24264068711928, 0.0, 0.0}};
2079+
2080+ // Compute reference derivatives.
2081+ // Declare pointer to array of derivatives on FIAT element.
2082+ double *derivatives = new double[num_derivatives];
2083+ for (unsigned int r = 0; r < num_derivatives; r++)
2084+ {
2085+ derivatives[r] = 0.0;
2086+ }// end loop over 'r'
2087+
2088+ // Declare derivative matrix (of polynomial basis).
2089+ double dmats[3][3] = \
2090+ {{1.0, 0.0, 0.0},
2091+ {0.0, 1.0, 0.0},
2092+ {0.0, 0.0, 1.0}};
2093+
2094+ // Declare (auxiliary) derivative matrix (of polynomial basis).
2095+ double dmats_old[3][3] = \
2096+ {{1.0, 0.0, 0.0},
2097+ {0.0, 1.0, 0.0},
2098+ {0.0, 0.0, 1.0}};
2099+
2100+ // Loop possible derivatives.
2101+ for (unsigned int r = 0; r < num_derivatives; r++)
2102+ {
2103+ // Resetting dmats values to compute next derivative.
2104+ for (unsigned int t = 0; t < 3; t++)
2105+ {
2106+ for (unsigned int u = 0; u < 3; u++)
2107+ {
2108+ dmats[t][u] = 0.0;
2109+ if (t == u)
2110+ {
2111+ dmats[t][u] = 1.0;
2112+ }
2113+
2114+ }// end loop over 'u'
2115+ }// end loop over 't'
2116+
2117+ // Looping derivative order to generate dmats.
2118+ for (unsigned int s = 0; s < n; s++)
2119+ {
2120+ // Updating dmats_old with new values and resetting dmats.
2121+ for (unsigned int t = 0; t < 3; t++)
2122+ {
2123+ for (unsigned int u = 0; u < 3; u++)
2124+ {
2125+ dmats_old[t][u] = dmats[t][u];
2126+ dmats[t][u] = 0.0;
2127+ }// end loop over 'u'
2128+ }// end loop over 't'
2129+
2130+ // Update dmats using an inner product.
2131+ if (combinations[r][s] == 0)
2132+ {
2133+ for (unsigned int t = 0; t < 3; t++)
2134+ {
2135+ for (unsigned int u = 0; u < 3; u++)
2136+ {
2137+ for (unsigned int tu = 0; tu < 3; tu++)
2138+ {
2139+ dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2140+ }// end loop over 'tu'
2141+ }// end loop over 'u'
2142+ }// end loop over 't'
2143+ }
2144+
2145+ if (combinations[r][s] == 1)
2146+ {
2147+ for (unsigned int t = 0; t < 3; t++)
2148+ {
2149+ for (unsigned int u = 0; u < 3; u++)
2150+ {
2151+ for (unsigned int tu = 0; tu < 3; tu++)
2152+ {
2153+ dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2154+ }// end loop over 'tu'
2155+ }// end loop over 'u'
2156+ }// end loop over 't'
2157+ }
2158+
2159+ }// end loop over 's'
2160+ for (unsigned int s = 0; s < 3; s++)
2161+ {
2162+ for (unsigned int t = 0; t < 3; t++)
2163+ {
2164+ derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2165+ }// end loop over 't'
2166+ }// end loop over 's'
2167+ }// end loop over 'r'
2168+
2169+ // Transform derivatives back to physical element
2170+ for (unsigned int r = 0; r < num_derivatives; r++)
2171+ {
2172+ for (unsigned int s = 0; s < num_derivatives; s++)
2173+ {
2174+ values[r] += transform[r][s]*derivatives[s];
2175+ }// end loop over 's'
2176+ }// end loop over 'r'
2177+
2178+ // Delete pointer to array of derivatives on FIAT element
2179+ delete [] derivatives;
2180+
2181+ // Delete pointer to array of combinations of derivatives and transform
2182+ for (unsigned int r = 0; r < num_derivatives; r++)
2183+ {
2184+ delete [] combinations[r];
2185+ }// end loop over 'r'
2186+ delete [] combinations;
2187+ for (unsigned int r = 0; r < num_derivatives; r++)
2188+ {
2189+ delete [] transform[r];
2190+ }// end loop over 'r'
2191+ delete [] transform;
2192+ break;
2193+ }
2194+ case 1:
2195+ {
2196+
2197+ // Array of basisvalues
2198+ double basisvalues[3] = {0.0, 0.0, 0.0};
2199+
2200+ // Declare helper variables
2201+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2202+
2203+ // Compute basisvalues
2204+ basisvalues[0] = 1.0;
2205+ basisvalues[1] = tmp0;
2206+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2207+ basisvalues[0] *= std::sqrt(0.5);
2208+ basisvalues[2] *= std::sqrt(1.0);
2209+ basisvalues[1] *= std::sqrt(3.0);
2210+
2211+ // Table(s) of coefficients
2212+ static const double coefficients0[3] = \
2213+ {0.471404520791032, 0.288675134594813, -0.166666666666667};
2214+
2215+ // Tables of derivatives of the polynomial base (transpose).
2216+ static const double dmats0[3][3] = \
2217+ {{0.0, 0.0, 0.0},
2218+ {4.89897948556636, 0.0, 0.0},
2219+ {0.0, 0.0, 0.0}};
2220+
2221+ static const double dmats1[3][3] = \
2222+ {{0.0, 0.0, 0.0},
2223+ {2.44948974278318, 0.0, 0.0},
2224+ {4.24264068711928, 0.0, 0.0}};
2225+
2226+ // Compute reference derivatives.
2227+ // Declare pointer to array of derivatives on FIAT element.
2228+ double *derivatives = new double[num_derivatives];
2229+ for (unsigned int r = 0; r < num_derivatives; r++)
2230+ {
2231+ derivatives[r] = 0.0;
2232+ }// end loop over 'r'
2233+
2234+ // Declare derivative matrix (of polynomial basis).
2235+ double dmats[3][3] = \
2236+ {{1.0, 0.0, 0.0},
2237+ {0.0, 1.0, 0.0},
2238+ {0.0, 0.0, 1.0}};
2239+
2240+ // Declare (auxiliary) derivative matrix (of polynomial basis).
2241+ double dmats_old[3][3] = \
2242+ {{1.0, 0.0, 0.0},
2243+ {0.0, 1.0, 0.0},
2244+ {0.0, 0.0, 1.0}};
2245+
2246+ // Loop possible derivatives.
2247+ for (unsigned int r = 0; r < num_derivatives; r++)
2248+ {
2249+ // Resetting dmats values to compute next derivative.
2250+ for (unsigned int t = 0; t < 3; t++)
2251+ {
2252+ for (unsigned int u = 0; u < 3; u++)
2253+ {
2254+ dmats[t][u] = 0.0;
2255+ if (t == u)
2256+ {
2257+ dmats[t][u] = 1.0;
2258+ }
2259+
2260+ }// end loop over 'u'
2261+ }// end loop over 't'
2262+
2263+ // Looping derivative order to generate dmats.
2264+ for (unsigned int s = 0; s < n; s++)
2265+ {
2266+ // Updating dmats_old with new values and resetting dmats.
2267+ for (unsigned int t = 0; t < 3; t++)
2268+ {
2269+ for (unsigned int u = 0; u < 3; u++)
2270+ {
2271+ dmats_old[t][u] = dmats[t][u];
2272+ dmats[t][u] = 0.0;
2273+ }// end loop over 'u'
2274+ }// end loop over 't'
2275+
2276+ // Update dmats using an inner product.
2277+ if (combinations[r][s] == 0)
2278+ {
2279+ for (unsigned int t = 0; t < 3; t++)
2280+ {
2281+ for (unsigned int u = 0; u < 3; u++)
2282+ {
2283+ for (unsigned int tu = 0; tu < 3; tu++)
2284+ {
2285+ dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2286+ }// end loop over 'tu'
2287+ }// end loop over 'u'
2288+ }// end loop over 't'
2289+ }
2290+
2291+ if (combinations[r][s] == 1)
2292+ {
2293+ for (unsigned int t = 0; t < 3; t++)
2294+ {
2295+ for (unsigned int u = 0; u < 3; u++)
2296+ {
2297+ for (unsigned int tu = 0; tu < 3; tu++)
2298+ {
2299+ dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2300+ }// end loop over 'tu'
2301+ }// end loop over 'u'
2302+ }// end loop over 't'
2303+ }
2304+
2305+ }// end loop over 's'
2306+ for (unsigned int s = 0; s < 3; s++)
2307+ {
2308+ for (unsigned int t = 0; t < 3; t++)
2309+ {
2310+ derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2311+ }// end loop over 't'
2312+ }// end loop over 's'
2313+ }// end loop over 'r'
2314+
2315+ // Transform derivatives back to physical element
2316+ for (unsigned int r = 0; r < num_derivatives; r++)
2317+ {
2318+ for (unsigned int s = 0; s < num_derivatives; s++)
2319+ {
2320+ values[r] += transform[r][s]*derivatives[s];
2321+ }// end loop over 's'
2322+ }// end loop over 'r'
2323+
2324+ // Delete pointer to array of derivatives on FIAT element
2325+ delete [] derivatives;
2326+
2327+ // Delete pointer to array of combinations of derivatives and transform
2328+ for (unsigned int r = 0; r < num_derivatives; r++)
2329+ {
2330+ delete [] combinations[r];
2331+ }// end loop over 'r'
2332+ delete [] combinations;
2333+ for (unsigned int r = 0; r < num_derivatives; r++)
2334+ {
2335+ delete [] transform[r];
2336+ }// end loop over 'r'
2337+ delete [] transform;
2338+ break;
2339+ }
2340+ case 2:
2341+ {
2342+
2343+ // Array of basisvalues
2344+ double basisvalues[3] = {0.0, 0.0, 0.0};
2345+
2346+ // Declare helper variables
2347+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2348+
2349+ // Compute basisvalues
2350+ basisvalues[0] = 1.0;
2351+ basisvalues[1] = tmp0;
2352+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2353+ basisvalues[0] *= std::sqrt(0.5);
2354+ basisvalues[2] *= std::sqrt(1.0);
2355+ basisvalues[1] *= std::sqrt(3.0);
2356+
2357+ // Table(s) of coefficients
2358+ static const double coefficients0[3] = \
2359+ {0.471404520791032, 0.0, 0.333333333333333};
2360+
2361+ // Tables of derivatives of the polynomial base (transpose).
2362+ static const double dmats0[3][3] = \
2363+ {{0.0, 0.0, 0.0},
2364+ {4.89897948556636, 0.0, 0.0},
2365+ {0.0, 0.0, 0.0}};
2366+
2367+ static const double dmats1[3][3] = \
2368+ {{0.0, 0.0, 0.0},
2369+ {2.44948974278318, 0.0, 0.0},
2370+ {4.24264068711928, 0.0, 0.0}};
2371+
2372+ // Compute reference derivatives.
2373+ // Declare pointer to array of derivatives on FIAT element.
2374+ double *derivatives = new double[num_derivatives];
2375+ for (unsigned int r = 0; r < num_derivatives; r++)
2376+ {
2377+ derivatives[r] = 0.0;
2378+ }// end loop over 'r'
2379+
2380+ // Declare derivative matrix (of polynomial basis).
2381+ double dmats[3][3] = \
2382+ {{1.0, 0.0, 0.0},
2383+ {0.0, 1.0, 0.0},
2384+ {0.0, 0.0, 1.0}};
2385+
2386+ // Declare (auxiliary) derivative matrix (of polynomial basis).
2387+ double dmats_old[3][3] = \
2388+ {{1.0, 0.0, 0.0},
2389+ {0.0, 1.0, 0.0},
2390+ {0.0, 0.0, 1.0}};
2391+
2392+ // Loop possible derivatives.
2393+ for (unsigned int r = 0; r < num_derivatives; r++)
2394+ {
2395+ // Resetting dmats values to compute next derivative.
2396+ for (unsigned int t = 0; t < 3; t++)
2397+ {
2398+ for (unsigned int u = 0; u < 3; u++)
2399+ {
2400+ dmats[t][u] = 0.0;
2401+ if (t == u)
2402+ {
2403+ dmats[t][u] = 1.0;
2404+ }
2405+
2406+ }// end loop over 'u'
2407+ }// end loop over 't'
2408+
2409+ // Looping derivative order to generate dmats.
2410+ for (unsigned int s = 0; s < n; s++)
2411+ {
2412+ // Updating dmats_old with new values and resetting dmats.
2413+ for (unsigned int t = 0; t < 3; t++)
2414+ {
2415+ for (unsigned int u = 0; u < 3; u++)
2416+ {
2417+ dmats_old[t][u] = dmats[t][u];
2418+ dmats[t][u] = 0.0;
2419+ }// end loop over 'u'
2420+ }// end loop over 't'
2421+
2422+ // Update dmats using an inner product.
2423+ if (combinations[r][s] == 0)
2424+ {
2425+ for (unsigned int t = 0; t < 3; t++)
2426+ {
2427+ for (unsigned int u = 0; u < 3; u++)
2428+ {
2429+ for (unsigned int tu = 0; tu < 3; tu++)
2430+ {
2431+ dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2432+ }// end loop over 'tu'
2433+ }// end loop over 'u'
2434+ }// end loop over 't'
2435+ }
2436+
2437+ if (combinations[r][s] == 1)
2438+ {
2439+ for (unsigned int t = 0; t < 3; t++)
2440+ {
2441+ for (unsigned int u = 0; u < 3; u++)
2442+ {
2443+ for (unsigned int tu = 0; tu < 3; tu++)
2444+ {
2445+ dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2446+ }// end loop over 'tu'
2447+ }// end loop over 'u'
2448+ }// end loop over 't'
2449+ }
2450+
2451+ }// end loop over 's'
2452+ for (unsigned int s = 0; s < 3; s++)
2453+ {
2454+ for (unsigned int t = 0; t < 3; t++)
2455+ {
2456+ derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2457+ }// end loop over 't'
2458+ }// end loop over 's'
2459+ }// end loop over 'r'
2460+
2461+ // Transform derivatives back to physical element
2462+ for (unsigned int r = 0; r < num_derivatives; r++)
2463+ {
2464+ for (unsigned int s = 0; s < num_derivatives; s++)
2465+ {
2466+ values[r] += transform[r][s]*derivatives[s];
2467+ }// end loop over 's'
2468+ }// end loop over 'r'
2469+
2470+ // Delete pointer to array of derivatives on FIAT element
2471+ delete [] derivatives;
2472+
2473+ // Delete pointer to array of combinations of derivatives and transform
2474+ for (unsigned int r = 0; r < num_derivatives; r++)
2475+ {
2476+ delete [] combinations[r];
2477+ }// end loop over 'r'
2478+ delete [] combinations;
2479+ for (unsigned int r = 0; r < num_derivatives; r++)
2480+ {
2481+ delete [] transform[r];
2482+ }// end loop over 'r'
2483+ delete [] transform;
2484+ break;
2485+ }
2486+ case 3:
2487+ {
2488+
2489+ // Array of basisvalues
2490+ double basisvalues[3] = {0.0, 0.0, 0.0};
2491+
2492+ // Declare helper variables
2493+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2494+
2495+ // Compute basisvalues
2496+ basisvalues[0] = 1.0;
2497+ basisvalues[1] = tmp0;
2498+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2499+ basisvalues[0] *= std::sqrt(0.5);
2500+ basisvalues[2] *= std::sqrt(1.0);
2501+ basisvalues[1] *= std::sqrt(3.0);
2502+
2503+ // Table(s) of coefficients
2504+ static const double coefficients0[3] = \
2505+ {0.471404520791032, -0.288675134594813, -0.166666666666667};
2506+
2507+ // Tables of derivatives of the polynomial base (transpose).
2508+ static const double dmats0[3][3] = \
2509+ {{0.0, 0.0, 0.0},
2510+ {4.89897948556636, 0.0, 0.0},
2511+ {0.0, 0.0, 0.0}};
2512+
2513+ static const double dmats1[3][3] = \
2514+ {{0.0, 0.0, 0.0},
2515+ {2.44948974278318, 0.0, 0.0},
2516+ {4.24264068711928, 0.0, 0.0}};
2517+
2518+ // Compute reference derivatives.
2519+ // Declare pointer to array of derivatives on FIAT element.
2520+ double *derivatives = new double[num_derivatives];
2521+ for (unsigned int r = 0; r < num_derivatives; r++)
2522+ {
2523+ derivatives[r] = 0.0;
2524+ }// end loop over 'r'
2525+
2526+ // Declare derivative matrix (of polynomial basis).
2527+ double dmats[3][3] = \
2528+ {{1.0, 0.0, 0.0},
2529+ {0.0, 1.0, 0.0},
2530+ {0.0, 0.0, 1.0}};
2531+
2532+ // Declare (auxiliary) derivative matrix (of polynomial basis).
2533+ double dmats_old[3][3] = \
2534+ {{1.0, 0.0, 0.0},
2535+ {0.0, 1.0, 0.0},
2536+ {0.0, 0.0, 1.0}};
2537+
2538+ // Loop possible derivatives.
2539+ for (unsigned int r = 0; r < num_derivatives; r++)
2540+ {
2541+ // Resetting dmats values to compute next derivative.
2542+ for (unsigned int t = 0; t < 3; t++)
2543+ {
2544+ for (unsigned int u = 0; u < 3; u++)
2545+ {
2546+ dmats[t][u] = 0.0;
2547+ if (t == u)
2548+ {
2549+ dmats[t][u] = 1.0;
2550+ }
2551+
2552+ }// end loop over 'u'
2553+ }// end loop over 't'
2554+
2555+ // Looping derivative order to generate dmats.
2556+ for (unsigned int s = 0; s < n; s++)
2557+ {
2558+ // Updating dmats_old with new values and resetting dmats.
2559+ for (unsigned int t = 0; t < 3; t++)
2560+ {
2561+ for (unsigned int u = 0; u < 3; u++)
2562+ {
2563+ dmats_old[t][u] = dmats[t][u];
2564+ dmats[t][u] = 0.0;
2565+ }// end loop over 'u'
2566+ }// end loop over 't'
2567+
2568+ // Update dmats using an inner product.
2569+ if (combinations[r][s] == 0)
2570+ {
2571+ for (unsigned int t = 0; t < 3; t++)
2572+ {
2573+ for (unsigned int u = 0; u < 3; u++)
2574+ {
2575+ for (unsigned int tu = 0; tu < 3; tu++)
2576+ {
2577+ dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2578+ }// end loop over 'tu'
2579+ }// end loop over 'u'
2580+ }// end loop over 't'
2581+ }
2582+
2583+ if (combinations[r][s] == 1)
2584+ {
2585+ for (unsigned int t = 0; t < 3; t++)
2586+ {
2587+ for (unsigned int u = 0; u < 3; u++)
2588+ {
2589+ for (unsigned int tu = 0; tu < 3; tu++)
2590+ {
2591+ dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2592+ }// end loop over 'tu'
2593+ }// end loop over 'u'
2594+ }// end loop over 't'
2595+ }
2596+
2597+ }// end loop over 's'
2598+ for (unsigned int s = 0; s < 3; s++)
2599+ {
2600+ for (unsigned int t = 0; t < 3; t++)
2601+ {
2602+ derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2603+ }// end loop over 't'
2604+ }// end loop over 's'
2605+ }// end loop over 'r'
2606+
2607+ // Transform derivatives back to physical element
2608+ for (unsigned int r = 0; r < num_derivatives; r++)
2609+ {
2610+ for (unsigned int s = 0; s < num_derivatives; s++)
2611+ {
2612+ values[num_derivatives + r] += transform[r][s]*derivatives[s];
2613+ }// end loop over 's'
2614+ }// end loop over 'r'
2615+
2616+ // Delete pointer to array of derivatives on FIAT element
2617+ delete [] derivatives;
2618+
2619+ // Delete pointer to array of combinations of derivatives and transform
2620+ for (unsigned int r = 0; r < num_derivatives; r++)
2621+ {
2622+ delete [] combinations[r];
2623+ }// end loop over 'r'
2624+ delete [] combinations;
2625+ for (unsigned int r = 0; r < num_derivatives; r++)
2626+ {
2627+ delete [] transform[r];
2628+ }// end loop over 'r'
2629+ delete [] transform;
2630+ break;
2631+ }
2632+ case 4:
2633+ {
2634+
2635+ // Array of basisvalues
2636+ double basisvalues[3] = {0.0, 0.0, 0.0};
2637+
2638+ // Declare helper variables
2639+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2640+
2641+ // Compute basisvalues
2642+ basisvalues[0] = 1.0;
2643+ basisvalues[1] = tmp0;
2644+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2645+ basisvalues[0] *= std::sqrt(0.5);
2646+ basisvalues[2] *= std::sqrt(1.0);
2647+ basisvalues[1] *= std::sqrt(3.0);
2648+
2649+ // Table(s) of coefficients
2650+ static const double coefficients0[3] = \
2651+ {0.471404520791032, 0.288675134594813, -0.166666666666667};
2652+
2653+ // Tables of derivatives of the polynomial base (transpose).
2654+ static const double dmats0[3][3] = \
2655+ {{0.0, 0.0, 0.0},
2656+ {4.89897948556636, 0.0, 0.0},
2657+ {0.0, 0.0, 0.0}};
2658+
2659+ static const double dmats1[3][3] = \
2660+ {{0.0, 0.0, 0.0},
2661+ {2.44948974278318, 0.0, 0.0},
2662+ {4.24264068711928, 0.0, 0.0}};
2663+
2664+ // Compute reference derivatives.
2665+ // Declare pointer to array of derivatives on FIAT element.
2666+ double *derivatives = new double[num_derivatives];
2667+ for (unsigned int r = 0; r < num_derivatives; r++)
2668+ {
2669+ derivatives[r] = 0.0;
2670+ }// end loop over 'r'
2671+
2672+ // Declare derivative matrix (of polynomial basis).
2673+ double dmats[3][3] = \
2674+ {{1.0, 0.0, 0.0},
2675+ {0.0, 1.0, 0.0},
2676+ {0.0, 0.0, 1.0}};
2677+
2678+ // Declare (auxiliary) derivative matrix (of polynomial basis).
2679+ double dmats_old[3][3] = \
2680+ {{1.0, 0.0, 0.0},
2681+ {0.0, 1.0, 0.0},
2682+ {0.0, 0.0, 1.0}};
2683+
2684+ // Loop possible derivatives.
2685+ for (unsigned int r = 0; r < num_derivatives; r++)
2686+ {
2687+ // Resetting dmats values to compute next derivative.
2688+ for (unsigned int t = 0; t < 3; t++)
2689+ {
2690+ for (unsigned int u = 0; u < 3; u++)
2691+ {
2692+ dmats[t][u] = 0.0;
2693+ if (t == u)
2694+ {
2695+ dmats[t][u] = 1.0;
2696+ }
2697+
2698+ }// end loop over 'u'
2699+ }// end loop over 't'
2700+
2701+ // Looping derivative order to generate dmats.
2702+ for (unsigned int s = 0; s < n; s++)
2703+ {
2704+ // Updating dmats_old with new values and resetting dmats.
2705+ for (unsigned int t = 0; t < 3; t++)
2706+ {
2707+ for (unsigned int u = 0; u < 3; u++)
2708+ {
2709+ dmats_old[t][u] = dmats[t][u];
2710+ dmats[t][u] = 0.0;
2711+ }// end loop over 'u'
2712+ }// end loop over 't'
2713+
2714+ // Update dmats using an inner product.
2715+ if (combinations[r][s] == 0)
2716+ {
2717+ for (unsigned int t = 0; t < 3; t++)
2718+ {
2719+ for (unsigned int u = 0; u < 3; u++)
2720+ {
2721+ for (unsigned int tu = 0; tu < 3; tu++)
2722+ {
2723+ dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2724+ }// end loop over 'tu'
2725+ }// end loop over 'u'
2726+ }// end loop over 't'
2727+ }
2728+
2729+ if (combinations[r][s] == 1)
2730+ {
2731+ for (unsigned int t = 0; t < 3; t++)
2732+ {
2733+ for (unsigned int u = 0; u < 3; u++)
2734+ {
2735+ for (unsigned int tu = 0; tu < 3; tu++)
2736+ {
2737+ dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2738+ }// end loop over 'tu'
2739+ }// end loop over 'u'
2740+ }// end loop over 't'
2741+ }
2742+
2743+ }// end loop over 's'
2744+ for (unsigned int s = 0; s < 3; s++)
2745+ {
2746+ for (unsigned int t = 0; t < 3; t++)
2747+ {
2748+ derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2749+ }// end loop over 't'
2750+ }// end loop over 's'
2751+ }// end loop over 'r'
2752+
2753+ // Transform derivatives back to physical element
2754+ for (unsigned int r = 0; r < num_derivatives; r++)
2755+ {
2756+ for (unsigned int s = 0; s < num_derivatives; s++)
2757+ {
2758+ values[num_derivatives + r] += transform[r][s]*derivatives[s];
2759+ }// end loop over 's'
2760+ }// end loop over 'r'
2761+
2762+ // Delete pointer to array of derivatives on FIAT element
2763+ delete [] derivatives;
2764+
2765+ // Delete pointer to array of combinations of derivatives and transform
2766+ for (unsigned int r = 0; r < num_derivatives; r++)
2767+ {
2768+ delete [] combinations[r];
2769+ }// end loop over 'r'
2770+ delete [] combinations;
2771+ for (unsigned int r = 0; r < num_derivatives; r++)
2772+ {
2773+ delete [] transform[r];
2774+ }// end loop over 'r'
2775+ delete [] transform;
2776+ break;
2777+ }
2778+ case 5:
2779+ {
2780+
2781+ // Array of basisvalues
2782+ double basisvalues[3] = {0.0, 0.0, 0.0};
2783+
2784+ // Declare helper variables
2785+ double tmp0 = (1.0 + Y + 2.0*X)/2.0;
2786+
2787+ // Compute basisvalues
2788+ basisvalues[0] = 1.0;
2789+ basisvalues[1] = tmp0;
2790+ basisvalues[2] = basisvalues[0]*(0.5 + 1.5*Y);
2791+ basisvalues[0] *= std::sqrt(0.5);
2792+ basisvalues[2] *= std::sqrt(1.0);
2793+ basisvalues[1] *= std::sqrt(3.0);
2794+
2795+ // Table(s) of coefficients
2796+ static const double coefficients0[3] = \
2797+ {0.471404520791032, 0.0, 0.333333333333333};
2798+
2799+ // Tables of derivatives of the polynomial base (transpose).
2800+ static const double dmats0[3][3] = \
2801+ {{0.0, 0.0, 0.0},
2802+ {4.89897948556636, 0.0, 0.0},
2803+ {0.0, 0.0, 0.0}};
2804+
2805+ static const double dmats1[3][3] = \
2806+ {{0.0, 0.0, 0.0},
2807+ {2.44948974278318, 0.0, 0.0},
2808+ {4.24264068711928, 0.0, 0.0}};
2809+
2810+ // Compute reference derivatives.
2811+ // Declare pointer to array of derivatives on FIAT element.
2812+ double *derivatives = new double[num_derivatives];
2813+ for (unsigned int r = 0; r < num_derivatives; r++)
2814+ {
2815+ derivatives[r] = 0.0;
2816+ }// end loop over 'r'
2817+
2818+ // Declare derivative matrix (of polynomial basis).
2819+ double dmats[3][3] = \
2820+ {{1.0, 0.0, 0.0},
2821+ {0.0, 1.0, 0.0},
2822+ {0.0, 0.0, 1.0}};
2823+
2824+ // Declare (auxiliary) derivative matrix (of polynomial basis).
2825+ double dmats_old[3][3] = \
2826+ {{1.0, 0.0, 0.0},
2827+ {0.0, 1.0, 0.0},
2828+ {0.0, 0.0, 1.0}};
2829+
2830+ // Loop possible derivatives.
2831+ for (unsigned int r = 0; r < num_derivatives; r++)
2832+ {
2833+ // Resetting dmats values to compute next derivative.
2834+ for (unsigned int t = 0; t < 3; t++)
2835+ {
2836+ for (unsigned int u = 0; u < 3; u++)
2837+ {
2838+ dmats[t][u] = 0.0;
2839+ if (t == u)
2840+ {
2841+ dmats[t][u] = 1.0;
2842+ }
2843+
2844+ }// end loop over 'u'
2845+ }// end loop over 't'
2846+
2847+ // Looping derivative order to generate dmats.
2848+ for (unsigned int s = 0; s < n; s++)
2849+ {
2850+ // Updating dmats_old with new values and resetting dmats.
2851+ for (unsigned int t = 0; t < 3; t++)
2852+ {
2853+ for (unsigned int u = 0; u < 3; u++)
2854+ {
2855+ dmats_old[t][u] = dmats[t][u];
2856+ dmats[t][u] = 0.0;
2857+ }// end loop over 'u'
2858+ }// end loop over 't'
2859+
2860+ // Update dmats using an inner product.
2861+ if (combinations[r][s] == 0)
2862+ {
2863+ for (unsigned int t = 0; t < 3; t++)
2864+ {
2865+ for (unsigned int u = 0; u < 3; u++)
2866+ {
2867+ for (unsigned int tu = 0; tu < 3; tu++)
2868+ {
2869+ dmats[t][u] += dmats0[t][tu]*dmats_old[tu][u];
2870+ }// end loop over 'tu'
2871+ }// end loop over 'u'
2872+ }// end loop over 't'
2873+ }
2874+
2875+ if (combinations[r][s] == 1)
2876+ {
2877+ for (unsigned int t = 0; t < 3; t++)
2878+ {
2879+ for (unsigned int u = 0; u < 3; u++)
2880+ {
2881+ for (unsigned int tu = 0; tu < 3; tu++)
2882+ {
2883+ dmats[t][u] += dmats1[t][tu]*dmats_old[tu][u];
2884+ }// end loop over 'tu'
2885+ }// end loop over 'u'
2886+ }// end loop over 't'
2887+ }
2888+
2889+ }// end loop over 's'
2890+ for (unsigned int s = 0; s < 3; s++)
2891+ {
2892+ for (unsigned int t = 0; t < 3; t++)
2893+ {
2894+ derivatives[r] += coefficients0[s]*dmats[s][t]*basisvalues[t];
2895+ }// end loop over 't'
2896+ }// end loop over 's'
2897+ }// end loop over 'r'
2898+
2899+ // Transform derivatives back to physical element
2900+ for (unsigned int r = 0; r < num_derivatives; r++)
2901+ {
2902+ for (unsigned int s = 0; s < num_derivatives; s++)
2903+ {
2904+ values[num_derivatives + r] += transform[r][s]*derivatives[s];
2905+ }// end loop over 's'
2906+ }// end loop over 'r'
2907+
2908+ // Delete pointer to array of derivatives on FIAT element
2909+ delete [] derivatives;
2910+
2911+ // Delete pointer to array of combinations of derivatives and transform
2912+ for (unsigned int r = 0; r < num_derivatives; r++)
2913+ {
2914+ delete [] combinations[r];
2915+ }// end loop over 'r'
2916+ delete [] combinations;
2917+ for (unsigned int r = 0; r < num_derivatives; r++)
2918+ {
2919+ delete [] transform[r];
2920+ }// end loop over 'r'
2921+ delete [] transform;
2922+ break;
2923+ }
2924+ }
2925+
2926+ }
2927+
2928+ /// Evaluate order n derivatives of all basis functions at given point x in cell
2929+ virtual void evaluate_basis_derivatives_all(std::size_t n,
2930+ double* values,
2931+ const double* x,
2932+ const double* vertex_coordinates,
2933+ int cell_orientation) const
2934+ {
2935+ // Compute number of derivatives.
2936+ unsigned int num_derivatives = 1;
2937+ for (unsigned int r = 0; r < n; r++)
2938+ {
2939+ num_derivatives *= 2;
2940+ }// end loop over 'r'
2941+
2942+ // Helper variable to hold values of a single dof.
2943+ double *dof_values = new double[2*num_derivatives];
2944+ for (unsigned int r = 0; r < 2*num_derivatives; r++)
2945+ {
2946+ dof_values[r] = 0.0;
2947+ }// end loop over 'r'
2948+
2949+ // Loop dofs and call evaluate_basis_derivatives.
2950+ for (unsigned int r = 0; r < 6; r++)
2951+ {
2952+ evaluate_basis_derivatives(r, n, dof_values, x, vertex_coordinates, cell_orientation);
2953+ for (unsigned int s = 0; s < 2*num_derivatives; s++)
2954+ {
2955+ values[r*2*num_derivatives + s] = dof_values[s];
2956+ }// end loop over 's'
2957+ }// end loop over 'r'
2958+
2959+ // Delete pointer.
2960+ delete [] dof_values;
2961+ }
2962+
2963+ /// Evaluate linear functional for dof i on the function f
2964+ virtual double evaluate_dof(std::size_t i,
2965+ const ufc::function& f,
2966+ const double* vertex_coordinates,
2967+ int cell_orientation,
2968+ const ufc::cell& c) const
2969+ {
2970+ // Declare variables for result of evaluation
2971+ double vals[2];
2972+
2973+ // Declare variable for physical coordinates
2974+ double y[2];
2975+ switch (i)
2976+ {
2977+ case 0:
2978+ {
2979+ y[0] = vertex_coordinates[0];
2980+ y[1] = vertex_coordinates[1];
2981+ f.evaluate(vals, y, c);
2982+ return vals[0];
2983+ break;
2984+ }
2985+ case 1:
2986+ {
2987+ y[0] = vertex_coordinates[2];
2988+ y[1] = vertex_coordinates[3];
2989+ f.evaluate(vals, y, c);
2990+ return vals[0];
2991+ break;
2992+ }
2993+ case 2:
2994+ {
2995+ y[0] = vertex_coordinates[4];
2996+ y[1] = vertex_coordinates[5];
2997+ f.evaluate(vals, y, c);
2998+ return vals[0];
2999+ break;
3000+ }
3001+ case 3:
3002+ {
3003+ y[0] = vertex_coordinates[0];
3004+ y[1] = vertex_coordinates[1];
3005+ f.evaluate(vals, y, c);
3006+ return vals[1];
3007+ break;
3008+ }
3009+ case 4:
3010+ {
3011+ y[0] = vertex_coordinates[2];
3012+ y[1] = vertex_coordinates[3];
3013+ f.evaluate(vals, y, c);
3014+ return vals[1];
3015+ break;
3016+ }
3017+ case 5:
3018+ {
3019+ y[0] = vertex_coordinates[4];
3020+ y[1] = vertex_coordinates[5];
3021+ f.evaluate(vals, y, c);
3022+ return vals[1];
3023+ break;
3024+ }
3025+ }
3026+
3027+ return 0.0;
3028+ }
3029+
3030+ /// Evaluate linear functionals for all dofs on the function f
3031+ virtual void evaluate_dofs(double* values,
3032+ const ufc::function& f,
3033+ const double* vertex_coordinates,
3034+ int cell_orientation,
3035+ const ufc::cell& c) const
3036+ {
3037+ // Declare variables for result of evaluation
3038+ double vals[2];
3039+
3040+ // Declare variable for physical coordinates
3041+ double y[2];
3042+ y[0] = vertex_coordinates[0];
3043+ y[1] = vertex_coordinates[1];
3044+ f.evaluate(vals, y, c);
3045+ values[0] = vals[0];
3046+ y[0] = vertex_coordinates[2];
3047+ y[1] = vertex_coordinates[3];
3048+ f.evaluate(vals, y, c);
3049+ values[1] = vals[0];
3050+ y[0] = vertex_coordinates[4];
3051+ y[1] = vertex_coordinates[5];
3052+ f.evaluate(vals, y, c);
3053+ values[2] = vals[0];
3054+ y[0] = vertex_coordinates[0];
3055+ y[1] = vertex_coordinates[1];
3056+ f.evaluate(vals, y, c);
3057+ values[3] = vals[1];
3058+ y[0] = vertex_coordinates[2];
3059+ y[1] = vertex_coordinates[3];
3060+ f.evaluate(vals, y, c);
3061+ values[4] = vals[1];
3062+ y[0] = vertex_coordinates[4];
3063+ y[1] = vertex_coordinates[5];
3064+ f.evaluate(vals, y, c);
3065+ values[5] = vals[1];
3066+ }
3067+
3068+ /// Interpolate vertex values from dof values
3069+ virtual void interpolate_vertex_values(double* vertex_values,
3070+ const double* dof_values,
3071+ const double* vertex_coordinates,
3072+ int cell_orientation,
3073+ const ufc::cell& c) const
3074+ {
3075+ // Evaluate function and change variables
3076+ vertex_values[0] = dof_values[0];
3077+ vertex_values[2] = dof_values[1];
3078+ vertex_values[4] = dof_values[2];
3079+ // Evaluate function and change variables
3080+ vertex_values[1] = dof_values[3];
3081+ vertex_values[3] = dof_values[4];
3082+ vertex_values[5] = dof_values[5];
3083+ }
3084+
3085+ /// Map coordinate xhat from reference cell to coordinate x in cell
3086+ virtual void map_from_reference_cell(double* x,
3087+ const double* xhat,
3088+ const ufc::cell& c) const
3089+ {
3090+ throw std::runtime_error("map_from_reference_cell not yet implemented.");
3091+ }
3092+
3093+ /// Map from coordinate x in cell to coordinate xhat in reference cell
3094+ virtual void map_to_reference_cell(double* xhat,
3095+ const double* x,
3096+ const ufc::cell& c) const
3097+ {
3098+ throw std::runtime_error("map_to_reference_cell not yet implemented.");
3099+ }
3100+
3101+ /// Return the number of sub elements (for a mixed element)
3102+ virtual std::size_t num_sub_elements() const
3103+ {
3104+ return 2;
3105+ }
3106+
3107+ /// Create a new finite element for sub element i (for a mixed element)
3108+ virtual ufc::finite_element* create_sub_element(std::size_t i) const
3109+ {
3110+ switch (i)
3111+ {
3112+ case 0:
3113+ {
3114+ return new elasticity_finite_element_1();
3115+ break;
3116+ }
3117+ case 1:
3118+ {
3119+ return new elasticity_finite_element_1();
3120+ break;
3121+ }
3122+ }
3123+
3124+ return 0;
3125+ }
3126+
3127+ /// Create a new class instance
3128+ virtual ufc::finite_element* create() const
3129+ {
3130+ return new elasticity_finite_element_2();
3131+ }
3132+
3133+};
3134+
3135+/// This class defines the interface for a local-to-global mapping of
3136+/// degrees of freedom (dofs).
3137+
3138+class elasticity_dofmap_0: public ufc::dofmap
3139+{
3140+public:
3141+
3142+ /// Constructor
3143+ elasticity_dofmap_0() : ufc::dofmap()
3144+ {
3145+ // Do nothing
3146+ }
3147+
3148+ /// Destructor
3149+ virtual ~elasticity_dofmap_0()
3150+ {
3151+ // Do nothing
3152+ }
3153+
3154+ /// Return a string identifying the dofmap
3155+ virtual const char* signature() const
3156+ {
3157+ return "FFC dofmap for FiniteElement('Real', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 0, None)";
3158+ }
3159+
3160+ /// Return true iff mesh entities of topological dimension d are needed
3161+ virtual bool needs_mesh_entities(std::size_t d) const
3162+ {
3163+ switch (d)
3164+ {
3165+ case 0:
3166+ {
3167+ return false;
3168+ break;
3169+ }
3170+ case 1:
3171+ {
3172+ return false;
3173+ break;
3174+ }
3175+ case 2:
3176+ {
3177+ return false;
3178+ break;
3179+ }
3180+ }
3181+
3182+ return false;
3183+ }
3184+
3185+ /// Return the topological dimension of the associated cell shape
3186+ virtual std::size_t topological_dimension() const
3187+ {
3188+ return 2;
3189+ }
3190+
3191+ /// Return the geometric dimension of the associated cell shape
3192+ virtual std::size_t geometric_dimension() const
3193+ {
3194+ return 2;
3195+ }
3196+
3197+ /// Return the dimension of the global finite element function space
3198+ virtual std::size_t global_dimension(const std::vector<std::size_t>&
3199+ num_global_entities) const
3200+ {
3201+ return 1;
3202+ }
3203+
3204+ /// Return the dimension of the local finite element function space for a cell
3205+ virtual std::size_t local_dimension(const ufc::cell& c) const
3206+ {
3207+ return 1;
3208+ }
3209+
3210+ /// Return the maximum dimension of the local finite element function space
3211+ virtual std::size_t max_local_dimension() const
3212+ {
3213+ return 1;
3214+ }
3215+
3216+ /// Return the number of dofs on each cell facet
3217+ virtual std::size_t num_facet_dofs() const
3218+ {
3219+ return 0;
3220+ }
3221+
3222+ /// Return the number of dofs associated with each cell entity of dimension d
3223+ virtual std::size_t num_entity_dofs(std::size_t d) const
3224+ {
3225+ switch (d)
3226+ {
3227+ case 0:
3228+ {
3229+ return 0;
3230+ break;
3231+ }
3232+ case 1:
3233+ {
3234+ return 0;
3235+ break;
3236+ }
3237+ case 2:
3238+ {
3239+ return 1;
3240+ break;
3241+ }
3242+ }
3243+
3244+ return 0;
3245+ }
3246+
3247+ /// Tabulate the local-to-global mapping of dofs on a cell
3248+ virtual void tabulate_dofs(std::size_t* dofs,
3249+ const std::vector<std::size_t>& num_global_entities,
3250+ const ufc::cell& c) const
3251+ {
3252+ dofs[0] = 0;
3253+ }
3254+
3255+ /// Tabulate the local-to-local mapping from facet dofs to cell dofs
3256+ virtual void tabulate_facet_dofs(std::size_t* dofs,
3257+ std::size_t facet) const
3258+ {
3259+ switch (facet)
3260+ {
3261+ case 0:
3262+ {
3263+
3264+ break;
3265+ }
3266+ case 1:
3267+ {
3268+
3269+ break;
3270+ }
3271+ case 2:
3272+ {
3273+
3274+ break;
3275+ }
3276+ }
3277+
3278+ }
3279+
3280+ /// Tabulate the local-to-local mapping of dofs on entity (d, i)
3281+ virtual void tabulate_entity_dofs(std::size_t* dofs,
3282+ std::size_t d, std::size_t i) const
3283+ {
3284+ if (d > 2)
3285+ {
3286+ throw std::runtime_error("d is larger than dimension (2)");
3287+ }
3288+
3289+ switch (d)
3290+ {
3291+ case 0:
3292+ {
3293+
3294+ break;
3295+ }
3296+ case 1:
3297+ {
3298+
3299+ break;
3300+ }
3301+ case 2:
3302+ {
3303+ if (i > 0)
3304+ {
3305+ throw std::runtime_error("i is larger than number of entities (0)");
3306+ }
3307+
3308+ dofs[0] = 0;
3309+ break;
3310+ }
3311+ }
3312+
3313+ }
3314+
3315+ /// Tabulate the coordinates of all dofs on a cell
3316+ virtual void tabulate_coordinates(double** dof_coordinates,
3317+ const double* vertex_coordinates) const
3318+ {
3319+ dof_coordinates[0][0] = 0.333333333333333*vertex_coordinates[0] + 0.333333333333333*vertex_coordinates[2] + 0.333333333333333*vertex_coordinates[4];
3320+ dof_coordinates[0][1] = 0.333333333333333*vertex_coordinates[1] + 0.333333333333333*vertex_coordinates[3] + 0.333333333333333*vertex_coordinates[5];
3321+ }
3322+
3323+ /// Return the number of sub dofmaps (for a mixed element)
3324+ virtual std::size_t num_sub_dofmaps() const
3325+ {
3326+ return 0;
3327+ }
3328+
3329+ /// Create a new dofmap for sub dofmap i (for a mixed element)
3330+ virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3331+ {
3332+ return 0;
3333+ }
3334+
3335+ /// Create a new class instance
3336+ virtual ufc::dofmap* create() const
3337+ {
3338+ return new elasticity_dofmap_0();
3339+ }
3340+
3341+};
3342+
3343+/// This class defines the interface for a local-to-global mapping of
3344+/// degrees of freedom (dofs).
3345+
3346+class elasticity_dofmap_1: public ufc::dofmap
3347+{
3348+public:
3349+
3350+ /// Constructor
3351+ elasticity_dofmap_1() : ufc::dofmap()
3352+ {
3353+ // Do nothing
3354+ }
3355+
3356+ /// Destructor
3357+ virtual ~elasticity_dofmap_1()
3358+ {
3359+ // Do nothing
3360+ }
3361+
3362+ /// Return a string identifying the dofmap
3363+ virtual const char* signature() const
3364+ {
3365+ return "FFC dofmap for FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)";
3366+ }
3367+
3368+ /// Return true iff mesh entities of topological dimension d are needed
3369+ virtual bool needs_mesh_entities(std::size_t d) const
3370+ {
3371+ switch (d)
3372+ {
3373+ case 0:
3374+ {
3375+ return true;
3376+ break;
3377+ }
3378+ case 1:
3379+ {
3380+ return false;
3381+ break;
3382+ }
3383+ case 2:
3384+ {
3385+ return false;
3386+ break;
3387+ }
3388+ }
3389+
3390+ return false;
3391+ }
3392+
3393+ /// Return the topological dimension of the associated cell shape
3394+ virtual std::size_t topological_dimension() const
3395+ {
3396+ return 2;
3397+ }
3398+
3399+ /// Return the geometric dimension of the associated cell shape
3400+ virtual std::size_t geometric_dimension() const
3401+ {
3402+ return 2;
3403+ }
3404+
3405+ /// Return the dimension of the global finite element function space
3406+ virtual std::size_t global_dimension(const std::vector<std::size_t>&
3407+ num_global_entities) const
3408+ {
3409+ return num_global_entities[0];
3410+ }
3411+
3412+ /// Return the dimension of the local finite element function space for a cell
3413+ virtual std::size_t local_dimension(const ufc::cell& c) const
3414+ {
3415+ return 3;
3416+ }
3417+
3418+ /// Return the maximum dimension of the local finite element function space
3419+ virtual std::size_t max_local_dimension() const
3420+ {
3421+ return 3;
3422+ }
3423+
3424+ /// Return the number of dofs on each cell facet
3425+ virtual std::size_t num_facet_dofs() const
3426+ {
3427+ return 2;
3428+ }
3429+
3430+ /// Return the number of dofs associated with each cell entity of dimension d
3431+ virtual std::size_t num_entity_dofs(std::size_t d) const
3432+ {
3433+ switch (d)
3434+ {
3435+ case 0:
3436+ {
3437+ return 1;
3438+ break;
3439+ }
3440+ case 1:
3441+ {
3442+ return 0;
3443+ break;
3444+ }
3445+ case 2:
3446+ {
3447+ return 0;
3448+ break;
3449+ }
3450+ }
3451+
3452+ return 0;
3453+ }
3454+
3455+ /// Tabulate the local-to-global mapping of dofs on a cell
3456+ virtual void tabulate_dofs(std::size_t* dofs,
3457+ const std::vector<std::size_t>& num_global_entities,
3458+ const ufc::cell& c) const
3459+ {
3460+ dofs[0] = c.entity_indices[0][0];
3461+ dofs[1] = c.entity_indices[0][1];
3462+ dofs[2] = c.entity_indices[0][2];
3463+ }
3464+
3465+ /// Tabulate the local-to-local mapping from facet dofs to cell dofs
3466+ virtual void tabulate_facet_dofs(std::size_t* dofs,
3467+ std::size_t facet) const
3468+ {
3469+ switch (facet)
3470+ {
3471+ case 0:
3472+ {
3473+ dofs[0] = 1;
3474+ dofs[1] = 2;
3475+ break;
3476+ }
3477+ case 1:
3478+ {
3479+ dofs[0] = 0;
3480+ dofs[1] = 2;
3481+ break;
3482+ }
3483+ case 2:
3484+ {
3485+ dofs[0] = 0;
3486+ dofs[1] = 1;
3487+ break;
3488+ }
3489+ }
3490+
3491+ }
3492+
3493+ /// Tabulate the local-to-local mapping of dofs on entity (d, i)
3494+ virtual void tabulate_entity_dofs(std::size_t* dofs,
3495+ std::size_t d, std::size_t i) const
3496+ {
3497+ if (d > 2)
3498+ {
3499+ throw std::runtime_error("d is larger than dimension (2)");
3500+ }
3501+
3502+ switch (d)
3503+ {
3504+ case 0:
3505+ {
3506+ if (i > 2)
3507+ {
3508+ throw std::runtime_error("i is larger than number of entities (2)");
3509+ }
3510+
3511+ switch (i)
3512+ {
3513+ case 0:
3514+ {
3515+ dofs[0] = 0;
3516+ break;
3517+ }
3518+ case 1:
3519+ {
3520+ dofs[0] = 1;
3521+ break;
3522+ }
3523+ case 2:
3524+ {
3525+ dofs[0] = 2;
3526+ break;
3527+ }
3528+ }
3529+
3530+ break;
3531+ }
3532+ case 1:
3533+ {
3534+
3535+ break;
3536+ }
3537+ case 2:
3538+ {
3539+
3540+ break;
3541+ }
3542+ }
3543+
3544+ }
3545+
3546+ /// Tabulate the coordinates of all dofs on a cell
3547+ virtual void tabulate_coordinates(double** dof_coordinates,
3548+ const double* vertex_coordinates) const
3549+ {
3550+ dof_coordinates[0][0] = vertex_coordinates[0];
3551+ dof_coordinates[0][1] = vertex_coordinates[1];
3552+ dof_coordinates[1][0] = vertex_coordinates[2];
3553+ dof_coordinates[1][1] = vertex_coordinates[3];
3554+ dof_coordinates[2][0] = vertex_coordinates[4];
3555+ dof_coordinates[2][1] = vertex_coordinates[5];
3556+ }
3557+
3558+ /// Return the number of sub dofmaps (for a mixed element)
3559+ virtual std::size_t num_sub_dofmaps() const
3560+ {
3561+ return 0;
3562+ }
3563+
3564+ /// Create a new dofmap for sub dofmap i (for a mixed element)
3565+ virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3566+ {
3567+ return 0;
3568+ }
3569+
3570+ /// Create a new class instance
3571+ virtual ufc::dofmap* create() const
3572+ {
3573+ return new elasticity_dofmap_1();
3574+ }
3575+
3576+};
3577+
3578+/// This class defines the interface for a local-to-global mapping of
3579+/// degrees of freedom (dofs).
3580+
3581+class elasticity_dofmap_2: public ufc::dofmap
3582+{
3583+public:
3584+
3585+ /// Constructor
3586+ elasticity_dofmap_2() : ufc::dofmap()
3587+ {
3588+ // Do nothing
3589+ }
3590+
3591+ /// Destructor
3592+ virtual ~elasticity_dofmap_2()
3593+ {
3594+ // Do nothing
3595+ }
3596+
3597+ /// Return a string identifying the dofmap
3598+ virtual const char* signature() const
3599+ {
3600+ return "FFC dofmap for VectorElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, 2, None)";
3601+ }
3602+
3603+ /// Return true iff mesh entities of topological dimension d are needed
3604+ virtual bool needs_mesh_entities(std::size_t d) const
3605+ {
3606+ switch (d)
3607+ {
3608+ case 0:
3609+ {
3610+ return true;
3611+ break;
3612+ }
3613+ case 1:
3614+ {
3615+ return false;
3616+ break;
3617+ }
3618+ case 2:
3619+ {
3620+ return false;
3621+ break;
3622+ }
3623+ }
3624+
3625+ return false;
3626+ }
3627+
3628+ /// Return the topological dimension of the associated cell shape
3629+ virtual std::size_t topological_dimension() const
3630+ {
3631+ return 2;
3632+ }
3633+
3634+ /// Return the geometric dimension of the associated cell shape
3635+ virtual std::size_t geometric_dimension() const
3636+ {
3637+ return 2;
3638+ }
3639+
3640+ /// Return the dimension of the global finite element function space
3641+ virtual std::size_t global_dimension(const std::vector<std::size_t>&
3642+ num_global_entities) const
3643+ {
3644+ return 2*num_global_entities[0];
3645+ }
3646+
3647+ /// Return the dimension of the local finite element function space for a cell
3648+ virtual std::size_t local_dimension(const ufc::cell& c) const
3649+ {
3650+ return 6;
3651+ }
3652+
3653+ /// Return the maximum dimension of the local finite element function space
3654+ virtual std::size_t max_local_dimension() const
3655+ {
3656+ return 6;
3657+ }
3658+
3659+ /// Return the number of dofs on each cell facet
3660+ virtual std::size_t num_facet_dofs() const
3661+ {
3662+ return 4;
3663+ }
3664+
3665+ /// Return the number of dofs associated with each cell entity of dimension d
3666+ virtual std::size_t num_entity_dofs(std::size_t d) const
3667+ {
3668+ switch (d)
3669+ {
3670+ case 0:
3671+ {
3672+ return 2;
3673+ break;
3674+ }
3675+ case 1:
3676+ {
3677+ return 0;
3678+ break;
3679+ }
3680+ case 2:
3681+ {
3682+ return 0;
3683+ break;
3684+ }
3685+ }
3686+
3687+ return 0;
3688+ }
3689+
3690+ /// Tabulate the local-to-global mapping of dofs on a cell
3691+ virtual void tabulate_dofs(std::size_t* dofs,
3692+ const std::vector<std::size_t>& num_global_entities,
3693+ const ufc::cell& c) const
3694+ {
3695+ unsigned int offset = 0;
3696+ dofs[0] = offset + c.entity_indices[0][0];
3697+ dofs[1] = offset + c.entity_indices[0][1];
3698+ dofs[2] = offset + c.entity_indices[0][2];
3699+ offset += num_global_entities[0];
3700+ dofs[3] = offset + c.entity_indices[0][0];
3701+ dofs[4] = offset + c.entity_indices[0][1];
3702+ dofs[5] = offset + c.entity_indices[0][2];
3703+ offset += num_global_entities[0];
3704+ }
3705+
3706+ /// Tabulate the local-to-local mapping from facet dofs to cell dofs
3707+ virtual void tabulate_facet_dofs(std::size_t* dofs,
3708+ std::size_t facet) const
3709+ {
3710+ switch (facet)
3711+ {
3712+ case 0:
3713+ {
3714+ dofs[0] = 1;
3715+ dofs[1] = 2;
3716+ dofs[2] = 4;
3717+ dofs[3] = 5;
3718+ break;
3719+ }
3720+ case 1:
3721+ {
3722+ dofs[0] = 0;
3723+ dofs[1] = 2;
3724+ dofs[2] = 3;
3725+ dofs[3] = 5;
3726+ break;
3727+ }
3728+ case 2:
3729+ {
3730+ dofs[0] = 0;
3731+ dofs[1] = 1;
3732+ dofs[2] = 3;
3733+ dofs[3] = 4;
3734+ break;
3735+ }
3736+ }
3737+
3738+ }
3739+
3740+ /// Tabulate the local-to-local mapping of dofs on entity (d, i)
3741+ virtual void tabulate_entity_dofs(std::size_t* dofs,
3742+ std::size_t d, std::size_t i) const
3743+ {
3744+ if (d > 2)
3745+ {
3746+ throw std::runtime_error("d is larger than dimension (2)");
3747+ }
3748+
3749+ switch (d)
3750+ {
3751+ case 0:
3752+ {
3753+ if (i > 2)
3754+ {
3755+ throw std::runtime_error("i is larger than number of entities (2)");
3756+ }
3757+
3758+ switch (i)
3759+ {
3760+ case 0:
3761+ {
3762+ dofs[0] = 0;
3763+ dofs[1] = 3;
3764+ break;
3765+ }
3766+ case 1:
3767+ {
3768+ dofs[0] = 1;
3769+ dofs[1] = 4;
3770+ break;
3771+ }
3772+ case 2:
3773+ {
3774+ dofs[0] = 2;
3775+ dofs[1] = 5;
3776+ break;
3777+ }
3778+ }
3779+
3780+ break;
3781+ }
3782+ case 1:
3783+ {
3784+
3785+ break;
3786+ }
3787+ case 2:
3788+ {
3789+
3790+ break;
3791+ }
3792+ }
3793+
3794+ }
3795+
3796+ /// Tabulate the coordinates of all dofs on a cell
3797+ virtual void tabulate_coordinates(double** dof_coordinates,
3798+ const double* vertex_coordinates) const
3799+ {
3800+ dof_coordinates[0][0] = vertex_coordinates[0];
3801+ dof_coordinates[0][1] = vertex_coordinates[1];
3802+ dof_coordinates[1][0] = vertex_coordinates[2];
3803+ dof_coordinates[1][1] = vertex_coordinates[3];
3804+ dof_coordinates[2][0] = vertex_coordinates[4];
3805+ dof_coordinates[2][1] = vertex_coordinates[5];
3806+ dof_coordinates[3][0] = vertex_coordinates[0];
3807+ dof_coordinates[3][1] = vertex_coordinates[1];
3808+ dof_coordinates[4][0] = vertex_coordinates[2];
3809+ dof_coordinates[4][1] = vertex_coordinates[3];
3810+ dof_coordinates[5][0] = vertex_coordinates[4];
3811+ dof_coordinates[5][1] = vertex_coordinates[5];
3812+ }
3813+
3814+ /// Return the number of sub dofmaps (for a mixed element)
3815+ virtual std::size_t num_sub_dofmaps() const
3816+ {
3817+ return 2;
3818+ }
3819+
3820+ /// Create a new dofmap for sub dofmap i (for a mixed element)
3821+ virtual ufc::dofmap* create_sub_dofmap(std::size_t i) const
3822+ {
3823+ switch (i)
3824+ {
3825+ case 0:
3826+ {
3827+ return new elasticity_dofmap_1();
3828+ break;
3829+ }
3830+ case 1:
3831+ {
3832+ return new elasticity_dofmap_1();
3833+ break;
3834+ }
3835+ }
3836+
3837+ return 0;
3838+ }
3839+
3840+ /// Create a new class instance
3841+ virtual ufc::dofmap* create() const
3842+ {
3843+ return new elasticity_dofmap_2();
3844+ }
3845+
3846+};
3847+
3848+/// This class defines the interface for the tabulation of the cell
3849+/// tensor corresponding to the local contribution to a form from
3850+/// the integral over a cell.
3851+
3852+class elasticity_cell_integral_0_otherwise: public ufc::cell_integral
3853+{
3854+public:
3855+
3856+ /// Constructor
3857+ elasticity_cell_integral_0_otherwise() : ufc::cell_integral()
3858+ {
3859+ // Do nothing
3860+ }
3861+
3862+ /// Destructor
3863+ virtual ~elasticity_cell_integral_0_otherwise()
3864+ {
3865+ // Do nothing
3866+ }
3867+
3868+ /// Tabulate the tensor for the contribution from a local cell
3869+ virtual void tabulate_tensor(double* A,
3870+ const double * const * w,
3871+ const double* vertex_coordinates,
3872+ int cell_orientation) const
3873+ {
3874+ // Compute Jacobian
3875+ double J[4];
3876+ compute_jacobian_triangle_2d(J, vertex_coordinates);
3877+
3878+ // Compute Jacobian inverse and determinant
3879+ double K[4];
3880+ double detJ;
3881+ compute_jacobian_inverse_triangle_2d(K, detJ, J);
3882+
3883+ // Set scale factor
3884+ const double det = std::abs(detJ);
3885+
3886+ // Cell volume
3887+
3888+ // Compute circumradius of triangle in 2D
3889+
3890+
3891+ // Facet area
3892+
3893+ // Array of quadrature weights.
3894+ static const double W1 = 0.5;
3895+ // Quadrature points on the UFC reference element: (0.333333333333333, 0.333333333333333)
3896+
3897+ // Value of basis functions at quadrature points.
3898+ static const double FE0_C0_D01[1][6] = \
3899+ {{-1.0, 0.0, 1.0, 0.0, 0.0, 0.0}};
3900+
3901+ static const double FE0_C0_D10[1][6] = \
3902+ {{-1.0, 1.0, 0.0, 0.0, 0.0, 0.0}};
3903+
3904+ static const double FE0_C1_D01[1][6] = \
3905+ {{0.0, 0.0, 0.0, -1.0, 0.0, 1.0}};
3906+
3907+ static const double FE0_C1_D10[1][6] = \
3908+ {{0.0, 0.0, 0.0, -1.0, 1.0, 0.0}};
3909+
3910+ // Reset values in the element tensor.
3911+ for (unsigned int r = 0; r < 36; r++)
3912+ {
3913+ A[r] = 0.0;
3914+ }// end loop over 'r'
3915+
3916+ // Compute element tensor using UFL quadrature representation
3917+ // Optimisations: ('eliminate zeros', False), ('ignore ones', False), ('ignore zero tables', False), ('optimisation', False), ('remove zero terms', False)
3918+
3919+ // Loop quadrature points for integral.
3920+ // Number of operations to compute element tensor for following IP loop = 3528
3921+ // Only 1 integration point, omitting IP loop.
3922+
3923+ // Number of operations for primary indices: 3528
3924+ for (unsigned int j = 0; j < 6; j++)
3925+ {
3926+ for (unsigned int k = 0; k < 6; k++)
3927+ {
3928+ // Number of operations to compute entry: 98
3929+ A[j*6 + k] += ((((((K[1]*FE0_C0_D10[0][j] + K[3]*FE0_C0_D01[0][j]) + (K[0]*FE0_C1_D10[0][j] + K[2]*FE0_C1_D01[0][j])))*0.5)*(((((K[1]*FE0_C0_D10[0][k] + K[3]*FE0_C0_D01[0][k]) + (K[0]*FE0_C1_D10[0][k] + K[2]*FE0_C1_D01[0][k])))*0.5)*2.0*w[0][0]) + ((2.0*((K[1]*FE0_C1_D10[0][j] + K[3]*FE0_C1_D01[0][j])))*0.5)*((((2.0*((K[1]*FE0_C1_D10[0][k] + K[3]*FE0_C1_D01[0][k])))*0.5)*2.0*w[0][0] + ((((2.0*((K[1]*FE0_C1_D10[0][k] + K[3]*FE0_C1_D01[0][k])))*0.5 + (2.0*((K[0]*FE0_C0_D10[0][k] + K[2]*FE0_C0_D01[0][k])))*0.5))*w[1][0])*1.0))) + (((((K[0]*FE0_C1_D10[0][j] + K[2]*FE0_C1_D01[0][j]) + (K[1]*FE0_C0_D10[0][j] + K[3]*FE0_C0_D01[0][j])))*0.5)*(((((K[0]*FE0_C1_D10[0][k] + K[2]*FE0_C1_D01[0][k]) + (K[1]*FE0_C0_D10[0][k] + K[3]*FE0_C0_D01[0][k])))*0.5)*2.0*w[0][0]) + ((2.0*((K[0]*FE0_C0_D10[0][j] + K[2]*FE0_C0_D01[0][j])))*0.5)*((((2.0*((K[0]*FE0_C0_D10[0][k] + K[2]*FE0_C0_D01[0][k])))*0.5)*2.0*w[0][0] + ((((2.0*((K[1]*FE0_C1_D10[0][k] + K[3]*FE0_C1_D01[0][k])))*0.5 + (2.0*((K[0]*FE0_C0_D10[0][k] + K[2]*FE0_C0_D01[0][k])))*0.5))*w[1][0])*1.0))))*W1*det;
3930+ }// end loop over 'k'
3931+ }// end loop over 'j'
3932+ }
3933+
3934+};
3935+
3936+/// This class defines the interface for the tabulation of the cell
3937+/// tensor corresponding to the local contribution to a form from
3938+/// the integral over a cell.
3939+
3940+class elasticity_cell_integral_1_otherwise: public ufc::cell_integral
3941+{
3942+public:
3943+
3944+ /// Constructor
3945+ elasticity_cell_integral_1_otherwise() : ufc::cell_integral()
3946+ {
3947+ // Do nothing
3948+ }
3949+
3950+ /// Destructor
3951+ virtual ~elasticity_cell_integral_1_otherwise()
3952+ {
3953+ // Do nothing
3954+ }
3955+
3956+ /// Tabulate the tensor for the contribution from a local cell
3957+ virtual void tabulate_tensor(double* A,
3958+ const double * const * w,
3959+ const double* vertex_coordinates,
3960+ int cell_orientation) const
3961+ {
3962+ // Number of operations (multiply-add pairs) for Jacobian data: 3
3963+ // Number of operations (multiply-add pairs) for geometry tensor: 6
3964+ // Number of operations (multiply-add pairs) for tensor contraction: 15
3965+ // Total number of operations (multiply-add pairs): 24
3966+
3967+ // Compute Jacobian
3968+ double J[4];
3969+ compute_jacobian_triangle_2d(J, vertex_coordinates);
3970+
3971+ // Compute Jacobian inverse and determinant
3972+ double K[4];
3973+ double detJ;
3974+ compute_jacobian_inverse_triangle_2d(K, detJ, J);
3975+
3976+ // Set scale factor
3977+ const double det = std::abs(detJ);
3978+
3979+ // Compute geometry tensor
3980+ const double G0_0 = det*w[0][0]*(1.0);
3981+ const double G0_1 = det*w[0][1]*(1.0);
3982+ const double G0_2 = det*w[0][2]*(1.0);
3983+ const double G0_3 = det*w[0][3]*(1.0);
3984+ const double G0_4 = det*w[0][4]*(1.0);
3985+ const double G0_5 = det*w[0][5]*(1.0);
3986+
3987+ // Compute element tensor
3988+ A[0] = 0.0833333333333334*G0_0 + 0.0416666666666667*G0_1 + 0.0416666666666667*G0_2;
3989+ A[1] = 0.0416666666666667*G0_0 + 0.0833333333333333*G0_1 + 0.0416666666666666*G0_2;
3990+ A[2] = 0.0416666666666667*G0_0 + 0.0416666666666666*G0_1 + 0.0833333333333333*G0_2;
3991+ A[3] = 0.0833333333333334*G0_3 + 0.0416666666666667*G0_4 + 0.0416666666666667*G0_5;
3992+ A[4] = 0.0416666666666667*G0_3 + 0.0833333333333333*G0_4 + 0.0416666666666666*G0_5;
3993+ A[5] = 0.0416666666666667*G0_3 + 0.0416666666666666*G0_4 + 0.0833333333333333*G0_5;
3994+ }
3995+
3996+};
3997+
3998+/// This class defines the interface for the assembly of the global
3999+/// tensor corresponding to a form with r + n arguments, that is, a
4000+/// mapping
4001+///
4002+/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
4003+///
4004+/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
4005+/// global tensor A is defined by
4006+///
4007+/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
4008+///
4009+/// where each argument Vj represents the application to the
4010+/// sequence of basis functions of Vj and w1, w2, ..., wn are given
4011+/// fixed functions (coefficients).
4012+
4013+class elasticity_form_0: public ufc::form
4014+{
4015+public:
4016+
4017+ /// Constructor
4018+ elasticity_form_0() : ufc::form()
4019+ {
4020+ // Do nothing
4021+ }
4022+
4023+ /// Destructor
4024+ virtual ~elasticity_form_0()
4025+ {
4026+ // Do nothing
4027+ }
4028+
4029+ /// Return a string identifying the form
4030+ virtual const char* signature() const
4031+ {
4032+ return "ab5e8e2538548572ca24be768ed54ee22ed28f77b71b9b99587eb0d3c0bc2ad1375883fbfe073a10a12e9152fcd5a65f57529736c8217b158f2ca24377e2ff1e";
4033+ }
4034+
4035+ /// Return the rank of the global tensor (r)
4036+ virtual std::size_t rank() const
4037+ {
4038+ return 2;
4039+ }
4040+
4041+ /// Return the number of coefficients (n)
4042+ virtual std::size_t num_coefficients() const
4043+ {
4044+ return 2;
4045+ }
4046+
4047+ /// Return the number of cell domains
4048+ virtual std::size_t num_cell_domains() const
4049+ {
4050+ return 0;
4051+ }
4052+
4053+ /// Return the number of exterior facet domains
4054+ virtual std::size_t num_exterior_facet_domains() const
4055+ {
4056+ return 0;
4057+ }
4058+
4059+ /// Return the number of interior facet domains
4060+ virtual std::size_t num_interior_facet_domains() const
4061+ {
4062+ return 0;
4063+ }
4064+
4065+ /// Return the number of point domains
4066+ virtual std::size_t num_point_domains() const
4067+ {
4068+ return 0;
4069+ }
4070+
4071+ /// Return whether the form has any cell integrals
4072+ virtual bool has_cell_integrals() const
4073+ {
4074+ return true;
4075+ }
4076+
4077+ /// Return whether the form has any exterior facet integrals
4078+ virtual bool has_exterior_facet_integrals() const
4079+ {
4080+ return false;
4081+ }
4082+
4083+ /// Return whether the form has any interior facet integrals
4084+ virtual bool has_interior_facet_integrals() const
4085+ {
4086+ return false;
4087+ }
4088+
4089+ /// Return whether the form has any point integrals
4090+ virtual bool has_point_integrals() const
4091+ {
4092+ return false;
4093+ }
4094+
4095+ /// Create a new finite element for argument function i
4096+ virtual ufc::finite_element* create_finite_element(std::size_t i) const
4097+ {
4098+ switch (i)
4099+ {
4100+ case 0:
4101+ {
4102+ return new elasticity_finite_element_2();
4103+ break;
4104+ }
4105+ case 1:
4106+ {
4107+ return new elasticity_finite_element_2();
4108+ break;
4109+ }
4110+ case 2:
4111+ {
4112+ return new elasticity_finite_element_0();
4113+ break;
4114+ }
4115+ case 3:
4116+ {
4117+ return new elasticity_finite_element_0();
4118+ break;
4119+ }
4120+ }
4121+
4122+ return 0;
4123+ }
4124+
4125+ /// Create a new dofmap for argument function i
4126+ virtual ufc::dofmap* create_dofmap(std::size_t i) const
4127+ {
4128+ switch (i)
4129+ {
4130+ case 0:
4131+ {
4132+ return new elasticity_dofmap_2();
4133+ break;
4134+ }
4135+ case 1:
4136+ {
4137+ return new elasticity_dofmap_2();
4138+ break;
4139+ }
4140+ case 2:
4141+ {
4142+ return new elasticity_dofmap_0();
4143+ break;
4144+ }
4145+ case 3:
4146+ {
4147+ return new elasticity_dofmap_0();
4148+ break;
4149+ }
4150+ }
4151+
4152+ return 0;
4153+ }
4154+
4155+ /// Create a new cell integral on sub domain i
4156+ virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
4157+ {
4158+ return 0;
4159+ }
4160+
4161+ /// Create a new exterior facet integral on sub domain i
4162+ virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
4163+ {
4164+ return 0;
4165+ }
4166+
4167+ /// Create a new interior facet integral on sub domain i
4168+ virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
4169+ {
4170+ return 0;
4171+ }
4172+
4173+ /// Create a new point integral on sub domain i
4174+ virtual ufc::point_integral* create_point_integral(std::size_t i) const
4175+ {
4176+ return 0;
4177+ }
4178+
4179+ /// Create a new cell integral on everywhere else
4180+ virtual ufc::cell_integral* create_default_cell_integral() const
4181+ {
4182+ return new elasticity_cell_integral_0_otherwise();
4183+ }
4184+
4185+ /// Create a new exterior facet integral on everywhere else
4186+ virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
4187+ {
4188+ return 0;
4189+ }
4190+
4191+ /// Create a new interior facet integral on everywhere else
4192+ virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
4193+ {
4194+ return 0;
4195+ }
4196+
4197+ /// Create a new point integral on everywhere else
4198+ virtual ufc::point_integral* create_default_point_integral() const
4199+ {
4200+ return 0;
4201+ }
4202+
4203+};
4204+
4205+/// This class defines the interface for the assembly of the global
4206+/// tensor corresponding to a form with r + n arguments, that is, a
4207+/// mapping
4208+///
4209+/// a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
4210+///
4211+/// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
4212+/// global tensor A is defined by
4213+///
4214+/// A = a(V1, V2, ..., Vr, w1, w2, ..., wn),
4215+///
4216+/// where each argument Vj represents the application to the
4217+/// sequence of basis functions of Vj and w1, w2, ..., wn are given
4218+/// fixed functions (coefficients).
4219+
4220+class elasticity_form_1: public ufc::form
4221+{
4222+public:
4223+
4224+ /// Constructor
4225+ elasticity_form_1() : ufc::form()
4226+ {
4227+ // Do nothing
4228+ }
4229+
4230+ /// Destructor
4231+ virtual ~elasticity_form_1()
4232+ {
4233+ // Do nothing
4234+ }
4235+
4236+ /// Return a string identifying the form
4237+ virtual const char* signature() const
4238+ {
4239+ return "186114385784b40ca08bfe5debc0035ccff00f9604a40db76b1933e3593f3d5db3a55dfb7312d1c8663114a4d391e1676b1eba5eaac07058fa8c31a5a452cced";
4240+ }
4241+
4242+ /// Return the rank of the global tensor (r)
4243+ virtual std::size_t rank() const
4244+ {
4245+ return 1;
4246+ }
4247+
4248+ /// Return the number of coefficients (n)
4249+ virtual std::size_t num_coefficients() const
4250+ {
4251+ return 1;
4252+ }
4253+
4254+ /// Return the number of cell domains
4255+ virtual std::size_t num_cell_domains() const
4256+ {
4257+ return 0;
4258+ }
4259+
4260+ /// Return the number of exterior facet domains
4261+ virtual std::size_t num_exterior_facet_domains() const
4262+ {
4263+ return 0;
4264+ }
4265+
4266+ /// Return the number of interior facet domains
4267+ virtual std::size_t num_interior_facet_domains() const
4268+ {
4269+ return 0;
4270+ }
4271+
4272+ /// Return the number of point domains
4273+ virtual std::size_t num_point_domains() const
4274+ {
4275+ return 0;
4276+ }
4277+
4278+ /// Return whether the form has any cell integrals
4279+ virtual bool has_cell_integrals() const
4280+ {
4281+ return true;
4282+ }
4283+
4284+ /// Return whether the form has any exterior facet integrals
4285+ virtual bool has_exterior_facet_integrals() const
4286+ {
4287+ return false;
4288+ }
4289+
4290+ /// Return whether the form has any interior facet integrals
4291+ virtual bool has_interior_facet_integrals() const
4292+ {
4293+ return false;
4294+ }
4295+
4296+ /// Return whether the form has any point integrals
4297+ virtual bool has_point_integrals() const
4298+ {
4299+ return false;
4300+ }
4301+
4302+ /// Create a new finite element for argument function i
4303+ virtual ufc::finite_element* create_finite_element(std::size_t i) const
4304+ {
4305+ switch (i)
4306+ {
4307+ case 0:
4308+ {
4309+ return new elasticity_finite_element_2();
4310+ break;
4311+ }
4312+ case 1:
4313+ {
4314+ return new elasticity_finite_element_2();
4315+ break;
4316+ }
4317+ }
4318+
4319+ return 0;
4320+ }
4321+
4322+ /// Create a new dofmap for argument function i
4323+ virtual ufc::dofmap* create_dofmap(std::size_t i) const
4324+ {
4325+ switch (i)
4326+ {
4327+ case 0:
4328+ {
4329+ return new elasticity_dofmap_2();
4330+ break;
4331+ }
4332+ case 1:
4333+ {
4334+ return new elasticity_dofmap_2();
4335+ break;
4336+ }
4337+ }
4338+
4339+ return 0;
4340+ }
4341+
4342+ /// Create a new cell integral on sub domain i
4343+ virtual ufc::cell_integral* create_cell_integral(std::size_t i) const
4344+ {
4345+ return 0;
4346+ }
4347+
4348+ /// Create a new exterior facet integral on sub domain i
4349+ virtual ufc::exterior_facet_integral* create_exterior_facet_integral(std::size_t i) const
4350+ {
4351+ return 0;
4352+ }
4353+
4354+ /// Create a new interior facet integral on sub domain i
4355+ virtual ufc::interior_facet_integral* create_interior_facet_integral(std::size_t i) const
4356+ {
4357+ return 0;
4358+ }
4359+
4360+ /// Create a new point integral on sub domain i
4361+ virtual ufc::point_integral* create_point_integral(std::size_t i) const
4362+ {
4363+ return 0;
4364+ }
4365+
4366+ /// Create a new cell integral on everywhere else
4367+ virtual ufc::cell_integral* create_default_cell_integral() const
4368+ {
4369+ return new elasticity_cell_integral_1_otherwise();
4370+ }
4371+
4372+ /// Create a new exterior facet integral on everywhere else
4373+ virtual ufc::exterior_facet_integral* create_default_exterior_facet_integral() const
4374+ {
4375+ return 0;
4376+ }
4377+
4378+ /// Create a new interior facet integral on everywhere else
4379+ virtual ufc::interior_facet_integral* create_default_interior_facet_integral() const
4380+ {
4381+ return 0;
4382+ }
4383+
4384+ /// Create a new point integral on everywhere else
4385+ virtual ufc::point_integral* create_default_point_integral() const
4386+ {
4387+ return 0;
4388+ }
4389+
4390+};
4391+
4392+// DOLFIN wrappers
4393+
4394+// Standard library includes
4395+#include <string>
4396+
4397+// DOLFIN includes
4398+#include <dolfin/common/NoDeleter.h>
4399+#include <dolfin/mesh/Restriction.h>
4400+#include <dolfin/fem/FiniteElement.h>
4401+#include <dolfin/fem/DofMap.h>
4402+#include <dolfin/fem/Form.h>
4403+#include <dolfin/function/FunctionSpace.h>
4404+#include <dolfin/function/GenericFunction.h>
4405+#include <dolfin/function/CoefficientAssigner.h>
4406+#include <dolfin/adaptivity/ErrorControl.h>
4407+#include <dolfin/adaptivity/GoalFunctional.h>
4408+
4409+namespace Elasticity
4410+{
4411+
4412+class CoefficientSpace_f: public dolfin::FunctionSpace
4413+{
4414+public:
4415+
4416+ //--- Constructors for standard function space, 2 different versions ---
4417+
4418+ // Create standard function space (reference version)
4419+ CoefficientSpace_f(const dolfin::Mesh& mesh):
4420+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4421+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4422+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), mesh)))
4423+ {
4424+ // Do nothing
4425+ }
4426+
4427+ // Create standard function space (shared pointer version)
4428+ CoefficientSpace_f(boost::shared_ptr<const dolfin::Mesh> mesh):
4429+ dolfin::FunctionSpace(mesh,
4430+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4431+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), *mesh)))
4432+ {
4433+ // Do nothing
4434+ }
4435+
4436+ //--- Constructors for constrained function space, 2 different versions ---
4437+
4438+ // Create standard function space (reference version)
4439+ CoefficientSpace_f(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
4440+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4441+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4442+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), mesh,
4443+ dolfin::reference_to_no_delete_pointer(constrained_domain))))
4444+ {
4445+ // Do nothing
4446+ }
4447+
4448+ // Create standard function space (shared pointer version)
4449+ CoefficientSpace_f(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
4450+ dolfin::FunctionSpace(mesh,
4451+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4452+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), *mesh, constrained_domain)))
4453+ {
4454+ // Do nothing
4455+ }
4456+
4457+ //--- Constructors for restricted function space, 2 different versions ---
4458+
4459+ // Create restricted function space (reference version)
4460+ CoefficientSpace_f(const dolfin::Restriction& restriction):
4461+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
4462+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4463+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()),
4464+ reference_to_no_delete_pointer(restriction))))
4465+ {
4466+ // Do nothing
4467+ }
4468+
4469+ // Create restricted function space (shared pointer version)
4470+ CoefficientSpace_f(boost::shared_ptr<const dolfin::Restriction> restriction):
4471+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
4472+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4473+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()),
4474+ restriction)))
4475+ {
4476+ // Do nothing
4477+ }
4478+
4479+ // Copy constructor
4480+ ~CoefficientSpace_f()
4481+ {
4482+ }
4483+
4484+};
4485+
4486+class CoefficientSpace_lmbda: public dolfin::FunctionSpace
4487+{
4488+public:
4489+
4490+ //--- Constructors for standard function space, 2 different versions ---
4491+
4492+ // Create standard function space (reference version)
4493+ CoefficientSpace_lmbda(const dolfin::Mesh& mesh):
4494+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4495+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4496+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()), mesh)))
4497+ {
4498+ // Do nothing
4499+ }
4500+
4501+ // Create standard function space (shared pointer version)
4502+ CoefficientSpace_lmbda(boost::shared_ptr<const dolfin::Mesh> mesh):
4503+ dolfin::FunctionSpace(mesh,
4504+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4505+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()), *mesh)))
4506+ {
4507+ // Do nothing
4508+ }
4509+
4510+ //--- Constructors for constrained function space, 2 different versions ---
4511+
4512+ // Create standard function space (reference version)
4513+ CoefficientSpace_lmbda(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
4514+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4515+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4516+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()), mesh,
4517+ dolfin::reference_to_no_delete_pointer(constrained_domain))))
4518+ {
4519+ // Do nothing
4520+ }
4521+
4522+ // Create standard function space (shared pointer version)
4523+ CoefficientSpace_lmbda(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
4524+ dolfin::FunctionSpace(mesh,
4525+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4526+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()), *mesh, constrained_domain)))
4527+ {
4528+ // Do nothing
4529+ }
4530+
4531+ //--- Constructors for restricted function space, 2 different versions ---
4532+
4533+ // Create restricted function space (reference version)
4534+ CoefficientSpace_lmbda(const dolfin::Restriction& restriction):
4535+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
4536+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4537+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()),
4538+ reference_to_no_delete_pointer(restriction))))
4539+ {
4540+ // Do nothing
4541+ }
4542+
4543+ // Create restricted function space (shared pointer version)
4544+ CoefficientSpace_lmbda(boost::shared_ptr<const dolfin::Restriction> restriction):
4545+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
4546+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4547+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()),
4548+ restriction)))
4549+ {
4550+ // Do nothing
4551+ }
4552+
4553+ // Copy constructor
4554+ ~CoefficientSpace_lmbda()
4555+ {
4556+ }
4557+
4558+};
4559+
4560+class CoefficientSpace_mu: public dolfin::FunctionSpace
4561+{
4562+public:
4563+
4564+ //--- Constructors for standard function space, 2 different versions ---
4565+
4566+ // Create standard function space (reference version)
4567+ CoefficientSpace_mu(const dolfin::Mesh& mesh):
4568+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4569+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4570+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()), mesh)))
4571+ {
4572+ // Do nothing
4573+ }
4574+
4575+ // Create standard function space (shared pointer version)
4576+ CoefficientSpace_mu(boost::shared_ptr<const dolfin::Mesh> mesh):
4577+ dolfin::FunctionSpace(mesh,
4578+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4579+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()), *mesh)))
4580+ {
4581+ // Do nothing
4582+ }
4583+
4584+ //--- Constructors for constrained function space, 2 different versions ---
4585+
4586+ // Create standard function space (reference version)
4587+ CoefficientSpace_mu(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
4588+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4589+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4590+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()), mesh,
4591+ dolfin::reference_to_no_delete_pointer(constrained_domain))))
4592+ {
4593+ // Do nothing
4594+ }
4595+
4596+ // Create standard function space (shared pointer version)
4597+ CoefficientSpace_mu(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
4598+ dolfin::FunctionSpace(mesh,
4599+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4600+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()), *mesh, constrained_domain)))
4601+ {
4602+ // Do nothing
4603+ }
4604+
4605+ //--- Constructors for restricted function space, 2 different versions ---
4606+
4607+ // Create restricted function space (reference version)
4608+ CoefficientSpace_mu(const dolfin::Restriction& restriction):
4609+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
4610+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4611+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()),
4612+ reference_to_no_delete_pointer(restriction))))
4613+ {
4614+ // Do nothing
4615+ }
4616+
4617+ // Create restricted function space (shared pointer version)
4618+ CoefficientSpace_mu(boost::shared_ptr<const dolfin::Restriction> restriction):
4619+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
4620+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_0()))),
4621+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_0()),
4622+ restriction)))
4623+ {
4624+ // Do nothing
4625+ }
4626+
4627+ // Copy constructor
4628+ ~CoefficientSpace_mu()
4629+ {
4630+ }
4631+
4632+};
4633+
4634+class Form_a_FunctionSpace_0: public dolfin::FunctionSpace
4635+{
4636+public:
4637+
4638+ //--- Constructors for standard function space, 2 different versions ---
4639+
4640+ // Create standard function space (reference version)
4641+ Form_a_FunctionSpace_0(const dolfin::Mesh& mesh):
4642+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4643+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4644+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), mesh)))
4645+ {
4646+ // Do nothing
4647+ }
4648+
4649+ // Create standard function space (shared pointer version)
4650+ Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh):
4651+ dolfin::FunctionSpace(mesh,
4652+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4653+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), *mesh)))
4654+ {
4655+ // Do nothing
4656+ }
4657+
4658+ //--- Constructors for constrained function space, 2 different versions ---
4659+
4660+ // Create standard function space (reference version)
4661+ Form_a_FunctionSpace_0(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
4662+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4663+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4664+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), mesh,
4665+ dolfin::reference_to_no_delete_pointer(constrained_domain))))
4666+ {
4667+ // Do nothing
4668+ }
4669+
4670+ // Create standard function space (shared pointer version)
4671+ Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
4672+ dolfin::FunctionSpace(mesh,
4673+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4674+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), *mesh, constrained_domain)))
4675+ {
4676+ // Do nothing
4677+ }
4678+
4679+ //--- Constructors for restricted function space, 2 different versions ---
4680+
4681+ // Create restricted function space (reference version)
4682+ Form_a_FunctionSpace_0(const dolfin::Restriction& restriction):
4683+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
4684+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4685+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()),
4686+ reference_to_no_delete_pointer(restriction))))
4687+ {
4688+ // Do nothing
4689+ }
4690+
4691+ // Create restricted function space (shared pointer version)
4692+ Form_a_FunctionSpace_0(boost::shared_ptr<const dolfin::Restriction> restriction):
4693+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
4694+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4695+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()),
4696+ restriction)))
4697+ {
4698+ // Do nothing
4699+ }
4700+
4701+ // Copy constructor
4702+ ~Form_a_FunctionSpace_0()
4703+ {
4704+ }
4705+
4706+};
4707+
4708+class Form_a_FunctionSpace_1: public dolfin::FunctionSpace
4709+{
4710+public:
4711+
4712+ //--- Constructors for standard function space, 2 different versions ---
4713+
4714+ // Create standard function space (reference version)
4715+ Form_a_FunctionSpace_1(const dolfin::Mesh& mesh):
4716+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4717+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4718+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), mesh)))
4719+ {
4720+ // Do nothing
4721+ }
4722+
4723+ // Create standard function space (shared pointer version)
4724+ Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh):
4725+ dolfin::FunctionSpace(mesh,
4726+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4727+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), *mesh)))
4728+ {
4729+ // Do nothing
4730+ }
4731+
4732+ //--- Constructors for constrained function space, 2 different versions ---
4733+
4734+ // Create standard function space (reference version)
4735+ Form_a_FunctionSpace_1(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
4736+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4737+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4738+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), mesh,
4739+ dolfin::reference_to_no_delete_pointer(constrained_domain))))
4740+ {
4741+ // Do nothing
4742+ }
4743+
4744+ // Create standard function space (shared pointer version)
4745+ Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
4746+ dolfin::FunctionSpace(mesh,
4747+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4748+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), *mesh, constrained_domain)))
4749+ {
4750+ // Do nothing
4751+ }
4752+
4753+ //--- Constructors for restricted function space, 2 different versions ---
4754+
4755+ // Create restricted function space (reference version)
4756+ Form_a_FunctionSpace_1(const dolfin::Restriction& restriction):
4757+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
4758+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4759+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()),
4760+ reference_to_no_delete_pointer(restriction))))
4761+ {
4762+ // Do nothing
4763+ }
4764+
4765+ // Create restricted function space (shared pointer version)
4766+ Form_a_FunctionSpace_1(boost::shared_ptr<const dolfin::Restriction> restriction):
4767+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
4768+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4769+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()),
4770+ restriction)))
4771+ {
4772+ // Do nothing
4773+ }
4774+
4775+ // Copy constructor
4776+ ~Form_a_FunctionSpace_1()
4777+ {
4778+ }
4779+
4780+};
4781+
4782+typedef CoefficientSpace_mu Form_a_FunctionSpace_2;
4783+
4784+typedef CoefficientSpace_lmbda Form_a_FunctionSpace_3;
4785+
4786+class Form_a: public dolfin::Form
4787+{
4788+public:
4789+
4790+ // Constructor
4791+ Form_a(const dolfin::FunctionSpace& V1, const dolfin::FunctionSpace& V0):
4792+ dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1)
4793+ {
4794+ _function_spaces[0] = reference_to_no_delete_pointer(V0);
4795+ _function_spaces[1] = reference_to_no_delete_pointer(V1);
4796+
4797+ _ufc_form = boost::shared_ptr<const ufc::form>(new elasticity_form_0());
4798+ }
4799+
4800+ // Constructor
4801+ Form_a(const dolfin::FunctionSpace& V1, const dolfin::FunctionSpace& V0, const dolfin::GenericFunction& mu, const dolfin::GenericFunction& lmbda):
4802+ dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1)
4803+ {
4804+ _function_spaces[0] = reference_to_no_delete_pointer(V0);
4805+ _function_spaces[1] = reference_to_no_delete_pointer(V1);
4806+
4807+ this->mu = mu;
4808+ this->lmbda = lmbda;
4809+
4810+ _ufc_form = boost::shared_ptr<const ufc::form>(new elasticity_form_0());
4811+ }
4812+
4813+ // Constructor
4814+ Form_a(const dolfin::FunctionSpace& V1, const dolfin::FunctionSpace& V0, boost::shared_ptr<const dolfin::GenericFunction> mu, boost::shared_ptr<const dolfin::GenericFunction> lmbda):
4815+ dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1)
4816+ {
4817+ _function_spaces[0] = reference_to_no_delete_pointer(V0);
4818+ _function_spaces[1] = reference_to_no_delete_pointer(V1);
4819+
4820+ this->mu = *mu;
4821+ this->lmbda = *lmbda;
4822+
4823+ _ufc_form = boost::shared_ptr<const ufc::form>(new elasticity_form_0());
4824+ }
4825+
4826+ // Constructor
4827+ Form_a(boost::shared_ptr<const dolfin::FunctionSpace> V1, boost::shared_ptr<const dolfin::FunctionSpace> V0):
4828+ dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1)
4829+ {
4830+ _function_spaces[0] = V0;
4831+ _function_spaces[1] = V1;
4832+
4833+ _ufc_form = boost::shared_ptr<const ufc::form>(new elasticity_form_0());
4834+ }
4835+
4836+ // Constructor
4837+ Form_a(boost::shared_ptr<const dolfin::FunctionSpace> V1, boost::shared_ptr<const dolfin::FunctionSpace> V0, const dolfin::GenericFunction& mu, const dolfin::GenericFunction& lmbda):
4838+ dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1)
4839+ {
4840+ _function_spaces[0] = V0;
4841+ _function_spaces[1] = V1;
4842+
4843+ this->mu = mu;
4844+ this->lmbda = lmbda;
4845+
4846+ _ufc_form = boost::shared_ptr<const ufc::form>(new elasticity_form_0());
4847+ }
4848+
4849+ // Constructor
4850+ Form_a(boost::shared_ptr<const dolfin::FunctionSpace> V1, boost::shared_ptr<const dolfin::FunctionSpace> V0, boost::shared_ptr<const dolfin::GenericFunction> mu, boost::shared_ptr<const dolfin::GenericFunction> lmbda):
4851+ dolfin::Form(2, 2), mu(*this, 0), lmbda(*this, 1)
4852+ {
4853+ _function_spaces[0] = V0;
4854+ _function_spaces[1] = V1;
4855+
4856+ this->mu = *mu;
4857+ this->lmbda = *lmbda;
4858+
4859+ _ufc_form = boost::shared_ptr<const ufc::form>(new elasticity_form_0());
4860+ }
4861+
4862+ // Destructor
4863+ ~Form_a()
4864+ {}
4865+
4866+ /// Return the number of the coefficient with this name
4867+ virtual std::size_t coefficient_number(const std::string& name) const
4868+ {
4869+ if (name == "mu")
4870+ return 0;
4871+ else if (name == "lmbda")
4872+ return 1;
4873+
4874+ dolfin::dolfin_error("generated code for class Form",
4875+ "access coefficient data",
4876+ "Invalid coefficient");
4877+ return 0;
4878+ }
4879+
4880+ /// Return the name of the coefficient with this number
4881+ virtual std::string coefficient_name(std::size_t i) const
4882+ {
4883+ switch (i)
4884+ {
4885+ case 0:
4886+ return "mu";
4887+ case 1:
4888+ return "lmbda";
4889+ }
4890+
4891+ dolfin::dolfin_error("generated code for class Form",
4892+ "access coefficient data",
4893+ "Invalid coefficient");
4894+ return "unnamed";
4895+ }
4896+
4897+ // Typedefs
4898+ typedef Form_a_FunctionSpace_0 TestSpace;
4899+ typedef Form_a_FunctionSpace_1 TrialSpace;
4900+ typedef Form_a_FunctionSpace_2 CoefficientSpace_mu;
4901+ typedef Form_a_FunctionSpace_3 CoefficientSpace_lmbda;
4902+
4903+ // Coefficients
4904+ dolfin::CoefficientAssigner mu;
4905+ dolfin::CoefficientAssigner lmbda;
4906+};
4907+
4908+class Form_L_FunctionSpace_0: public dolfin::FunctionSpace
4909+{
4910+public:
4911+
4912+ //--- Constructors for standard function space, 2 different versions ---
4913+
4914+ // Create standard function space (reference version)
4915+ Form_L_FunctionSpace_0(const dolfin::Mesh& mesh):
4916+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4917+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4918+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), mesh)))
4919+ {
4920+ // Do nothing
4921+ }
4922+
4923+ // Create standard function space (shared pointer version)
4924+ Form_L_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh):
4925+ dolfin::FunctionSpace(mesh,
4926+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4927+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), *mesh)))
4928+ {
4929+ // Do nothing
4930+ }
4931+
4932+ //--- Constructors for constrained function space, 2 different versions ---
4933+
4934+ // Create standard function space (reference version)
4935+ Form_L_FunctionSpace_0(const dolfin::Mesh& mesh, const dolfin::SubDomain& constrained_domain):
4936+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(mesh),
4937+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4938+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), mesh,
4939+ dolfin::reference_to_no_delete_pointer(constrained_domain))))
4940+ {
4941+ // Do nothing
4942+ }
4943+
4944+ // Create standard function space (shared pointer version)
4945+ Form_L_FunctionSpace_0(boost::shared_ptr<const dolfin::Mesh> mesh, boost::shared_ptr<const dolfin::SubDomain> constrained_domain):
4946+ dolfin::FunctionSpace(mesh,
4947+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4948+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()), *mesh, constrained_domain)))
4949+ {
4950+ // Do nothing
4951+ }
4952+
4953+ //--- Constructors for restricted function space, 2 different versions ---
4954+
4955+ // Create restricted function space (reference version)
4956+ Form_L_FunctionSpace_0(const dolfin::Restriction& restriction):
4957+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction.mesh()),
4958+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4959+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()),
4960+ reference_to_no_delete_pointer(restriction))))
4961+ {
4962+ // Do nothing
4963+ }
4964+
4965+ // Create restricted function space (shared pointer version)
4966+ Form_L_FunctionSpace_0(boost::shared_ptr<const dolfin::Restriction> restriction):
4967+ dolfin::FunctionSpace(dolfin::reference_to_no_delete_pointer(restriction->mesh()),
4968+ boost::shared_ptr<const dolfin::FiniteElement>(new dolfin::FiniteElement(boost::shared_ptr<ufc::finite_element>(new elasticity_finite_element_2()))),
4969+ boost::shared_ptr<const dolfin::DofMap>(new dolfin::DofMap(boost::shared_ptr<ufc::dofmap>(new elasticity_dofmap_2()),
4970+ restriction)))
4971+ {
4972+ // Do nothing
4973+ }
4974+
4975+ // Copy constructor
4976+ ~Form_L_FunctionSpace_0()
4977+ {
4978+ }
4979+
4980+};
4981+
4982+typedef CoefficientSpace_f Form_L_FunctionSpace_1;
4983+
4984+class Form_L: public dolfin::Form
4985+{
4986+public:
4987+
4988+ // Constructor
4989+ Form_L(const dolfin::FunctionSpace& V0):
4990+ dolfin::Form(1, 1), f(*this, 0)
4991+ {
4992+ _function_spaces[0] = reference_to_no_delete_pointer(V0);
4993+
4994+ _ufc_form = boost::shared_ptr<const ufc::form>(new elasticity_form_1());
4995+ }
4996+
4997+ // Constructor
4998+ Form_L(const dolfin::FunctionSpace& V0, const dolfin::GenericFunction& f):
4999+ dolfin::Form(1, 1), f(*this, 0)
5000+ {
The diff has been truncated for viewing.