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