--- /dev/null
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.9)
+INCLUDE(../setup_testsubproject.cmake)
+PROJECT(testsuite CXX)
+IF(DEAL_II_WITH_SYMENGINE AND DEAL_II_WITH_CXX11)
+ DEAL_II_PICKUP_TESTS()
+ENDIF()
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that all fundamental operations with primitive SymEngine number
+// types work as expected
+
+// References:
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_basic.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_number.cpp
+
+#include "../tests.h"
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wextra-semi"
+#include <symengine/complex_double.h>
+#include <symengine/eval_double.h>
+#include <symengine/integer.h>
+#include <symengine/rational.h>
+#include <symengine/real_double.h>
+#include <symengine/visitor.h>
+#pragma clang diagnostic pop
+
+#include <fstream>
+#include <iomanip>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const int &val)
+{
+ return SE::integer(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const float &val)
+{
+ return SE::real_double(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const double &val)
+{
+ return SE::real_double(val);
+}
+
+template <typename NumberType>
+SE::RCP<const SE::Number>
+make_symengine_rcp(const std::complex<NumberType> &val)
+{
+ // Build complex from two SymEngine numbers
+ return SE::Complex::from_two_nums(*make_symengine_rcp(val.real()),
+ *make_symengine_rcp(val.imag()));
+}
+
+template <typename NumberType>
+void
+test_number()
+{
+ SE::RCP<const SE::Number> a = make_symengine_rcp(NumberType(4.2));
+ deallog << "a: " << *a << std::endl;
+
+ SE::RCP<const SE::Basic> b(make_symengine_rcp(NumberType(2.1)));
+ deallog << "b: " << *b << std::endl;
+
+ SE::RCP<const SE::Basic> c;
+ // deallog << "c (no constructor): " << *c << std::endl; // Fails
+
+ c = SE::add(a, b);
+ deallog << "c = a+b: " << *c << std::endl;
+
+ c = SE::sub(a, b);
+ deallog << "c = a-b: " << *c << std::endl;
+
+ c = SE::mul(a, b);
+ deallog << "c = a*b: " << *c << std::endl;
+
+ c = SE::div(a, b);
+ deallog << "c = a/b: " << *c << std::endl;
+
+ c = a;
+ deallog << "c = a: " << *c << std::endl;
+
+ c = SE::add(c, a);
+ deallog << "c += a: " << *c << std::endl;
+
+ c = SE::sub(c, a);
+ deallog << "c -= a: " << *c << std::endl;
+
+ c = SE::mul(c, a);
+ deallog << "c *= a: " << *c << std::endl;
+
+ c = SE::div(c, a);
+ deallog << "c /= a: " << *c << std::endl;
+}
+
+int
+main()
+{
+ initlog();
+
+ deallog.push("Integer");
+ test_number<int>();
+ deallog.pop();
+
+ deallog.push("Float");
+ test_number<float>();
+ deallog.pop();
+
+ deallog.push("Double");
+ test_number<double>();
+ deallog.pop();
+
+ // Not available yet
+ // SymEngine::SymEngineException: Invalid Format: Expected Integer or Rational
+ // deallog << "Complex double" << std::endl;
+ // test_number<std::complex<double>>();
+ // deallog.pop();
+
+ deallog << "OK" << std::endl;
+}
--- /dev/null
+
+DEAL:Integer::a: 4
+DEAL:Integer::b: 2
+DEAL:Integer::c = a+b: 6
+DEAL:Integer::c = a-b: 2
+DEAL:Integer::c = a*b: 8
+DEAL:Integer::c = a/b: 2
+DEAL:Integer::c = a: 4
+DEAL:Integer::c += a: 8
+DEAL:Integer::c -= a: 4
+DEAL:Integer::c *= a: 16
+DEAL:Integer::c /= a: 4
+DEAL:Float::a: 4.19999980926514
+DEAL:Float::b: 2.09999990463257
+DEAL:Float::c = a+b: 6.29999971389771
+DEAL:Float::c = a-b: 2.09999990463257
+DEAL:Float::c = a*b: 8.81999919891359
+DEAL:Float::c = a/b: 2.
+DEAL:Float::c = a: 4.19999980926514
+DEAL:Float::c += a: 8.39999961853027
+DEAL:Float::c -= a: 4.19999980926514
+DEAL:Float::c *= a: 17.6399983978272
+DEAL:Float::c /= a: 4.19999980926514
+DEAL:Double::a: 4.2
+DEAL:Double::b: 2.1
+DEAL:Double::c = a+b: 6.3
+DEAL:Double::c = a-b: 2.1
+DEAL:Double::c = a*b: 8.82
+DEAL:Double::c = a/b: 2.
+DEAL:Double::c = a: 4.2
+DEAL:Double::c += a: 8.4
+DEAL:Double::c -= a: 4.2
+DEAL:Double::c *= a: 17.64
+DEAL:Double::c /= a: 4.2
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that SymEngine can perform symbolic differentiation correctly
+
+// References:
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_basic.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_number.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_subs.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/symengine_casts.h
+// https://github.com/symengine/symengine/blob/master/symengine/symengine_rcp.h
+// https://github.com/symengine/symengine/blob/master/symengine/derivative.h
+// https://github.com/symengine/symengine/blob/master/symengine/subs.h
+
+#include "../tests.h"
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wextra-semi"
+#include <symengine/complex_double.h>
+#include <symengine/derivative.h>
+#include <symengine/eval_double.h>
+#include <symengine/integer.h>
+#include <symengine/rational.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+#include <symengine/symengine_casts.h>
+#include <symengine/visitor.h>
+#pragma clang diagnostic pop
+
+#include <fstream>
+#include <iomanip>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const int &val)
+{
+ return SE::integer(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const float &val)
+{
+ return SE::real_double(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const double &val)
+{
+ return SE::real_double(val);
+}
+
+template <typename NumberType>
+SE::RCP<const SE::Number>
+make_symengine_rcp(const std::complex<NumberType> &val)
+{
+ // Build complex from two SymEngine numbers
+ return SE::Complex::from_two_nums(*make_symengine_rcp(val.real()),
+ *make_symengine_rcp(val.imag()));
+}
+
+SE::RCP<const SE::Symbol>
+make_symengine_rcp(const std::string &name)
+{
+ return SE::symbol(name);
+}
+
+template <typename NumberType>
+void
+test_number()
+{
+ SE::RCP<const SE::Number> a = make_symengine_rcp(NumberType(3.1));
+ SE::RCP<const SE::Basic> b = make_symengine_rcp(NumberType(7.5));
+ deallog << "a: " << *a << std::endl;
+ deallog << "b: " << *b << std::endl;
+
+ SE::RCP<const SE::Symbol> x(make_symengine_rcp("x"));
+ SE::RCP<const SE::Basic> y(make_symengine_rcp("y"));
+ deallog << "x: " << *x << std::endl;
+ deallog << "y: " << *y << std::endl;
+
+ SE::RCP<const SE::Basic> c;
+ // deallog << "c (no constructor): " << *c << std::endl; // Fails
+
+ // Construction of symbolic function
+ c = SE::mul(y, SE::mul(SE::sub(y, b), SE::add(a, x)));
+ deallog << "c = y*(y-b)*(a+x): " << *c << std::endl;
+
+ // Perform symbolic differentiation
+ SE::RCP<const SE::Basic> dc_dx = c->diff(x);
+ // SE::RCP<const SE::Basic> dc_dy = diff(c,SE::implicit_cast<const
+ // SE::RCP<const SE::Symbol> &>(y));
+ SE::RCP<const SE::Basic> dc_dy =
+ diff(c, SE::rcp_static_cast<const SE::Symbol>(y));
+ SE::RCP<const SE::Basic> dc_dy_2 = sdiff(c, y);
+
+ deallog << "dc_dx = a*y*(y-b): " << *dc_dx << std::endl;
+ deallog << "dc_dy = 2*y*(a+x): " << *dc_dy << std::endl;
+ deallog << "dc_dy = 2*y*(a+x): " << *dc_dy_2 << std::endl;
+
+ // Substitute values
+ SE::map_basic_basic sub_vals;
+ sub_vals[x] = make_symengine_rcp(NumberType(1));
+ sub_vals[y] = make_symengine_rcp(NumberType(2.2));
+ deallog << "dc_dx(x=1,y=2.2): " << *(dc_dx->subs(sub_vals)) << std::endl;
+
+ sub_vals[x] = make_symengine_rcp(NumberType(0.75));
+ sub_vals[y] = make_symengine_rcp(NumberType(-2.08));
+ deallog << "dc_dy(x=0.75,y=-1.08): " << *(dc_dy->subs(sub_vals)) << std::endl;
+
+ deallog << std::endl;
+}
+
+int
+main()
+{
+ initlog();
+
+ deallog.push("Integer");
+ test_number<int>();
+ deallog.pop();
+
+ deallog.push("Float");
+ test_number<float>();
+ deallog.pop();
+
+ deallog.push("Double");
+ test_number<double>();
+ deallog.pop();
+
+ // Not available yet
+ // SymEngine::SymEngineException: Invalid Format: Expected Integer or Rational
+ // deallog << "Complex double" << std::endl;
+ // test_number<std::complex<double>>();
+ // deallog.pop();
+
+ deallog << "OK" << std::endl;
+}
--- /dev/null
+
+DEAL:Integer::a: 3
+DEAL:Integer::b: 7
+DEAL:Integer::x: x
+DEAL:Integer::y: y
+DEAL:Integer::c = y*(y-b)*(a+x): y*(-7 + y)*(3 + x)
+DEAL:Integer::dc_dx = a*y*(y-b): y*(-7 + y)
+DEAL:Integer::dc_dy = 2*y*(a+x): y*(3 + x) + (-7 + y)*(3 + x)
+DEAL:Integer::dc_dy = 2*y*(a+x): y*(3 + x) + (-7 + y)*(3 + x)
+DEAL:Integer::dc_dx(x=1,y=2.2): -10
+DEAL:Integer::dc_dy(x=0.75,y=-1.08): -33
+DEAL:Integer::
+DEAL:Float::a: 3.09999990463257
+DEAL:Float::b: 7.5
+DEAL:Float::x: x
+DEAL:Float::y: y
+DEAL:Float::c = y*(y-b)*(a+x): y*(3.09999990463257 + x)*(-7.5 + y)
+DEAL:Float::dc_dx = a*y*(y-b): y*(-7.5 + y)
+DEAL:Float::dc_dy = 2*y*(a+x): y*(3.09999990463257 + x) + (3.09999990463257 + x)*(-7.5 + y)
+DEAL:Float::dc_dy = 2*y*(a+x): y*(3.09999990463257 + x) + (3.09999990463257 + x)*(-7.5 + y)
+DEAL:Float::dc_dx(x=1,y=2.2): -11.6600001478195
+DEAL:Float::dc_dy(x=0.75,y=-1.08): -44.8909983005524
+DEAL:Float::
+DEAL:Double::a: 3.1
+DEAL:Double::b: 7.5
+DEAL:Double::x: x
+DEAL:Double::y: y
+DEAL:Double::c = y*(y-b)*(a+x): y*(3.1 + x)*(-7.5 + y)
+DEAL:Double::dc_dx = a*y*(y-b): y*(-7.5 + y)
+DEAL:Double::dc_dy = 2*y*(a+x): y*(3.1 + x) + (3.1 + x)*(-7.5 + y)
+DEAL:Double::dc_dy = 2*y*(a+x): y*(3.1 + x) + (3.1 + x)*(-7.5 + y)
+DEAL:Double::dc_dx(x=1,y=2.2): -11.66
+DEAL:Double::dc_dy(x=0.75,y=-1.08): -44.891
+DEAL:Double::
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that SymEngine can perform symbolic differentiation using function
+// symbols correctly
+
+// References:
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_functions.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/function.h
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_basic.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_number.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/tests/basic/test_subs.cpp
+// https://github.com/symengine/symengine/blob/master/symengine/symengine_casts.h
+// https://github.com/symengine/symengine/blob/master/symengine/symengine_rcp.h
+// https://github.com/symengine/symengine/blob/master/symengine/derivative.h
+// https://github.com/symengine/symengine/blob/master/symengine/subs.h
+
+#include "../tests.h"
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wextra-semi"
+#include <symengine/complex_double.h>
+#include <symengine/derivative.h>
+#include <symengine/eval_double.h>
+#include <symengine/integer.h>
+#include <symengine/rational.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+#include <symengine/symengine_casts.h>
+#include <symengine/visitor.h>
+#pragma clang diagnostic pop
+
+#include <fstream>
+#include <iomanip>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const int &val)
+{
+ return SE::integer(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const float &val)
+{
+ return SE::real_double(val);
+}
+
+SE::RCP<const SE::Number>
+make_symengine_rcp(const double &val)
+{
+ return SE::real_double(val);
+}
+
+template <typename NumberType>
+SE::RCP<const SE::Number>
+make_symengine_rcp(const std::complex<NumberType> &val)
+{
+ // Build complex from two SymEngine numbers
+ return SE::Complex::from_two_nums(*make_symengine_rcp(val.real()),
+ *make_symengine_rcp(val.imag()));
+}
+
+SE::RCP<const SE::Symbol>
+make_symengine_rcp(const std::string &name)
+{
+ return SE::symbol(name);
+}
+
+SE::RCP<const SE::Basic>
+make_symengine_rcp(const std::string & name,
+ const std::vector<SE::RCP<const SE::Basic>> &args)
+{
+ return SE::function_symbol(name, args);
+}
+
+SE::RCP<const SE::Basic>
+make_symengine_rcp(const std::string &name, const SE::RCP<const SE::Basic> &arg)
+{
+ return SE::function_symbol(name, arg);
+}
+
+template <typename NumberType>
+void
+test_number()
+{
+ // Modified snippet from SymEngine test
+ // {
+ // SE::RCP<const SE::Symbol> x = SE::symbol("x");
+ // SE::RCP<const SE::Symbol> _xi_1 = SE::symbol("_xi_1");
+ // SE::RCP<const SE::Basic> f = SE::function_symbol("f", x);
+ // SE::RCP<const SE::Basic> r1, r2;
+ //
+ // f = function_symbol("f", SE::pow(x, SE::integer(2)));
+ // deallog << "f: " << *f << std::endl;
+ // r1 = f->diff(x);
+ // deallog << "r1: " << *r1 << std::endl;
+ // r2 = SE::Derivative::create(SE::function_symbol("f", _xi_1), {_xi_1});
+ // deallog << "r2: " << *r2 << std::endl;
+ // r2 = SE::Subs::create(r2, {{_xi_1, SE::pow(x, SE::integer(2))}});
+ // deallog << "r2: " << *r2 << std::endl;
+ // }
+
+ // Normal symbols
+ SE::RCP<const SE::Symbol> x(make_symengine_rcp("x"));
+ SE::RCP<const SE::Basic> y(make_symengine_rcp("y"));
+ deallog << "x: " << *x << std::endl;
+ deallog << "y: " << *y << std::endl;
+
+ // Function symbols
+ SE::RCP<const SE::Basic> f(make_symengine_rcp("f", x));
+ SE::RCP<const SE::Basic> g(make_symengine_rcp("g", {x, y}));
+ SE::RCP<const SE::Basic> h(
+ make_symengine_rcp("h", {SE::add(x, y), SE::mul(x, y)}));
+ deallog << "f: " << *f << std::endl;
+ deallog << "g: " << *g << std::endl;
+ deallog << "h: " << *h << std::endl;
+
+ // Perform symbolic differentiation
+ SE::RCP<const SE::Basic> df_dx = f->diff(x);
+ SE::RCP<const SE::Basic> dg_dx = g->diff(x);
+ SE::RCP<const SE::Basic> dh_dx = h->diff(x);
+ deallog << "df_dx: " << *df_dx << std::endl;
+ deallog << "dg_dx: " << *dg_dx << std::endl;
+ deallog << "dh_dx: " << *dh_dx << std::endl;
+ SE::RCP<const SE::Basic> df_dy =
+ f->diff(SE::rcp_static_cast<const SE::Symbol>(y));
+ SE::RCP<const SE::Basic> dg_dy =
+ g->diff(SE::rcp_static_cast<const SE::Symbol>(y));
+ SE::RCP<const SE::Basic> dh_dy =
+ h->diff(SE::rcp_static_cast<const SE::Symbol>(y));
+ deallog << "df_dy: " << *df_dy << std::endl;
+ deallog << "dg_dy: " << *dg_dy << std::endl;
+ deallog << "dh_dy: " << *dh_dy << std::endl;
+
+ // Substitute values
+ SE::map_basic_basic sub_vals;
+ sub_vals[x] = make_symengine_rcp(NumberType(1));
+ sub_vals[y] = make_symengine_rcp(NumberType(2.2));
+ sub_vals[f] = SE::pow(x, SE::integer(2));
+ sub_vals[g] = SE::add(x, y);
+ sub_vals[h] = SE::sub(x, y);
+
+ deallog << "f(x=1,y=2.2): " << *(f->subs(sub_vals))
+ << std::endl; // ->subs(sub_vals)
+ deallog << "g(x=1,y=2.2): " << *(g->subs(sub_vals))
+ << std::endl; // ->subs(sub_vals)
+ deallog << "h(x=1,y=2.2): " << *(h->subs(sub_vals))
+ << std::endl; // ->subs(sub_vals)
+ deallog << "Eval: f(x=1,y=2.2): "
+ << eval_double(*(f->subs(sub_vals)->subs(sub_vals))) << std::endl;
+ deallog << "Eval: g(x=1,y=2.2): "
+ << eval_double(*(g->subs(sub_vals)->subs(sub_vals))) << std::endl;
+ deallog << "Eval: h(x=1,y=2.2): "
+ << eval_double(*(h->subs(sub_vals)->subs(sub_vals))) << std::endl;
+
+ // Not yet implemented
+ deallog << "df_dx(x=1,y=2.2): " << *(df_dx->subs(sub_vals))
+ << std::endl; // ->subs(sub_vals)
+ deallog << "dg_dx(x=1,y=2.2): " << *(dg_dx->subs(sub_vals))
+ << std::endl; // ->subs(sub_vals)
+ deallog << "dh_dx(x=1,y=2.2): " << *(dh_dx->subs(sub_vals))
+ << std::endl; // ->subs(sub_vals)
+ // deallog << "Eval: df_dx(x=1,y=2.2): " <<
+ // eval_double(*(df_dx->subs(sub_vals)->subs(sub_vals))) << std::endl; deallog
+ // << "Eval: dg_dx(x=1,y=2.2): " <<
+ // eval_double(*(dg_dx->subs(sub_vals)->subs(sub_vals))) << std::endl; deallog
+ // << "Eval: dh_dx(x=1,y=2.2): " <<
+ // eval_double(*(dh_dx->subs(sub_vals)->subs(sub_vals))) << std::endl;
+}
+
+int
+main()
+{
+ initlog();
+
+ deallog.push("Integer");
+ test_number<int>();
+ deallog.pop();
+
+ deallog.push("Float");
+ test_number<float>();
+ deallog.pop();
+
+ deallog.push("Double");
+ test_number<double>();
+ deallog.pop();
+
+ // Not available yet
+ // SymEngine::SymEngineException: Invalid Format: Expected Integer or Rational
+ // deallog << "Complex double" << std::endl;
+ // test_number<std::complex<double>>();
+ // deallog.pop();
+
+ deallog << "OK" << std::endl;
+}
--- /dev/null
+
+DEAL:Integer::x: x
+DEAL:Integer::y: y
+DEAL:Integer::f: f(x)
+DEAL:Integer::g: g(x, y)
+DEAL:Integer::h: h(x + y, x*y)
+DEAL:Integer::df_dx: Derivative(f(x), x)
+DEAL:Integer::dg_dx: Derivative(g(x, y), x)
+DEAL:Integer::dh_dx: y*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Integer::df_dy: 0
+DEAL:Integer::dg_dy: Derivative(g(x, y), y)
+DEAL:Integer::dh_dy: x*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Integer::f(x=1,y=2.2): x**2
+DEAL:Integer::g(x=1,y=2.2): x + y
+DEAL:Integer::h(x=1,y=2.2): x - y
+DEAL:Integer::Eval: f(x=1,y=2.2): 1.00
+DEAL:Integer::Eval: g(x=1,y=2.2): 3.00
+DEAL:Integer::Eval: h(x=1,y=2.2): -1.00
+DEAL:Integer::df_dx(x=1,y=2.2): 2*x
+DEAL:Integer::dg_dx(x=1,y=2.2): 1
+DEAL:Integer::dh_dx(x=1,y=2.2): 2*Subs(Derivative(h(3, _xi_2), _xi_2), (_xi_2), (2)) + Subs(Derivative(h(_xi_1, 2), _xi_1), (_xi_1), (3))
+DEAL:Float::x: x
+DEAL:Float::y: y
+DEAL:Float::f: f(x)
+DEAL:Float::g: g(x, y)
+DEAL:Float::h: h(x + y, x*y)
+DEAL:Float::df_dx: Derivative(f(x), x)
+DEAL:Float::dg_dx: Derivative(g(x, y), x)
+DEAL:Float::dh_dx: y*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Float::df_dy: 0
+DEAL:Float::dg_dy: Derivative(g(x, y), y)
+DEAL:Float::dh_dy: x*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Float::f(x=1,y=2.2): x**2
+DEAL:Float::g(x=1,y=2.2): x + y
+DEAL:Float::h(x=1,y=2.2): x - y
+DEAL:Float::Eval: f(x=1,y=2.2): 1.00
+DEAL:Float::Eval: g(x=1,y=2.2): 3.20
+DEAL:Float::Eval: h(x=1,y=2.2): -1.20
+DEAL:Float::df_dx(x=1,y=2.2): 2*x
+DEAL:Float::dg_dx(x=1,y=2.2): 1
+DEAL:Float::dh_dx(x=1,y=2.2): 2.20000004768372*Subs(Derivative(h(3.20000004768372, _xi_2), _xi_2), (_xi_2), (2.20000004768372)) + Subs(Derivative(h(_xi_1, 2.20000004768372), _xi_1), (_xi_1), (3.20000004768372))
+DEAL:Double::x: x
+DEAL:Double::y: y
+DEAL:Double::f: f(x)
+DEAL:Double::g: g(x, y)
+DEAL:Double::h: h(x + y, x*y)
+DEAL:Double::df_dx: Derivative(f(x), x)
+DEAL:Double::dg_dx: Derivative(g(x, y), x)
+DEAL:Double::dh_dx: y*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Double::df_dy: 0
+DEAL:Double::dg_dy: Derivative(g(x, y), y)
+DEAL:Double::dh_dy: x*Subs(Derivative(h(x + y, _xi_2), _xi_2), (_xi_2), (x*y)) + Subs(Derivative(h(_xi_1, x*y), _xi_1), (_xi_1), (x + y))
+DEAL:Double::f(x=1,y=2.2): x**2
+DEAL:Double::g(x=1,y=2.2): x + y
+DEAL:Double::h(x=1,y=2.2): x - y
+DEAL:Double::Eval: f(x=1,y=2.2): 1.00
+DEAL:Double::Eval: g(x=1,y=2.2): 3.20
+DEAL:Double::Eval: h(x=1,y=2.2): -1.20
+DEAL:Double::df_dx(x=1,y=2.2): 2*x
+DEAL:Double::dg_dx(x=1,y=2.2): 1
+DEAL:Double::dh_dx(x=1,y=2.2): 2.2*Subs(Derivative(h(3.2, _xi_2), _xi_2), (_xi_2), (2.2)) + Subs(Derivative(h(_xi_1, 2.2), _xi_1), (_xi_1), (3.2))
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that SymEngine can do some optimsation using lambda functions
+
+// References:
+// https://github.com/symengine/symengine/blob/master/symengine/tests/eval/test_lambda_double.cpp
+
+#include "../tests.h"
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wextra-semi"
+
+#include <symengine/add.h>
+#include <symengine/basic.h>
+#include <symengine/constants.h>
+#include <symengine/dict.h>
+#include <symengine/integer.h>
+#include <symengine/lambda_double.h>
+#include <symengine/llvm_double.h>
+#include <symengine/mul.h>
+#include <symengine/parser.h>
+#include <symengine/pow.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+
+#pragma GCC diagnostic pop
+
+#include <chrono>
+#include <iostream>
+#include <vector>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+int
+main(int argc, char *argv[])
+{
+ initlog();
+
+ const unsigned int n_runs = 1000;
+
+ SE::RCP<const SE::Symbol> t_11 = SE::symbol("t_11");
+ SE::RCP<const SE::Symbol> t_12 = SE::symbol("t_12");
+ SE::RCP<const SE::Symbol> s = SE::symbol("s");
+ SE::RCP<const SE::Symbol> t_02 = SE::symbol("t_02");
+ SE::RCP<const SE::Symbol> v_0 = SE::symbol("v_0");
+ SE::RCP<const SE::Symbol> v_1 = SE::symbol("v_1");
+ SE::RCP<const SE::Symbol> v_2 = SE::symbol("v_2");
+ SE::RCP<const SE::Symbol> t_10 = SE::symbol("t_10");
+ SE::RCP<const SE::Symbol> t_00 = SE::symbol("t_00");
+ SE::RCP<const SE::Symbol> t_01 = SE::symbol("t_01");
+ SE::RCP<const SE::Symbol> t_20 = SE::symbol("t_20");
+ SE::RCP<const SE::Symbol> t_21 = SE::symbol("t_21");
+ SE::RCP<const SE::Symbol> t_22 = SE::symbol("t_22");
+
+ SE::vec_basic v = {
+ t_11, t_12, s, t_02, v_0, v_1, v_2, t_10, t_00, t_01, t_20, t_21, t_22};
+
+ SE::RCP<const SE::Basic> h =
+ SE::parse("2.*s**2.2*(v_0**2 + v_1**2 + v_2**2)**3*"
+ "(-1.*t_20*(1.*t_01*t_12 - 1.*t_02*t_11) "
+ "+ 1.*t_21*(1.*t_00*t_12 - 1.*t_02*t_10) "
+ "- 1.*t_22*(1.*t_00*t_11 - 1.*t_01*t_10))"
+ "*(1.*t_21*t_12 - 1.*t_22*t_11)");
+
+ SE::map_basic_basic dict;
+ std::vector<double> vals(v.size());
+ for (unsigned i = 0; i < v.size(); i++)
+ {
+ dict[v[i]] = SE::real_double(i);
+ vals[i] = i;
+ }
+
+ double res = 0.0, res1 = 0.0, res2 = 0.0;
+
+
+ // Standard substitution
+
+ auto t1 = std::chrono::high_resolution_clock::now();
+ for (unsigned j = 0; j < n_runs; j++)
+ {
+ res += static_cast<const SE::RealDouble &>(*h->subs(dict)).i;
+ }
+ auto t2 = std::chrono::high_resolution_clock::now();
+ const double diff_symm_subs =
+ std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+ std::cout
+ << "Subs " << n_runs << " calls :"
+ << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+ << " us" << std::endl;
+
+ deallog << "Value (subs): " << res << std::endl;
+
+ // Optimisation via lambda functions
+
+ t1 = std::chrono::high_resolution_clock::now();
+ SE::LambdaRealDoubleVisitor l;
+ l.init(v, *h);
+ t2 = std::chrono::high_resolution_clock::now();
+ const double diff_lambda_real_double_setup =
+ std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+ std::cout
+ << "LambdaDoubleVisitor setup-time :"
+ << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+ << " us" << std::endl;
+
+ t1 = std::chrono::high_resolution_clock::now();
+ for (unsigned j = 0; j < n_runs; j++)
+ {
+ res1 += l.call(vals);
+ }
+ t2 = std::chrono::high_resolution_clock::now();
+ const double diff_lambda_real_double_call =
+ std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+ std::cout
+ << "LambdaDoubleVisitor run-time :"
+ << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+ << " us" << std::endl;
+
+ deallog << "Value (lambda): " << res1 << std::endl;
+
+ std::cout << "Timing ratio (subs/lambda_double): "
+ << (diff_symm_subs) /
+ (diff_lambda_real_double_setup + diff_lambda_real_double_call)
+ << std::endl;
+
+#ifdef HAVE_SYMENGINE_LLVM
+ // Optimisation via LLVM JIT optimiser
+
+ t1 = std::chrono::high_resolution_clock::now();
+ SE::LLVMDoubleVisitor l2;
+ l2.init(v, *h);
+ t2 = std::chrono::high_resolution_clock::now();
+ std::cout
+ << "LLVMDoubleVisitor setup-time :"
+ << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+ << " us" << std::endl;
+
+ t1 = std::chrono::high_resolution_clock::now();
+ for (unsigned j = 0; j < n_runs; j++)
+ {
+ res2 += l2.call(vals);
+ }
+ t2 = std::chrono::high_resolution_clock::now();
+ std::cout
+ << "LLVMDoubleVisitor run-time :"
+ << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+ << " us : Value: " << res2 << std::endl;
+
+ deallog << "Value (llvm): " << res2 << std::endl;
+#endif
+
+ deallog << "OK" << std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::Value (subs): 2.41358e+13
+DEAL::Value (lambda): 2.41358e+13
+DEAL::OK
--- /dev/null
+
+DEAL::Value (subs): 2.41358e+13
+DEAL::Value (lambda): 2.41358e+13
+DEAL::Value (llvm): 2.41358e+13
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check that SymEngine can do some optimsation on diffentiated expressions
+// using lambda functions
+// This is an extension of basic_04.cc
+
+// References:
+// https://github.com/symengine/symengine/blob/master/symengine/tests/eval/test_lambda_double.cpp
+
+#include "../tests.h"
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wextra-semi"
+
+#include <symengine/add.h>
+#include <symengine/basic.h>
+#include <symengine/dict.h>
+#include <symengine/integer.h>
+#include <symengine/lambda_double.h>
+#include <symengine/llvm_double.h>
+#include <symengine/mul.h>
+#include <symengine/parser.h>
+#include <symengine/pow.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+
+#pragma GCC diagnostic pop
+
+#include <chrono>
+#include <iostream>
+#include <vector>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+int
+main(int argc, char *argv[])
+{
+ initlog();
+
+ const unsigned int n_runs = 1000;
+
+ SE::RCP<const SE::Symbol> t_11 = SE::symbol("t_11");
+ SE::RCP<const SE::Symbol> t_12 = SE::symbol("t_12");
+ SE::RCP<const SE::Symbol> s = SE::symbol("s");
+ SE::RCP<const SE::Symbol> t_02 = SE::symbol("t_02");
+ SE::RCP<const SE::Symbol> v_0 = SE::symbol("v_0");
+ SE::RCP<const SE::Symbol> v_1 = SE::symbol("v_1");
+ SE::RCP<const SE::Symbol> v_2 = SE::symbol("v_2");
+ SE::RCP<const SE::Symbol> t_10 = SE::symbol("t_10");
+ SE::RCP<const SE::Symbol> t_00 = SE::symbol("t_00");
+ SE::RCP<const SE::Symbol> t_01 = SE::symbol("t_01");
+ SE::RCP<const SE::Symbol> t_20 = SE::symbol("t_20");
+ SE::RCP<const SE::Symbol> t_21 = SE::symbol("t_21");
+ SE::RCP<const SE::Symbol> t_22 = SE::symbol("t_22");
+
+ SE::vec_basic v = {
+ t_11, t_12, s, t_02, v_0, v_1, v_2, t_10, t_00, t_01, t_20, t_21, t_22};
+
+ SE::RCP<const SE::Basic> h =
+ SE::parse("2.*s**2.2*(v_0**2 + v_1**2 + v_2**2)**3*"
+ "(-1.*t_20*(1.*t_01*t_12 - 1.*t_02*t_11) "
+ "+ 1.*t_21*(1.*t_00*t_12 - 1.*t_02*t_10) "
+ "- 1.*t_22*(1.*t_00*t_11 - 1.*t_01*t_10))"
+ "*(1.*t_21*t_12 - 1.*t_22*t_11)");
+
+ SE::vec_basic diffs = {h->diff(t_20),
+ h->diff(t_21),
+ h->diff(t_22),
+ h->diff(t_10),
+ h->diff(t_11),
+ h->diff(t_12),
+ h->diff(t_00),
+ h->diff(t_01),
+ h->diff(t_02),
+ h->diff(s)};
+
+ std::cout << *h << std::endl;
+
+ SE::map_basic_basic dict;
+ std::vector<double> vals(v.size());
+ for (unsigned i = 0; i < v.size(); i++)
+ {
+ dict[v[i]] = SE::real_double(i);
+ vals[i] = i;
+ }
+
+ double res[10] = {}, res1[10] = {}, res2[10] = {};
+ double r = 0.0, r2 = 0.0, r3 = 0.0;
+
+ // Standard substitution
+
+ auto t1 = std::chrono::high_resolution_clock::now();
+ for (unsigned j = 0; j < n_runs; j++)
+ {
+ for (unsigned k = 0; k < diffs.size(); k++)
+ {
+ res[k] = static_cast<const SE::RealDouble &>(*diffs[k]->subs(dict)).i;
+ r += res[k];
+ }
+ }
+ auto t2 = std::chrono::high_resolution_clock::now();
+ const double diff_symm_subs =
+ std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+ std::cout
+ << "Subs " << n_runs << " calls :"
+ << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+ << " us" << std::endl;
+
+ deallog << "Value (subs): " << r << std::endl;
+
+ // Optimisation via lambda functions
+
+ t1 = std::chrono::high_resolution_clock::now();
+ SE::LambdaRealDoubleVisitor l;
+ l.init(v, diffs);
+ t2 = std::chrono::high_resolution_clock::now();
+ const double diff_lambda_real_double_setup =
+ std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+ std::cout
+ << "LambdaDoubleVisitor setup-time :"
+ << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+ << " us" << std::endl;
+
+ t1 = std::chrono::high_resolution_clock::now();
+ for (unsigned j = 0; j < n_runs; j++)
+ {
+ l.call(res1, vals.data());
+ for (unsigned k = 0; k < diffs.size(); k++)
+ {
+ r2 += res1[k];
+ }
+ }
+ t2 = std::chrono::high_resolution_clock::now();
+ const double diff_lambda_real_double_call =
+ std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
+ std::cout
+ << "LambdaDoubleVisitor run-time :"
+ << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+ << " us" << std::endl;
+
+ deallog << "Value (lambda): " << r2 << std::endl;
+
+#ifdef HAVE_SYMENGINE_LLVM
+ // Optimisation via LLVM JIT optimiser
+
+ t1 = std::chrono::high_resolution_clock::now();
+ SE::LLVMDoubleVisitor l2;
+ l2.init(v, diffs);
+ t2 = std::chrono::high_resolution_clock::now();
+ std::cout
+ << "LLVMDoubleVisitor setup-time :"
+ << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+ << " us" << std::endl;
+
+ t1 = std::chrono::high_resolution_clock::now();
+ for (unsigned j = 0; j < 1000; j++)
+ {
+ l2.call(res2, vals.data());
+ for (unsigned k = 0; k < diffs.size(); k++)
+ {
+ r3 += res2[k];
+ }
+ }
+ t2 = std::chrono::high_resolution_clock::now();
+ std::cout
+ << "LLVMDoubleVisitor run-time :"
+ << std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count()
+ << " us : Value: " << r3 << std::endl;
+
+ deallog << "Value (llvm): " << r3 << std::endl;
+#endif
+
+ deallog << "OK" << std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::Value (subs): 2.91338e+13
+DEAL::Value (lambda): 2.91338e+13
+DEAL::OK
--- /dev/null
+
+DEAL::Value (subs): 2.91338e+13
+DEAL::Value (lambda): 2.91338e+13
+DEAL::Value (llvm): 2.91338e+13
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check that SymEngine can do delayed substitution of explicitly dependent
+// symbolic functions and evaluate their derivatives as well.
+// This is inline with what is required to implement some material laws that
+// have internal variables
+
+#include "../tests.h"
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wextra-semi"
+
+#include <symengine/add.h>
+#include <symengine/basic.h>
+#include <symengine/dict.h>
+#include <symengine/eval_double.h>
+#include <symengine/functions.h>
+#include <symengine/integer.h>
+#include <symengine/mul.h>
+#include <symengine/pow.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+
+#pragma GCC diagnostic pop
+
+#include <iostream>
+#include <vector>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+int
+main(int argc, char *argv[])
+{
+ initlog();
+
+ // 0. Define some independent variables
+ // Here we have two definitions of Q: One it is treated as a completely
+ // indepdent variable, and the other where its dependence on C is
+ // explicitly defined
+ SE::RCP<const SE::Symbol> C = SE::symbol("C"); // Primary variable
+ SE::RCP<const SE::Symbol> Qi = SE::symbol("Qi"); // Some internal variable
+ SE::RCP<const SE::Basic> Q =
+ SE::function_symbol("Q", C); // This sets up the dependence of Q on C
+
+ deallog << *C << ", " << *Qi << ", " << *Q << std::endl;
+
+ // 1. Define a function with Q being independent of C
+ SE::RCP<const SE::Basic> f_CQ_symb =
+ SE::mul(SE::real_double(0.5), SE::mul(Qi, SE::pow(C, SE::integer(2))));
+
+ deallog << "f_CQ_symb: " << *f_CQ_symb << std::endl;
+
+ // 2. Compute the partial derivative of f wrt C. This is a "partial
+ // derivative" as Q != Q(C).
+ SE::RCP<const SE::Basic> df_CQ_dC_symb = f_CQ_symb->diff(C);
+
+ deallog << "df_CQ_dC_symb: " << *df_CQ_dC_symb << std::endl;
+
+ // 3. Substitute Qi -> Q=Q(C) to now make SymEngine aware of Q's dependence
+ // on C.
+ SE::map_basic_basic int_var_dict;
+ int_var_dict[Qi] = Q;
+ SE::RCP<const SE::Basic> df_CQ_dC_symb_subs =
+ df_CQ_dC_symb->subs(int_var_dict);
+
+ deallog << "df_CQ_dC_symb_subs: " << *df_CQ_dC_symb_subs << std::endl;
+
+ // 4. Compute total derivative of df_dC wrt C.
+ SE::RCP<const SE::Basic> D2f_CQ_DC_dC_symb = df_CQ_dC_symb_subs->diff(C);
+
+ deallog << "D2f_CQ_DC_dC_symb: " << *D2f_CQ_DC_dC_symb << std::endl;
+
+ // 5. Perform some numerical substitutions
+ SE::map_basic_basic all_var_dict;
+ all_var_dict[C] = SE::real_double(5);
+ all_var_dict[Qi] = SE::real_double(3);
+ all_var_dict[Q] =
+ SE::mul(SE::real_double(0.6),
+ C); // This explicitly defines the relationship between Q and C
+
+ SE::RCP<const SE::Basic> f_CQ_subs = f_CQ_symb->subs(all_var_dict);
+ deallog << "f_CQ_subs: " << *f_CQ_subs << std::endl;
+
+ SE::RCP<const SE::Basic> df_CQ_dC_subs = df_CQ_dC_symb->subs(all_var_dict);
+ deallog << "df_CQ_dC: subs: " << *df_CQ_dC_subs << std::endl
+ << "df_CQ_dC: eval: " << SE::eval_double(*df_CQ_dC_subs) << std::endl;
+
+ // This first substitution should presumably convert Q(C)->0.6C and thus
+ // dQ(C)_dC into 0.6
+ SE::RCP<const SE::Basic> D2f_CQ_DC_dC_subs_1 =
+ D2f_CQ_DC_dC_symb->subs(all_var_dict);
+ // This second substitution should fill out the remaining entries for C
+ SE::RCP<const SE::Basic> D2f_CQ_DC_dC_subs =
+ D2f_CQ_DC_dC_subs_1->subs(all_var_dict);
+ deallog << "D2f_CQ_DC_dC: subs 1: " << *D2f_CQ_DC_dC_subs_1 << std::endl
+ << "D2f_CQ_DC_dC: subs 2: " << *D2f_CQ_DC_dC_subs << std::endl
+ << "D2f_CQ_DC_dC: eval: " << SE::eval_double(*D2f_CQ_DC_dC_subs)
+ << std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::C, Qi, Q(C)
+DEAL::f_CQ_symb: 0.5*C**2*Qi
+DEAL::df_CQ_dC_symb: 1.0*C*Qi
+DEAL::df_CQ_dC_symb_subs: 1.0*C*Q(C)
+DEAL::D2f_CQ_DC_dC_symb: 1.0*C*Derivative(Q(C), C) + 1.0*Q(C)
+DEAL::f_CQ_subs: 37.5
+DEAL::df_CQ_dC: subs: 15.
+DEAL::df_CQ_dC: eval: 15.0000
+DEAL::D2f_CQ_DC_dC: subs 1: 3.0 + 0.6*C
+DEAL::D2f_CQ_DC_dC: subs 2: 6.0
+DEAL::D2f_CQ_DC_dC: eval: 6.00000
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check that SymEngine can do delayed substitution of implicitly dependent
+// symbolic functions and evaluate their derivatives as well.
+// This is inline with what is required to implement some material laws that
+// have internal variables
+
+#include "../tests.h"
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wextra-semi"
+
+#include <symengine/add.h>
+#include <symengine/basic.h>
+#include <symengine/dict.h>
+#include <symengine/eval_double.h>
+#include <symengine/functions.h>
+#include <symengine/integer.h>
+#include <symengine/mul.h>
+#include <symengine/pow.h>
+#include <symengine/real_double.h>
+#include <symengine/symbol.h>
+
+#pragma GCC diagnostic pop
+
+#include <iostream>
+#include <vector>
+
+using namespace dealii;
+namespace SE = SymEngine;
+
+int
+main(int argc, char *argv[])
+{
+ initlog();
+
+ // 0. Define some independent variables
+ // Here we have two definitions of Q: One it is treated as a completely
+ // independent variable, and the other where its dependence on C is
+ // implicitly defined
+ // That is to say, we'd have to compute Q=Q(C) using a nonlinear solution
+ // scheme, as might be the case for a viscoelastic material with a
+ // nonlinear evolution law.
+ SE::RCP<const SE::Symbol> C = SE::symbol("C"); // Primary variable
+ SE::RCP<const SE::Symbol> Qi = SE::symbol("Qi"); // Some internal variable
+ SE::RCP<const SE::Basic> Q =
+ SE::function_symbol("Q", C); // This sets up the dependence of Q on C
+
+ deallog << *C << ", " << *Qi << ", " << *Q << std::endl;
+
+ // 1. Define a function with Q being independent of C
+ SE::RCP<const SE::Basic> f_CQ_symb =
+ SE::mul(SE::real_double(0.5), SE::mul(Qi, SE::pow(C, SE::integer(2))));
+
+ deallog << "f_CQ_symb: " << *f_CQ_symb << std::endl;
+
+ // 2. Compute the partial derivative of f wrt C. This is a "partial
+ // derivative" as Q != Q(C).
+ SE::RCP<const SE::Basic> df_CQ_dC_symb = f_CQ_symb->diff(C);
+
+ deallog << "df_CQ_dC_symb: " << *df_CQ_dC_symb << std::endl;
+
+ // 3. Substitute Qi -> Q=Q(C) to now make SymEngine aware of Q's dependence
+ // on C.
+ SE::map_basic_basic int_var_dict;
+ int_var_dict[Qi] = Q;
+ SE::RCP<const SE::Basic> df_CQ_dC_symb_subs =
+ df_CQ_dC_symb->subs(int_var_dict);
+
+ deallog << "df_CQ_dC_symb_subs: " << *df_CQ_dC_symb_subs << std::endl;
+
+ // 4. Compute total derivative of df_dC wrt C.
+ SE::RCP<const SE::Basic> D2f_CQ_DC_dC_symb = df_CQ_dC_symb_subs->diff(C);
+
+ deallog << "D2f_CQ_DC_dC_symb: " << *D2f_CQ_DC_dC_symb << std::endl;
+
+ // 5. Perform some numerical substitutions
+ // Lets assume Q = 0.6*C*C
+ // Then dQ_dC = 1.2*C
+ SE::map_basic_basic all_var_dict;
+ all_var_dict[C] = SE::real_double(5);
+ all_var_dict[Qi] =
+ SE::real_double(0.6 * 5 * 5); // This would be the numerical end result of
+ // an internal nonlinear solver for Q
+ all_var_dict[Q] = Qi; // These are numerically the same thing
+ all_var_dict[Q->diff(C)] =
+ SE::mul(SE::real_double(1.2),
+ C); // This resolves the implicit relationship between Q and C
+
+ SE::RCP<const SE::Basic> f_CQ_subs = f_CQ_symb->subs(all_var_dict);
+ deallog << "f_CQ_subs: " << *f_CQ_subs << std::endl;
+
+ SE::RCP<const SE::Basic> df_CQ_dC_subs = df_CQ_dC_symb->subs(all_var_dict);
+ deallog << "df_CQ_dC: subs: " << *df_CQ_dC_subs << std::endl
+ << "df_CQ_dC: eval: " << SE::eval_double(*df_CQ_dC_subs) << std::endl;
+
+ // This first substitution should presumably convert Q(C)->0.6C and thus
+ // dQ(C)_dC into 0.6
+ SE::RCP<const SE::Basic> D2f_CQ_DC_dC_subs_1 =
+ D2f_CQ_DC_dC_symb->subs(all_var_dict);
+ // This second substitution should fill out the remaining entries for C
+ SE::RCP<const SE::Basic> D2f_CQ_DC_dC_subs =
+ D2f_CQ_DC_dC_subs_1->subs(all_var_dict);
+ deallog << "D2f_CQ_DC_dC: subs 1: " << *D2f_CQ_DC_dC_subs_1 << std::endl
+ << "D2f_CQ_DC_dC: subs 2: " << *D2f_CQ_DC_dC_subs << std::endl
+ << "D2f_CQ_DC_dC: eval: " << SE::eval_double(*D2f_CQ_DC_dC_subs)
+ << std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::C, Qi, Q(C)
+DEAL::f_CQ_symb: 0.5*C**2*Qi
+DEAL::df_CQ_dC_symb: 1.0*C*Qi
+DEAL::df_CQ_dC_symb_subs: 1.0*C*Q(C)
+DEAL::D2f_CQ_DC_dC_symb: 1.0*C*Derivative(Q(C), C) + 1.0*Q(C)
+DEAL::f_CQ_subs: 187.5
+DEAL::df_CQ_dC: subs: 75.
+DEAL::df_CQ_dC: eval: 75.0000
+DEAL::D2f_CQ_DC_dC: subs 1: 6.0*C + 1.0*Qi
+DEAL::D2f_CQ_DC_dC: subs 2: 45.0
+DEAL::D2f_CQ_DC_dC: eval: 45.0000
projects / dealii.git / commitdiff