python-ffc 2016.2.0-1 / usr / lib / python2.7 / dist-packages / ffc / backends / dolfin / goalfunctional.py

# -*- coding: utf-8 -*-
# Copyright (C) 2010 Marie E. Rognes
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can 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 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
#
# Last changed: 2011-07-06

__all__ = ["generate_update_ec"]


attach_coefficient_template = """
    // Attach coefficients from %(from)s to %(to)s
    for (std::size_t i = 0; i < %(from)s.num_coefficients(); i++)
    {
      name = %(from)s.coefficient_name(i);
      // Don't attach discrete primal solution here (not computed).
      if (name == "__discrete_primal_solution")
        continue;

      // Test whether %(to)s has coefficient named 'name'
      try {
        %(to)s->coefficient_number(name);
      } catch (...) {
        std::cout << "Attaching coefficient named: " << name << " to %(to)s";
        std::cout << " failed! But this might be expected." << std::endl;
        continue;
      }
      %(to)s->set_coefficient(name, %(from)s.coefficient(i));
    }
    """


attach_domains_template = """
    // Attach subdomains from %(from)s to %(to)s
    %(to)s->dx = %(from)s.cell_domains();
    %(to)s->ds = %(from)s.exterior_facet_domains();
    %(to)s->dS = %(from)s.interior_facet_domains();
"""


update_ec_template = """
  /// Initialize all error control forms, attach coefficients and
  /// (re-)set error control
  virtual void update_ec(const dolfin::Form& a, const dolfin::Form& L)
  {
    // This stuff is created here and shipped elsewhere
    std::shared_ptr<dolfin::Form> a_star;           // Dual lhs
    std::shared_ptr<dolfin::Form> L_star;           // Dual rhs
    std::shared_ptr<dolfin::FunctionSpace> V_Ez_h;  // Extrapolation space
    std::shared_ptr<dolfin::Function> Ez_h;         // Extrapolated dual
    std::shared_ptr<dolfin::Form> residual;         // Residual (as functional)
    std::shared_ptr<dolfin::FunctionSpace> V_R_T;   // Trial space for cell residual
    std::shared_ptr<dolfin::Form> a_R_T;            // Cell residual lhs
    std::shared_ptr<dolfin::Form> L_R_T;            // Cell residual rhs
    std::shared_ptr<dolfin::FunctionSpace> V_b_T;   // Function space for cell bubble
    std::shared_ptr<dolfin::Function> b_T;          // Cell bubble
    std::shared_ptr<dolfin::FunctionSpace> V_R_dT;  // Trial space for facet residual
    std::shared_ptr<dolfin::Form> a_R_dT;           // Facet residual lhs
    std::shared_ptr<dolfin::Form> L_R_dT;           // Facet residual rhs
    std::shared_ptr<dolfin::FunctionSpace> V_b_e;   // Function space for cell cone
    std::shared_ptr<dolfin::Function> b_e;          // Cell cone
    std::shared_ptr<dolfin::FunctionSpace> V_eta_T; // Function space for indicators
    std::shared_ptr<dolfin::Form> eta_T;            // Indicator form

    // Some handy views
    assert(a.function_space(0));
    const auto Vhat = a.function_space(0); // Primal test

    assert(a.function_space(0));
    const auto V = a.function_space(1);    // Primal trial

    assert(V->mesh());
    const auto mesh = V->mesh();

    std::string name;

    // Initialize dual forms
    a_star.reset(new %(a_star)s(V, Vhat));
    L_star.reset(new %(L_star)s(V));

    %(attach_a_star)s
    %(attach_L_star)s

    // Initialize residual
    residual.reset(new %(residual)s(mesh));
    %(attach_residual)s

    // Initialize extrapolation space and (fake) extrapolation
    V_Ez_h.reset(new %(V_Ez_h)s(mesh));
    Ez_h.reset(new dolfin::Function(V_Ez_h));
    residual->set_coefficient("__improved_dual", Ez_h);

    // Create bilinear and linear form for computing cell residual R_T
    V_R_T.reset(new %(V_R_T)s(mesh));
    a_R_T.reset(new %(a_R_T)s(V_R_T, V_R_T));
    L_R_T.reset(new %(L_R_T)s(V_R_T));

    // Initialize bubble and attach to a_R_T and L_R_T
    V_b_T.reset(new %(V_b_T)s(mesh));
    b_T.reset(new dolfin::Function(V_b_T));
    *b_T->vector() = 1.0;
    %(attach_L_R_T)s

    // Attach bubble function to _a_R_T and _L_R_T
    a_R_T->set_coefficient("__cell_bubble", b_T);
    L_R_T->set_coefficient("__cell_bubble", b_T);

    // Create bilinear and linear form for computing facet residual R_dT
    V_R_dT.reset(new %(V_R_dT)s(mesh));
    a_R_dT.reset(new %(a_R_dT)s(V_R_dT, V_R_dT));
    L_R_dT.reset(new %(L_R_dT)s(V_R_dT));
    %(attach_L_R_dT)s

    // Initialize (fake) cone and attach to a_R_dT and L_R_dT
    V_b_e.reset(new %(V_b_e)s(mesh));
    b_e.reset(new dolfin::Function(V_b_e));
    a_R_dT->set_coefficient("__cell_cone", b_e);
    L_R_dT->set_coefficient("__cell_cone", b_e);

    // Create error indicator form
    V_eta_T.reset(new %(V_eta_T)s(mesh));
    eta_T.reset(new %(eta_T)s(V_eta_T));

    // Update error control
    _ec.reset(new dolfin::ErrorControl(a_star, L_star, residual,
                                       a_R_T, L_R_T, a_R_dT, L_R_dT, eta_T,
                                       %(linear)s));
  }
"""


def _attach(tos, froms):

    if not isinstance(froms, tuple):
        return attach_coefficient_template % {"to": tos, "from": froms} \
            + attach_domains_template % {"to": tos, "from": froms}

    # NB: If multiple forms involved, attach domains from last form.
    coeffs = "\n".join([attach_coefficient_template % {"to": to, "from": fro}
                        for (to, fro) in zip(tos, froms)])
    domains = attach_domains_template % {"to": tos[-1], "from": froms[-1]}
    return coeffs + domains


def generate_maps(linear):
    """
    NB: This depends on the ordering of the forms
    """
    maps = {"a_star": "Form_%d" % 0,
            "L_star": "Form_%d" % 1,
            "residual": "Form_%d" % 2,
            "a_R_T": "Form_%d" % 3,
            "L_R_T": "Form_%d" % 4,
            "a_R_dT": "Form_%d" % 5,
            "L_R_dT": "Form_%d" % 6,
            "eta_T": "Form_%d" % 7,
            "V_Ez_h": "CoefficientSpace_%s" % "__improved_dual",
            "V_R_T": "Form_%d::TestSpace" % 4,
            "V_b_T": "CoefficientSpace_%s" % "__cell_bubble",
            "V_R_dT": "Form_%d::TestSpace" % 6,
            "V_b_e": "CoefficientSpace_%s" % "__cell_cone",
            "V_eta_T": "Form_%d::TestSpace" % 7,
            "attach_a_star": _attach("a_star", "a"),
            "attach_L_star": _attach("L_star", "(*this)"),
            "attach_residual": _attach(("residual",) * 2, ("a", "L")),
            "attach_L_R_T": _attach(("L_R_T",) * 2, ("a", "L")),
            "attach_L_R_dT": _attach(("L_R_dT",) * 2, ("a", "L")),
            "linear": "true" if linear else "false"
            }
    return maps


def generate_update_ec(form):

    linear = "__discrete_primal_solution" in form.coefficient_names
    maps = generate_maps(linear)
    code = update_ec_template % maps
    return code