[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Wed, 4 Jul 2018 07:40:56 -0400 (EDT) |
branch: devel-yves-double-integ
commit 36a8a027f9d45dea53eee7e0f1dc709f9a391602
Author: Yves Renard <address@hidden>
Date: Wed Jul 4 13:40:37 2018 +0200
A first version of working assembly on the direct product of two domains
---
doc/sphinx/source/userdoc/gasm_high.rst | 86 ++-
interface/src/gf_asm.cc | 89 +--
interface/src/gf_model_set.cc | 88 ++-
interface/tests/python/Makefile.am | 2 +
interface/tests/python/check_secondary_domain.py | 91 +++
src/getfem/getfem_generic_assembly.h | 37 +-
.../getfem_generic_assembly_compile_and_exec.h | 22 +-
src/getfem/getfem_generic_assembly_tree.h | 19 +-
src/getfem/getfem_models.h | 60 +-
src/getfem_generic_assembly_compile_and_exec.cc | 692 +++++++++++++++++----
src/getfem_generic_assembly_interpolation.cc | 35 ++
src/getfem_generic_assembly_semantic.cc | 30 +-
src/getfem_generic_assembly_tree.cc | 15 +-
src/getfem_generic_assembly_workspace.cc | 11 +-
src/getfem_models.cc | 124 +++-
15 files changed, 1142 insertions(+), 259 deletions(-)
diff --git a/doc/sphinx/source/userdoc/gasm_high.rst
b/doc/sphinx/source/userdoc/gasm_high.rst
index e3ef022..94c791e 100644
--- a/doc/sphinx/source/userdoc/gasm_high.rst
+++ b/doc/sphinx/source/userdoc/gasm_high.rst
@@ -8,8 +8,8 @@
.. _ud-gasm-high:
-Compute arbitrary terms - high-level generic assembly procedures
-================================================================
+Compute arbitrary terms - high-level generic assembly procedures - Weak-form
language
+=====================================================================================
This section presents what is now the main generic assembly of |gf|. It is a
high-level generic assembly in the sense that it is based on a weak form
language to describe the weak formulation of boundary value problems of partial
differential equations. It mainly has been developed to circumvent the
difficulties with the previous low-level generic assembly (see
:ref:`ud-gasm-low`) for which nonlinear terms were quite difficult to describe.
Conversely, a symbolic differentiation algorith [...]
@@ -72,6 +72,7 @@ A specific weak form language has been developed to describe
the weak formulatio
- ``Elementary_transformation(variable, transformation)``: Allow a linear
tranformation defined at the element level (i.e. not possible to define at the
gauss point level). This feature has been added mostly for defining a reduction
for plate elements (projection onto low-level vector element such as rotated
RT0). ``transformation`` is an object stored by the workspace or model object
which describes the trasformation for a particular element.
+ - Possibility of integration on the direct product of two-domains for double
integral computation or coupling of two variables with a Kernel / convolution /
exchange integral. This allows terms like
:math:`\displaystyle\int_{\Omega_1}\int_{\Omega_2}k(x,y)u(x)v(y)dydx` with
:math:`\Omega_1` and :math:`\Omega_2` two domains, different or not, having
their own meshes, integration methods and with :math:`u` a variable defined on
:math:`\Omega_1` and :math:`v` a variable defined on :math:`\ [...]
Some basic examples
-------------------
@@ -161,9 +162,9 @@ The generic assembly is driven by the object
``getfem::ga_workspace`` defined in
There is two ways to define a ``getfem::ga_workspace`` object. It can depend
on a model (see :ref:`ud-model`) and should be declared as::
- getfem::ga_workspace workspace(md);
+ getfem::ga_workspace workspace(model);
-with ``md`` a previously define ``getfem::model`` object. In that case the
variable and constant considered are the one of the model. The second way it to
define an independent ``getfem::ga_workspace`` object by::
+with ``model`` a previously define ``getfem::model`` object. In that case the
variable and constant considered are the one of the model. The second way it to
define an independent ``getfem::ga_workspace`` object by::
getfem::ga_workspace workspace;
@@ -356,6 +357,10 @@ The last example is the assembly of the stiffness matrix
of an order four proble
with ``D`` the flexion modulus and ``nu`` the Poisson ratio.
+Script languages call of the assembly
+-------------------------------------
+
+For the use with Python, Scilab or Matlab interfaces, see the respective
documentation, in particular the ``gf_asm`` command and the ``model`` object.
The tensors
@@ -529,7 +534,7 @@ Constant expressions
- Floating points with standards notations (for instance ``3``, ``1.456``,
``1E-6``)
- ``pi``: the constant Pi.
- ``meshdim``: the dimension of the current mesh (i.e. size of geometrical
nodes)
- - ``timestep``: the main time step of the model on which this assembly
string is evaluated (defined by ``md.set_time_step(dt)``). Do not work on pure
workspaces.
+ - ``timestep``: the main time step of the model on which this assembly
string is evaluated (defined by ``model.set_time_step(dt)``). Do not work on
pure workspaces.
- ``Id(n)``: the identity matrix of size :math:`n\times n`. `n` should be an
integer expression. For instance ``Id(meshdim)`` is allowed.
- ``qdim(u)``: the total dimension of the variable ``u`` (i.e. the size for
fixed size variables and the total dimension of the vector/tensor field for FEM
variables)
- ``qdims(u)``: the dimensions of the variable ``u`` (i.e. the size for
fixed size variables and the vector of dimensions of the vector/tensor field
for FEM variables)
@@ -721,24 +726,24 @@ The transformation defined by an expression can be added
to the workspace or the
or::
add_interpolate_transformation_from_expression
- (md, transname, source_mesh, target_mesh, expr);
+ (model, transname, source_mesh, target_mesh, expr);
-where ``workspace`` is a workspace object, ``md`` a model object,
``transname`` is the name given to the transformation, ``source_mesh`` the mesh
on which the integration occurs, ``target_mesh`` the mesh on which the
interpolation is performed and ``expr`` is a regular expression of the
high-level generic weak form language which may contains reference to the
variables of the workspace/model.
+where ``workspace`` is a workspace object, ``model`` a model object,
``transname`` is the name given to the transformation, ``source_mesh`` the mesh
on which the integration occurs, ``target_mesh`` the mesh on which the
interpolation is performed and ``expr`` is a regular expression of the
high-level generic weak form language which may contains reference to the
variables of the workspace/model.
For instance, an expression::
add_interpolate_transformation_from_expression
- (md, "my_transformation", my_mesh, my_mesh, "X-[1;0]");
+ (model, "my_transformation", my_mesh, my_mesh, "X-[1;0]");
will allow to integrate some expressions at the current position with a shift
of -1 with respect to the first coordinate. This simple kind of transformation
can be used to prescribe a periodic condition.
Of course, one may used more complex expressions such as::
add_interpolate_transformation_from_expression
- (md, "my_transformation", my_mesh, my_second_mesh, "[X[1]cos(X[2]);
X[1]sin(X[2])]");
+ (model, "my_transformation", my_mesh, my_second_mesh, "[X[1]cos(X[2]);
X[1]sin(X[2])]");
add_interpolate_transformation_from_expression
- (md, "my_transformation", my_mesh, my_mesh, "X+u");
+ (model, "my_transformation", my_mesh, my_mesh, "X+u");
where ``u`` is a vector variable of the workspace/model.
@@ -783,7 +788,7 @@ Element extrapolation transformation
A specific transformation (see previous section) is defined in order to allows
the evaluation of certain quantities by extrapolation with respect to another
element (in general a neighbour element). This is not strictly speaking a
transformation since the point location remain unchanged, but the evaluation is
made on another element extrapolating the shape functions outside it. This
transformation is used for stabilization term in fictitious domain applications
(with cut elements) where [...]
add_element_extrapolation_transformation
- (md, transname, my_mesh, std::map<size_type, size_type> &elt_corr);
+ (model, transname, my_mesh, std::map<size_type, size_type> &elt_corr);
add_element_extrapolation_transformation
(workspace, transname, my_mesh, std::map<size_type, size_type> &elt_corr);
@@ -793,7 +798,7 @@ The map elt_corr should contain the correspondances between
the elements where t
The following functions allow to change the element correspondance of a
previously added element extrapolation transformation::
set_element_extrapolation_correspondance
- (md, transname, std::map<size_type, size_type> &elt_corr);
+ (model, transname, std::map<size_type, size_type> &elt_corr);
set_element_extrapolation_correspondance
(workspace, transname, std::map<size_type, size_type> &elt_corr);
@@ -837,6 +842,63 @@ See for instance
:file:`interface/tests/python/demo_laplacian_DG.py` or :file:`i
Compared to other interpolate transformations, this transformation is more
optimized and benefits from finite element and geometric transformation
pre-computations.
+.. _ud-gasm-high-secondary-dom:
+
+Double domain integrals or terms (convolution - Kernel - Exchange integrals)
+----------------------------------------------------------------------------
+
+In some very special cases, it can be interesting to compute an integral on
the direct product of two domains, i.e. a double integral such as for instance
+
+.. math::
+
+ \int_{\Omega_1}\int_{\Omega_2}k(x,y)u(x)v(y)dydx,
+
+where :math:`k(x,y)` is a given kernel, :math:`u` a quantity defined on
:math:`\Omega_1` and :math:`v` a quantity defined on :math:`\Omega_2`,
eventually with :math:`\Omega_1` and :math:`\Omega_2` the same domain. This
can be interesting either to compute such an integral or to define an
interaction term between two variables defined on two different domains.
+
+CAUTION: Of course, this kind of term have to be used with great care, since
it naturally leads to fully populated stiffness or tangent matrices.
+
+
+The weak form language of |gf| furnishes a mecanism to compute such a term.
First, the secondary domain has to be declared in the workspace/model with its
integration methods. The addition of a standard secondary domain can be done
with one of the two following functions::
+
+ add_standard_secondary_domain(model, domain_name, mim, region);
+
+ add_standard_secondary_domain(workspace, domain_name, mim, region);
+
+where ``model`` or ``workspace`` is the model or workspace where the secondary
domain has to be declared, ``domain_name`` is a string for the identification
of this domain together with the mesh region and integration method, ``mim``
the integration method and ``region`` a mesh region. Note that with these
standard secondary domains, the integration is done on the whole region for
each element of the primary domain. It can be interesting to implement specific
secondary domains restrictin [...]
+
+Once a secondary domain has been declared, it can be specified that a weak
form language expression has to be assembled on the direct product of a current
domain and a secondary domain, adding the name of the secondary domain to the
``add_expression`` method of the workspace object or using
``add_linear_twodomain_term``, ``add_nonlinear_twodomain_term`` or
``add_twodomain_source_term`` functions::
+ workspace.add_expression(expr, mim, region, derivative_order,
secondary_domain)
+ add_twodomain_source_term(model, mim, expr, region, secondary_domain)
+ add_linear_twodomain_term(model, mim, expr, region, secondary_domain)
+ add_nonlinear_twodomain_term(model, mim, expr, region, secondary_domain)
+
+For the utilisation with the Python/Scilab/Matlab interface, see the
documentation on ``gf_asm`` command and the ``model`` object.
+
+
+Inside an expression of the weak form language, one can refer to the unit
normal vector to a boundary, to the current position or to the value of a
variable thanks to the expressions::
+
+ Secondary_domain(Normal)
+ Secondary_domain(X)
+ Secondary_domain(u)
+ Secondary_domain(Grad_u)
+ Secondary_domain(Div_u)
+ Secondary_domain(Hess_u)
+ Secondary_domain(Test_u)
+ Secondary_domain(Grad_Test_u)
+ Secondary_domain(Div_Test_u)
+ Secondary_domain(Hess_Test_u)
+
+For instance, a term like
+
+.. math::
+
+ \int_{\Omega_1}\int_{\Omega_1}e^{-\|x-y\|}u(x)u(y)dydx,
+
+would correspond to the following weak for language expression::
+
+ exp(Norm(X-Secondary_domain(X)))*u*Secondary_domain(u)
+
+
.. _ud-gasm-high-elem-trans:
Elementary transformations
diff --git a/interface/src/gf_asm.cc b/interface/src/gf_asm.cc
index eb77267..5792ecf 100644
--- a/interface/src/gf_asm.cc
+++ b/interface/src/gf_asm.cc
@@ -427,6 +427,9 @@ static void do_high_level_generic_assembly(mexargs_in& in,
mexargs_out& out) {
getfem::model dummy_md;
bool with_model = in.remaining() && is_model_object(in.front());
const getfem::model &md = with_model ? *to_model_object(in.pop()) : dummy_md;
+ bool with_secondary = false;
+ std::string secondary_domain;
+
getfem::ga_workspace workspace2(md);
getfem::ga_workspace &workspace = with_model ? workspace2 : workspace1;
@@ -436,46 +439,54 @@ static void do_high_level_generic_assembly(mexargs_in&
in, mexargs_out& out) {
while (in.remaining()) {
std::string varname = in.pop().to_string();
- bool is_cte = (in.pop().to_integer() == 0);
- const getfem::mesh_fem *mf(0);
- const getfem::im_data *mimd(0);
- if (is_meshfem_object(in.front()))
- mf = to_meshfem_object(in.pop());
- else if (is_meshimdata_object(in.front()))
- mimd = to_meshimdata_object(in.pop());
- darray U = in.pop().to_darray();
- GMM_ASSERT1(vectors.find(varname) == vectors.end(),
- "The same variable/constant name is repeated twice: "
- << varname)
- GMM_ASSERT1(!with_model || !md.variable_exists(varname),
- "The same variable/constant name is already defined in "
- "the model: " << varname)
- gmm::resize(vectors[varname], U.size());
- gmm::copy(U, vectors[varname]);
- if (is_cte) {
- if (mf)
- workspace.add_fem_constant(varname, *mf, vectors[varname]);
- else if (mimd)
- workspace.add_im_data(varname, *mimd, vectors[varname]);
- else
- workspace.add_fixed_size_constant(varname, vectors[varname]);
+ if (varname.compare("Secondary_domain") == 0 ||
+ varname.compare("Secondary_Domain") == 0) {
+ GMM_ASSERT1(!with_secondary,
+ "Only one secondary domain can be specified");
+ secondary_domain = in.pop().to_string();
+ with_secondary = true;
} else {
- if (mf) {
- gmm::sub_interval I(nbdof, mf->nb_dof());
- nbdof += mf->nb_dof();
- workspace.add_fem_variable(varname, *mf, I, vectors[varname]);
- } else if (mimd) {
- THROW_BADARG("Data defined on integration points can not be a
variable");
- } else {
- gmm::sub_interval I(nbdof, U.size());
- nbdof += U.size();
- workspace.add_fixed_size_variable(varname, I, vectors[varname]);
+ bool is_cte = (in.pop().to_integer() == 0);
+ const getfem::mesh_fem *mf(0);
+ const getfem::im_data *mimd(0);
+ if (is_meshfem_object(in.front()))
+ mf = to_meshfem_object(in.pop());
+ else if (is_meshimdata_object(in.front()))
+ mimd = to_meshimdata_object(in.pop());
+ darray U = in.pop().to_darray();
+ GMM_ASSERT1(vectors.find(varname) == vectors.end(),
+ "The same variable/constant name is repeated twice: "
+ << varname);
+ GMM_ASSERT1(!with_model || !md.variable_exists(varname),
+ "The same variable/constant name is already defined in "
+ "the model: " << varname);
+ gmm::resize(vectors[varname], U.size());
+ gmm::copy(U, vectors[varname]);
+ if (is_cte) {
+ if (mf)
+ workspace.add_fem_constant(varname, *mf, vectors[varname]);
+ else if (mimd)
+ workspace.add_im_data(varname, *mimd, vectors[varname]);
+ else
+ workspace.add_fixed_size_constant(varname, vectors[varname]);
+ } else {
+ if (mf) {
+ gmm::sub_interval I(nbdof, mf->nb_dof());
+ nbdof += mf->nb_dof();
+ workspace.add_fem_variable(varname, *mf, I, vectors[varname]);
+ } else if (mimd) {
+ THROW_BADARG("Data defined on integration points can not be a
variable");
+ } else {
+ gmm::sub_interval I(nbdof, U.size());
+ nbdof += U.size();
+ workspace.add_fixed_size_variable(varname, I, vectors[varname]);
+ }
}
}
}
-
- size_type add_derivative_order = (order != 0) ? 2 : 0;
- workspace.add_expression(expr, *mim, region, add_derivative_order);
+
+ size_type der_order = (order != 0) ? 2 : 0;
+ workspace.add_expression(expr, *mim, region, der_order, secondary_domain);
switch (order) {
case 0:
@@ -727,7 +738,7 @@ void gf_asm(getfemint::mexargs_in& m_in,
getfemint::mexargs_out& m_out) {
if (subc_tab.size() == 0) {
- /address@hidden @CELL{...} = ('generic', @tmim mim, @int order, @str
expression, @int region, address@hidden model,] address@hidden varname, @int
is_variable[, address@hidden mf, @tmimd mimd}], value], ...)
+ /address@hidden @CELL{...} = ('generic', @tmim mim, @int order, @str
expression, @int region, address@hidden model, ['Secondary_domain', 'name',]]
address@hidden varname, @int is_variable[, address@hidden mf, @tmimd mimd}],
value], ...)
High-level generic assembly procedure for volumic or boundary assembly.
Performs the generic assembly of `expression` with the integration
@@ -741,7 +752,9 @@ void gf_asm(getfemint::mexargs_in& m_in,
getfemint::mexargs_out& m_out) {
tangent (matrix) (order = 2) is to be computed.
`model` is an optional parameter allowing to take into account
- all variables and data of a model.
+ all variables and data of a model. Optionnally, for the integration
+ on the product of two domains, a secondary domain of the model can
+ be specified after a 'Secondary_domain' string.
The variables and constant (data) are listed after the
region number (or optionally the model).
diff --git a/interface/src/gf_model_set.cc b/interface/src/gf_model_set.cc
index 3680dc4..67e05e4 100644
--- a/interface/src/gf_model_set.cc
+++ b/interface/src/gf_model_set.cc
@@ -474,7 +474,21 @@ void gf_model_set(getfemint::mexargs_in& m_in,
elt_corr);
);
- /address@hidden ('set element extrapolation correspondance', @str
transname, @mat elt_corr)
+ /address@hidden ('add standard secondary domain', @str name, @tmim mim,
@int region = -1)
+ Add a secondary domain to the model which can be used in a weak-form
language expression for integration on the product of two domains. `name` is
the name
+ of the secondary domain, `mim` is an integration method on this domain
+ and `region` the region on which the integration is to be performed. @*/
+ sub_command
+ ("add standard secondary domain", 3, 3, 0, 0,
+ std::string name = in.pop().to_string();
+ getfem::mesh_im *mim = to_meshim_object(in.pop());
+ size_type region = size_type(-1);
+ if (in.remaining()) region = in.pop().to_integer();
+ add_standard_secondary_domain(*md, name, *mim, region);
+ );
+
+
+ /address@hidden ('set element extrapolation correspondance', @str
transname, @mat elt_corr)
Change the correspondance map of an element extrapolation interpolate
transformation. @*/
sub_command
@@ -632,6 +646,30 @@ void gf_model_set(getfemint::mexargs_in& m_in,
out.pop().from_integer(int(ind));
);
+ /address@hidden ind = ('add linear twodomain term', @tmim mim, @str
expression, @int region, @str secondary_domain[, @int is_symmetric[, @int
is_coercive]])
+ Adds a linear term given by a weak form language expression like
+ MODEL:SET('add linear term') but for an integration on a direct
+ product of two domains, a first specfied by ``mim`` and ``region``
+ and a second one by ``secondary_domain`` which has to be declared
+ first into the address@hidden/
+ sub_command
+ ("add linear twodomain term", 4, 6, 0, 1,
+ getfem::mesh_im *mim = to_meshim_object(in.pop());
+ std::string expr = in.pop().to_string();
+ size_type region = in.pop().to_integer();
+ std::string secdom = in.pop().to_string();
+ int is_symmetric = 0;
+ if (in.remaining()) is_symmetric = in.pop().to_integer();
+ int is_coercive = 0;
+ if (in.remaining()) is_coercive = in.pop().to_integer();
+
+ size_type ind = getfem::add_linear_twodomain_term
+ (*md, *mim, expr, region, secdom, is_symmetric,
+ is_coercive) + config::base_index();
+ workspace().set_dependence(md, mim);
+ out.pop().from_integer(int(ind));
+ );
+
/address@hidden ind = ('add linear generic assembly brick', @tmim mim,
@str expression[, @int region[, @int is_symmetric[, @int is_coercive]]])
Deprecated. Use MODEL:SET('add linear term') instead. @*/
sub_command
@@ -683,6 +721,32 @@ void gf_model_set(getfemint::mexargs_in& m_in,
out.pop().from_integer(int(ind));
);
+ /address@hidden ind = ('add nonlinear twodomain term', @tmim mim, @str
expression, @int region, @str secondary_domain[, @int is_symmetric[, @int
is_coercive]])
+ Adds a nonlinear term given by a weak form language expression like
+ MODEL:SET('add nonlinear term') but for an integration on a direct
+ product of two domains, a first specfied by ``mim`` and ``region``
+ and a second one by ``secondary_domain`` which has to be declared
+ first into the address@hidden/
+ sub_command
+ ("add nonlinear twodomain term", 4, 6, 0, 1,
+ getfem::mesh_im *mim = to_meshim_object(in.pop());
+ std::string expr = in.pop().to_string();
+ size_type region = in.pop().to_integer();
+ std::string secdom = in.pop().to_string();
+ int is_symmetric = 0;
+ if (in.remaining()) is_symmetric = in.pop().to_integer();
+ int is_coercive = 0;
+ if (in.remaining()) is_coercive = in.pop().to_integer();
+
+ size_type ind
+ = getfem::add_nonlinear_twodomain_term
+ (*md, *mim, expr, region, secdom, is_symmetric,
+ is_coercive) + config::base_index();
+ workspace().set_dependence(md, mim);
+ out.pop().from_integer(int(ind));
+ );
+
+
/address@hidden ind = ('add nonlinear generic assembly brick', @tmim mim,
@str expression[, @int region[, @int is_symmetric[, @int is_coercive]]])
Deprecated. Use MODEL:SET('add nonlinear term') address@hidden/
sub_command
@@ -720,12 +784,32 @@ void gf_model_set(getfemint::mexargs_in& m_in,
if (in.remaining()) region = in.pop().to_integer();
size_type ind
- = getfem::add_source_term_generic_assembly_brick
+ = getfem::add_source_term
(*md, *mim, expr, region) + config::base_index();
workspace().set_dependence(md, mim);
out.pop().from_integer(int(ind));
);
+ /address@hidden ind = ('add twodomain source term', @tmim mim, @str
expression, @int region, @str secondary_domain)
+ Adds a source term given by a weak form language expression like
+ MODEL:SET('add source term') but for an integration on a direct
+ product of two domains, a first specfied by ``mim`` and ``region``
+ and a second one by ``secondary_domain`` which has to be declared
+ first into the address@hidden/
+ sub_command
+ ("add twodomain source term", 4, 4, 0, 1,
+ getfem::mesh_im *mim = to_meshim_object(in.pop());
+ std::string expr = in.pop().to_string();
+ size_type region = in.pop().to_integer();
+ std::string secdom = in.pop().to_string();
+
+ size_type ind
+ = getfem::add_twodomain_source_term
+ (*md, *mim, expr, region, secdom) + config::base_index();
+ workspace().set_dependence(md, mim);
+ out.pop().from_integer(int(ind));
+ );
+
/address@hidden ind = ('add source term generic assembly brick', @tmim
mim, @str expression[, @int region])
Deprecated. Use MODEL:SET('add source term') instead. @*/
sub_command
diff --git a/interface/tests/python/Makefile.am
b/interface/tests/python/Makefile.am
index bd3c0ba..9c2b1f0 100644
--- a/interface/tests/python/Makefile.am
+++ b/interface/tests/python/Makefile.am
@@ -27,6 +27,7 @@ EXTRA_DIST= \
check_global_functions.py \
check_levelset.py \
check_asm.py \
+ check_secondary_domain.py \
check_mixed_mesh.py \
demo_crack.py \
demo_fictitious_domains.py \
@@ -61,6 +62,7 @@ TESTS = \
check_export.py \
check_global_functions.py \
check_asm.py \
+ check_secondary_domain.py \
check_mixed_mesh.py \
demo_wave.py \
demo_laplacian.py \
diff --git a/interface/tests/python/check_secondary_domain.py
b/interface/tests/python/check_secondary_domain.py
new file mode 100644
index 0000000..28ee9dd
--- /dev/null
+++ b/interface/tests/python/check_secondary_domain.py
@@ -0,0 +1,91 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Python GetFEM++ interface
+#
+# Copyright (C) 2018-2018 Yves Renard.
+#
+# This file is a part of GetFEM++
+#
+# GetFEM++ 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 2.1 of the License, or
+# (at your option) any later version.
+# This program 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 this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+#
+############################################################################
+""" Test of the assembly on the direct product of two domains.
+
+ This program is used to check that Python-GetFEM interface, and more
+ generally GetFEM are working. It focuses on testing some operations
+ of the high generic assembly language.
+
+ $Id$
+"""
+import numpy as np
+import getfem as gf
+import os
+
+
+NX = 4
+
+m1 = gf.Mesh('cartesian', np.arange(0,1+1./NX,1./NX)) # Structured 1D mesh
+mfu1 = gf.MeshFem(m1, 1); mfu1.set_fem(gf.Fem('FEM_PK(1,1)'))
+mim1 = gf.MeshIm(m1, gf.Integ('IM_GAUSS1D(4)'))
+U1 = mfu1.eval('x')
+
+
+m2 = gf.Mesh('triangles grid', np.arange(0,1+1./NX,1./NX),
+ np.arange(0,1+1./NX,1./NX)) # Structured 2D mesh
+mfu2 = gf.MeshFem(m2, 1); mfu2.set_fem(gf.Fem('FEM_PK(2,1)'))
+mim2 = gf.MeshIm(m2, gf.Integ('IM_TRIANGLE(4)'))
+U2 = mfu2.eval('x+y')
+
+
+
+md = gf.Model('real')
+
+md.add_fem_variable('u1', mfu1)
+md.set_variable('u1', U1)
+md.add_fem_variable('u2', mfu2)
+md.set_variable('u2', U2)
+
+md.add_standard_secondary_domain('secdom', mim2);
+
+
+result = gf.asm('generic', mim1, 0, "u1*Secondary_Domain(u2)", -1, md,
+ "Secondary_domain", "secdom")
+if (abs(result-0.5) > 1e-8) : print("Bad value"); exit(1)
+
+result = gf.asm('generic', mim1, 0, "u1*(Secondary_Domain(Grad_u2)(1))", -1,
md,
+ "Secondary_domain", "secdom")
+if (abs(result-0.5) > 1e-8) : print("Bad value"); exit(1)
+
+result = gf.asm('generic', mim1, 0,
+ "u1*(Secondary_Domain(X)(1)+Secondary_Domain(X)(2))", -1, md,
+ "Secondary_domain", "secdom")
+if (abs(result-0.5) > 1e-8) : print("Bad value"); exit(1)
+
+result = gf.asm('generic', mim1, 0, "u1*sqr(Secondary_Domain(u2))", -1, md,
+ "Secondary_domain", "secdom")
+if (abs(result-7./12.) > 1e-8) : print("Bad value"); exit(1)
+
+
+
+
+
+M = gf.asm('generic', mim1, 2, "Test_u1*Secondary_Domain(Test2_u2)", -1, md,
+ "Secondary_domain", "secdom")
+V1 = gf.asm('generic', mim1, 1, "Test_u1", -1, md)
+V2 = gf.asm('generic', mim2, 1, "Test_u2", -1, md)
+
+for i in range(0, V1.size):
+ for j in range(0, V2.size):
+ if (abs(M[i,j] - V1[i]*V2[j]) > 1E-8):
+ print("Bad value for matrix assembly"); exit(1)
+
diff --git a/src/getfem/getfem_generic_assembly.h
b/src/getfem/getfem_generic_assembly.h
index 1072a47..265ba11 100644
--- a/src/getfem/getfem_generic_assembly.h
+++ b/src/getfem/getfem_generic_assembly.h
@@ -133,23 +133,24 @@ namespace getfem {
pelementary_transformation;
//=========================================================================
- // Virtual secundary_domain object.
+ // Virtual secondary_domain object.
//=========================================================================
class APIDECL virtual_secondary_domain {
- const mesh_im &mim;
- size_type region;
+ protected:
+ const mesh_im &mim_;
+ const mesh_region region;
public:
- // interface to be reconsidered
- virtual void give_secondary_elements(const mesh &m, size_type cv,
- std::set<size_type> &elt_lst) const = 0;
- virtual void init(const ga_workspace &workspace) const = 0;
- virtual void finalize() const = 0;
+ const mesh_im &mim(void) const { return mim_; }
+ virtual const mesh_region &give_region(const mesh &m,
+ size_type cv, short_type f) const = 0;
+ // virtual void init(const ga_workspace &workspace) const = 0;
+ // virtual void finalize() const = 0;
- virtual_secondary_domain(const mesh_im &mim_, size_type region_)
- : mim(mim_), region(region_) {}
+ virtual_secondary_domain(const mesh_im &mim__, const mesh_region ®ion_)
+ : mim_(mim__), region(region_) {}
virtual ~virtual_secondary_domain() {}
};
@@ -310,6 +311,7 @@ namespace getfem {
std::string varname_interpolation; // Where to interpolate
std::string name_test1, name_test2;
std::string interpolate_name_test1, interpolate_name_test2;
+ std::string secondary_domain;
const mesh_im *mim;
const mesh *m;
const mesh_region *rg;
@@ -406,7 +408,8 @@ namespace getfem {
*/
size_type add_expression(const std::string &expr, const mesh_im &mim,
const mesh_region &rg=mesh_region::all_convexes(),
- size_type add_derivative_order = 2);
+ size_type add_derivative_order = 2,
+ const std::string &secondary_dom = "");
/* Internal use */
void add_function_expression(const std::string &expr);
/* Internal use */
@@ -773,7 +776,17 @@ namespace getfem {
(ga_workspace &workspace, const std::string &name,
std::map<size_type, size_type> &elt_corr);
-
+ //=========================================================================
+ // Secondary domains
+ //=========================================================================
+
+ void add_standard_secondary_domain
+ (model &md, const std::string &name, const mesh_im &mim,
+ const mesh_region &rg=mesh_region::all_convexes());
+
+ void add_standard_secondary_domain
+ (ga_workspace &workspace, const std::string &name, const mesh_im &mim,
+ const mesh_region &rg=mesh_region::all_convexes());
} /* end of namespace getfem. */
diff --git a/src/getfem/getfem_generic_assembly_compile_and_exec.h
b/src/getfem/getfem_generic_assembly_compile_and_exec.h
index a994cfd..a9eb562 100644
--- a/src/getfem/getfem_generic_assembly_compile_and_exec.h
+++ b/src/getfem/getfem_generic_assembly_compile_and_exec.h
@@ -100,12 +100,14 @@ namespace getfem {
std::map<gauss_pt_corresp, bgeot::pstored_point_tab> neighbour_corresp;
struct region_mim
- : std::tuple<const mesh_im *, const mesh_region *, char *> {
+ : std::tuple<const mesh_im *, const mesh_region *, psecondary_domain> {
const mesh_im* mim() const { return std::get<0>(*this); }
const mesh_region* region() const { return std::get<1>(*this); }
- region_mim(const mesh_im *mim_, const mesh_region *region_) :
- std::tuple<const mesh_im *, const mesh_region *,char *>(mim_, region_, 0)
- {}
+ psecondary_domain psd() const { return std::get<2>(*this); }
+ region_mim(const mesh_im *mim_, const mesh_region *region_,
+ psecondary_domain psd) :
+ std::tuple<const mesh_im *, const mesh_region *, psecondary_domain>
+ (mim_, region_, psd) {}
};
std::map<std::string, const base_vector *> extended_vars;
@@ -135,13 +137,21 @@ namespace getfem {
std::map<const mesh_fem *, pfem_precomp> pfps;
};
+
struct secondary_domain_info {
+ // const mesh *m;
papprox_integration pai;
- bool has_ctx;
fem_interpolation_context ctx;
base_small_vector Normal;
+
+ std::map<std::string, base_vector> local_dofs;
+ std::map<const mesh_fem *, pfem_precomp> pfps;
+ std::map<const mesh_fem *, base_tensor> base;
+ std::map<const mesh_fem *, base_tensor> grad;
+ std::map<const mesh_fem *, base_tensor> hess;
};
+
struct elementary_trans_info {
base_matrix M;
const mesh_fem *mf;
@@ -176,7 +186,7 @@ namespace getfem {
std::set<std::string> transformations_der;
std::map<std::string, interpolate_info> interpolate_infos;
std::map<std::string, elementary_trans_info> elementary_trans_infos;
- std::map<std::string, secondary_domain_info> secondary_domain_infos;
+ secondary_domain_info secondary_domain_infos;
// Instructions being executed at the first Gauss point after
// a change of integration method only.
diff --git a/src/getfem/getfem_generic_assembly_tree.h
b/src/getfem/getfem_generic_assembly_tree.h
index 2095d21..9f54f21 100644
--- a/src/getfem/getfem_generic_assembly_tree.h
+++ b/src/getfem/getfem_generic_assembly_tree.h
@@ -416,6 +416,7 @@ namespace getfem {
struct ga_tree {
pga_tree_node root, current_node;
+ std::string secondary_domain;
void add_scalar(scalar_type val, size_type pos, pstring expr);
void add_allindices(size_type pos, pstring expr);
@@ -439,16 +440,22 @@ namespace getfem {
{ duplicate_with_operation(pnode, GA_MINUS); }
void insert_node(pga_tree_node pnode, GA_NODE_TYPE node_type);
void add_child(pga_tree_node pnode, GA_NODE_TYPE node_type = GA_NODE_VOID);
- void swap(ga_tree &tree)
- { std::swap(root, tree.root); std::swap(current_node, tree.current_node); }
+ void swap(ga_tree &tree) {
+ std::swap(root, tree.root);
+ std::swap(current_node, tree.current_node);
+ std::swap(secondary_domain, tree.secondary_domain);
+ }
- ga_tree() : root(nullptr), current_node(nullptr) {}
+ ga_tree() : root(nullptr), current_node(nullptr), secondary_domain() {}
- ga_tree(const ga_tree &tree) : root(nullptr), current_node(nullptr)
+ ga_tree(const ga_tree &tree) : root(nullptr), current_node(nullptr),
+ secondary_domain(tree.secondary_domain)
{ if (tree.root) copy_node(tree.root, nullptr, root); }
- ga_tree &operator = (const ga_tree &tree)
- { clear(); if (tree.root) copy_node(tree.root,nullptr,root); return *this;
}
+ ga_tree &operator =(const ga_tree &tree) {
+ clear(); secondary_domain = tree.secondary_domain;
+ if (tree.root) copy_node(tree.root,nullptr,root); return *this;
+ }
~ga_tree() { clear(); }
};
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index 8704f36..11b397e 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -351,8 +351,10 @@ namespace getfem {
std::string expr;
const mesh_im &mim;
size_type region;
+ std::string secondary_domain;
gen_expr(const std::string &expr_, const mesh_im &mim_,
- size_type region_) : expr(expr_), mim(mim_), region(region_) {}
+ size_type region_, const std::string &secdom)
+ : expr(expr_), mim(mim_), region(region_), secondary_domain(secdom) {}
};
// Structure for assignment in assembly
@@ -397,8 +399,11 @@ namespace getfem {
public:
void add_generic_expression(const std::string &expr, const mesh_im &mim,
- size_type region) const
- { generic_expressions.push_back(gen_expr(expr, mim, region)); }
+ size_type region,
+ const std::string &secondary_domain = "") const
{
+ generic_expressions.push_back(gen_expr(expr, mim, region,
+ secondary_domain));
+ }
void add_external_load(size_type ib, scalar_type e) const
{ bricks[ib].external_load = e; }
scalar_type approx_external_load() { return approx_external_load_; }
@@ -1536,7 +1541,7 @@ namespace getfem {
//
//=========================================================================
- /** Add a matrix term given by the assembly string `expr` which will
+ /** Add a term given by the weak form language expression `expr` which will
be assembled in region `region` and with the integration method `mim`.
Only the matrix term will be taken into account, assuming that it is
linear.
@@ -1568,8 +1573,9 @@ namespace getfem {
is_coercive, brickname, return_if_nonlin);
}
- /** Add a nonlinear term given by the assembly string `expr` which will
- be assembled in region `region` and with the integration method `mim`.
+ /** Add a nonlinear term given by the weak form language expression `expr`
+ which will be assembled in region `region` and with the integration
+ method `mim`.
The expression can describe a potential or a weak form. Second order
terms (i.e. containing second order test functions, Test2) are not
allowed.
@@ -1577,12 +1583,12 @@ namespace getfem {
If you are not sure, the better is to declare the term not symmetric
and not coercive. But some solvers (conjugate gradient for instance)
are not allowed for non-coercive problems.
- `brickname` is an otpional name for the brick.
+ `brickname` is an optional name for the brick.
*/
size_type APIDECL add_nonlinear_term
(model &md, const mesh_im &mim, const std::string &expr,
size_type region = size_type(-1), bool is_sym = false,
- bool is_coercive = false, std::string brickname = "");
+ bool is_coercive = false, const std::string &brickname = "");
inline size_type APIDECL add_nonlinear_generic_assembly_brick
(model &md, const mesh_im &mim, const std::string &expr,
@@ -1617,6 +1623,44 @@ namespace getfem {
directvarname, directdataname, return_if_nonlin);
}
+ /** Adds a linear term given by a weak form language expression like
+ ``add_linear_term`` function but for an integration on a direct
+ product of two domains, a first specfied by ``mim`` and ``region``
+ and a second one by ``secondary_domain`` which has to be declared
+ first into the model.
+ */
+ size_type APIDECL add_linear_twodomain_term
+ (model &md, const mesh_im &mim, const std::string &expr,
+ size_type region, const std::string &secondary_domain,
+ bool is_sym = false, bool is_coercive = false, std::string brickname = "",
+ bool return_if_nonlin = false);
+
+ /** Adds a nonlinear term given by a weak form language expression like
+ ``add_nonlinear_term`` function but for an integration on a direct
+ product of two domains, a first specfied by ``mim`` and ``region``
+ and a second one by ``secondary_domain`` which has to be declared
+ first into the model.
+ */
+ size_type APIDECL add_nonlinear_twodomain_term
+ (model &md, const mesh_im &mim, const std::string &expr,
+ size_type region, const std::string &secondary_domain,
+ bool is_sym = false, bool is_coercive = false,
+ const std::string &brickname = "");
+
+ /** Adds a source term given by a weak form language expression like
+ ``add_source_term`` function but for an integration on a direct
+ product of two domains, a first specfied by ``mim`` and ``region``
+ and a second one by ``secondary_domain`` which has to be declared
+ first into the model.
+ */
+ size_type APIDECL add_twodomain_source_term
+ (model &md, const mesh_im &mim, const std::string &expr,
+ size_type region, const std::string &secondary_domain,
+ std::string brickname = "", std::string directvarname = std::string(),
+ const std::string &directdataname = std::string(),
+ bool return_if_nonlin = false);
+
+
/** Add a Laplacian term on the variable `varname` (in fact with a minus :
:math:`-\text{div}(\nabla u)`). If it is a vector
valued variable, the Laplacian term is componentwise. `region` is an
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index 1828176..d408997 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -4765,13 +4765,17 @@ namespace getfem {
pctx1 = &(gis.ctx);
const std::string &intn1 = pnode->interpolate_name_test1;
if (intn1.size()) {
- tensor_to_adapt = true;
- pctx1 = &(rmi.interpolate_infos[intn1].ctx);
- if (workspace.variable_group_exists(pnode->name_test1)) {
- ga_instruction_set::variable_group_info &vgi =
- rmi.interpolate_infos[intn1].groups_info[pnode->name_test1];
- mfg1 = &(vgi.mf);
- mf1 = 0;
+ if (workspace.secondary_domain_exists(intn1)) {
+ pctx1 = &(rmi.secondary_domain_infos.ctx);
+ } else {
+ tensor_to_adapt = true;
+ pctx1 = &(rmi.interpolate_infos[intn1].ctx);
+ if (workspace.variable_group_exists(pnode->name_test1)) {
+ ga_instruction_set::variable_group_info &vgi =
+ rmi.interpolate_infos[intn1].groups_info[pnode->name_test1];
+ mfg1 = &(vgi.mf);
+ mf1 = 0;
+ }
}
}
}
@@ -4781,14 +4785,18 @@ namespace getfem {
pctx2 = &(gis.ctx);
const std::string &intn2 = pnode->interpolate_name_test2;
if (intn2.size()) {
- tensor_to_adapt = true;
- pctx2 = &(rmi.interpolate_infos[intn2].ctx);
- if (workspace.variable_group_exists(pnode->name_test2)) {
- ga_instruction_set::variable_group_info &vgi =
- rmi.interpolate_infos[intn2].groups_info[pnode->name_test2];
- mfg2 = &(vgi.mf);
- mf2 = 0;
- }
+ if (workspace.secondary_domain_exists(intn2)) {
+ pctx2 = &(rmi.secondary_domain_infos.ctx);
+ } else {
+ tensor_to_adapt = true;
+ pctx2 = &(rmi.interpolate_infos[intn2].ctx);
+ if (workspace.variable_group_exists(pnode->name_test2)) {
+ ga_instruction_set::variable_group_info &vgi =
+ rmi.interpolate_infos[intn2].groups_info[pnode->name_test2];
+ mfg2 = &(vgi.mf);
+ mf2 = 0;
+ }
+ }
}
}
}
@@ -5005,6 +5013,25 @@ namespace getfem {
rmi.instructions.push_back(std::move(pgai));
break;
+ case GA_NODE_SECONDARY_DOMAIN_X:
+ case GA_NODE_SECONDARY_DOMAIN_NORMAL:
+ {
+ GMM_ASSERT1(!function_case,
+ "No use of Secondary_domain is allowed in functions");
+ auto psd = workspace.secondary_domain(pnode->interpolate_name);
+ size_type sddim = psd->mim().linked_mesh().dim();
+ if (pnode->tensor().size() != sddim)
+ pnode->init_vector_tensor(sddim);
+ if (pnode->node_type == GA_NODE_SECONDARY_DOMAIN_X)
+ pgai = std::make_shared<ga_instruction_X>
+ (pnode->tensor(), rmi.secondary_domain_infos.ctx);
+ else if (pnode->node_type == GA_NODE_SECONDARY_DOMAIN_NORMAL)
+ pgai = std::make_shared<ga_instruction_copy_Normal>
+ (pnode->tensor(), rmi.secondary_domain_infos.Normal);
+ rmi.instructions.push_back(std::move(pgai));
+ }
+ break;
+
case GA_NODE_VAL: case GA_NODE_GRAD:
case GA_NODE_HESS: case GA_NODE_DIVERG:
case GA_NODE_ELEMENTARY_VAL: case GA_NODE_ELEMENTARY_GRAD:
@@ -5282,6 +5309,118 @@ namespace getfem {
}
break;
+ case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
+ case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
+ {
+ GMM_ASSERT1(!function_case, "internal error");
+ const mesh_fem *mf = workspace.associated_mf(pnode->name);
+ const im_data *imd = workspace.associated_im_data(pnode->name);
+ const std::string &intn = pnode->interpolate_name;
+ auto &sdi = rmi.secondary_domain_infos;
+
+ fem_interpolation_context *pctx = &(sdi.ctx);
+ papprox_integration pai = sdi.pai;
+ psecondary_domain psd = workspace.secondary_domain(intn);
+
+ if (imd) {
+ pgai = std::make_shared<ga_instruction_extract_local_im_data>
+ (pnode->tensor(), *imd, workspace.value(pnode->name),
+ pai, *pctx, workspace.qdim(pnode->name));
+ rmi.instructions.push_back(std::move(pgai));
+ } else {
+ GMM_ASSERT1(mf, "Internal error");
+ GMM_ASSERT1(&(mf->linked_mesh()) == &(psd->mim().linked_mesh()),
+ "The finite element of variable " << pnode->name <<
+ " has to be defined on the same mesh than the "
+ "integration method or interpolation used on the "
+ "secondary domain");
+
+ // An instruction for extracting local dofs of the variable.
+ if (sdi.local_dofs.count(pnode->name) == 0) {
+ sdi.local_dofs[pnode->name] = base_vector(1);
+ extend_variable_in_gis(workspace, pnode->name, gis);
+ size_type qmult2 = mf->get_qdim();
+ if (qmult2 > 1 && !(mf->is_uniformly_vectorized()))
+ qmult2 = size_type(-1);
+ pgai = std::make_shared<ga_instruction_slice_local_dofs>
+ (*mf, *(gis.extended_vars[pnode->name]), *pctx,
+ sdi.local_dofs[pnode->name],
+ workspace.qdim(pnode->name) / mf->get_qdim(), qmult2);
+ rmi.elt_instructions.push_back(std::move(pgai));
+ }
+
+ // An instruction for pfp update
+ if (sdi.pfps.count(mf) == 0) {
+ sdi.pfps[mf] = 0;
+ pgai = std::make_shared<ga_instruction_update_pfp>
+ (*mf, sdi.pfps[mf], *pctx, gis.fp_pool);
+ if (mf->is_uniform())
+ rmi.begin_instructions.push_back(std::move(pgai));
+ else
+ rmi.instructions.push_back(std::move(pgai));
+ }
+
+ // An instruction for the base value
+ pgai = pga_instruction();
+ switch (pnode->node_type) {
+ case GA_NODE_SECONDARY_DOMAIN_VAL:
+ if (sdi.base.count(mf) == 0 ||
+ !(if_hierarchy.is_compatible(rmi.base_hierarchy[mf]))) {
+ rmi.base_hierarchy[mf].push_back(if_hierarchy);
+ pgai = std::make_shared<ga_instruction_val_base>
+ (sdi.base[mf], *pctx, *mf, sdi.pfps[mf]);
+ }
+ break;
+ case GA_NODE_SECONDARY_DOMAIN_GRAD:
+ case GA_NODE_SECONDARY_DOMAIN_DIVERG:
+ if (sdi.grad.count(mf) == 0 ||
+ !(if_hierarchy.is_compatible(rmi.grad_hierarchy[mf]))) {
+ rmi.grad_hierarchy[mf].push_back(if_hierarchy);
+ pgai = std::make_shared<ga_instruction_grad_base>
+ (sdi.grad[mf], *pctx, *mf, sdi.pfps[mf]);
+ }
+ break;
+ case GA_NODE_SECONDARY_DOMAIN_HESS:
+ if (sdi.hess.count(mf) == 0 ||
+ !(if_hierarchy.is_compatible(rmi.hess_hierarchy[mf]))) {
+ rmi.hess_hierarchy[mf].push_back(if_hierarchy);
+ pgai = std::make_shared<ga_instruction_hess_base>
+ (sdi.hess[mf], *pctx, *mf, sdi.pfps[mf]);
+ }
+ break;
+ default : GMM_ASSERT1(false, "Internal error");
+ }
+ if (pgai) rmi.instructions.push_back(std::move(pgai));
+
+ // The eval instruction
+ switch (pnode->node_type) {
+ case GA_NODE_SECONDARY_DOMAIN_VAL: // --> t(target_dim*Qmult)
+ pgai = std::make_shared<ga_instruction_val>
+ (pnode->tensor(), sdi.base[mf], sdi.local_dofs[pnode->name],
+ workspace.qdim(pnode->name));
+ break;
+ case GA_NODE_SECONDARY_DOMAIN_GRAD: // --> t(target_dim*Qmult,N)
+ pgai = std::make_shared<ga_instruction_grad>
+ (pnode->tensor(), sdi.grad[mf],
+ sdi.local_dofs[pnode->name], workspace.qdim(pnode->name));
+ break;
+ case GA_NODE_SECONDARY_DOMAIN_HESS: // --> t(target_dim*Qmult,N,N)
+ pgai = std::make_shared<ga_instruction_hess>
+ (pnode->tensor(), sdi.hess[mf],
+ sdi.local_dofs[pnode->name], workspace.qdim(pnode->name));
+ break;
+ case GA_NODE_SECONDARY_DOMAIN_DIVERG: // --> t(1)
+ pgai = std::make_shared<ga_instruction_diverg>
+ (pnode->tensor(), sdi.grad[mf],
+ sdi.local_dofs[pnode->name], workspace.qdim(pnode->name));
+ break;
+ default: break;
+ }
+ rmi.instructions.push_back(std::move(pgai));
+ }
+ }
+ break;
+
case GA_NODE_INTERPOLATE_VAL: case GA_NODE_INTERPOLATE_GRAD:
case GA_NODE_INTERPOLATE_HESS: case GA_NODE_INTERPOLATE_DIVERG:
{
@@ -5591,6 +5730,114 @@ namespace getfem {
}
break;
+ case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
+ case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
+ case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
+ case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
+ {
+ GMM_ASSERT1(!function_case, "internal error");
+ const mesh_fem *mf = workspace.associated_mf(pnode->name);
+ const std::string &intn = pnode->interpolate_name;
+ auto &sdi = rmi.secondary_domain_infos;
+
+ fem_interpolation_context *pctx = &(sdi.ctx);
+ papprox_integration pai = sdi.pai;
+ psecondary_domain psd = workspace.secondary_domain(intn);
+ if (mf) {
+ GMM_ASSERT1(&(mf->linked_mesh()) == &(psd->mim().linked_mesh()),
+ "The finite element of variable " << pnode->name <<
+ " and the applied integration method have to be"
+ " defined on the same mesh for secondary domain");
+
+ // An instruction for pfp update
+ if (sdi.pfps.count(mf) == 0) {
+ sdi.pfps[mf] = 0;
+ pgai = std::make_shared<ga_instruction_update_pfp>
+ (*mf, sdi.pfps[mf], *pctx, gis.fp_pool);
+ if (is_uniform)
+ rmi.begin_instructions.push_back(std::move(pgai));
+ else
+ rmi.instructions.push_back(std::move(pgai));
+ }
+
+ // An instruction for the base value
+ pgai = pga_instruction();
+ switch (pnode->node_type) {
+ case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
+ if (sdi.base.count(mf) == 0 ||
+ !(if_hierarchy.is_compatible(rmi.base_hierarchy[mf]))) {
+ rmi.base_hierarchy[mf].push_back(if_hierarchy);
+ pgai = std::make_shared<ga_instruction_val_base>
+ (sdi.base[mf], *pctx, *mf, sdi.pfps[mf]);
+ }
+ break;
+ case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
+ case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
+ if (sdi.grad.count(mf) == 0 ||
+ !(if_hierarchy.is_compatible(rmi.grad_hierarchy[mf]))) {
+ rmi.grad_hierarchy[mf].push_back(if_hierarchy);
+ pgai = std::make_shared<ga_instruction_grad_base>
+ (sdi.grad[mf], *pctx, *mf, sdi.pfps[mf]);
+ }
+ break;
+ case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
+ if (sdi.hess.count(mf) == 0 ||
+ !(if_hierarchy.is_compatible(rmi.hess_hierarchy[mf]))) {
+ rmi.hess_hierarchy[mf].push_back(if_hierarchy);
+ pgai = std::make_shared<ga_instruction_hess_base>
+ (sdi.hess[mf], *pctx, *mf, sdi.pfps[mf]);
+ }
+ break;
+ default : GMM_ASSERT1(false, "Internal error");
+ }
+ if (pgai) rmi.instructions.push_back(std::move(pgai));
+
+ // The copy of the real_base_value
+ switch(pnode->node_type) {
+ case GA_NODE_SECONDARY_DOMAIN_VAL_TEST:
+ // --> t(Qmult*ndof,Qmult*target_dim)
+ if (mf->get_qdim() > 1 && mf->is_uniformly_vectorized()) {
+ pnode->t.set_sparsity(1, mf->get_qdim());
+ tensor_to_clear = true;
+ pgai = std::make_shared<ga_instruction_copy_vect_val_base>
+ (pnode->tensor(), sdi.base[mf], mf->get_qdim());
+ } else {
+ pgai = std::make_shared<ga_instruction_copy_val_base>
+ (pnode->tensor(), sdi.base[mf], mf->get_qdim());
+ }
+ break;
+ case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
+ // --> t(Qmult*ndof,Qmult*target_dim,N)
+ if (mf->get_qdim() > 1 && mf->is_uniformly_vectorized()) {
+ pnode->t.set_sparsity(2, mf->get_qdim());
+ tensor_to_clear = true;
+ pgai = std::make_shared<ga_instruction_copy_vect_grad_base>
+ (pnode->tensor(), sdi.grad[mf], mf->get_qdim());
+ } else {
+ pgai = std::make_shared<ga_instruction_copy_grad_base>
+ (pnode->tensor(), sdi.grad[mf], mf->get_qdim());
+ }
+ break;
+ case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
+ // --> t(Qmult*ndof,Qmult*target_dim,N,N)
+ pgai = std::make_shared<ga_instruction_copy_hess_base>
+ (pnode->tensor(), sdi.hess[mf], mf->get_qdim());
+ if (mf->get_qdim() > 1 && mf->is_uniformly_vectorized())
+ pnode->t.set_sparsity(3, mf->get_qdim());
+ break;
+ case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
+ // --> t(Qmult*ndof)
+ pgai = std::make_shared<ga_instruction_copy_diverg_base>
+ (pnode->tensor(), sdi.grad[mf], mf->get_qdim());
+ break;
+ default: break;
+ }
+ if (pgai) rmi.instructions.push_back(std::move(pgai));
+ }
+ add_interval_to_gis(workspace, pnode->name, gis);
+ }
+ break;
+
case GA_NODE_INTERPOLATE_VAL_TEST: case GA_NODE_INTERPOLATE_GRAD_TEST:
case GA_NODE_INTERPOLATE_HESS_TEST: case GA_NODE_INTERPOLATE_DIVERG_TEST:
{
@@ -6275,7 +6522,7 @@ namespace getfem {
if (root) {
GMM_ASSERT1(!scalar || (root->tensor().size() == 1),
"The result of the given expression is not a scalar");
- ga_instruction_set::region_mim rm(td.mim, td.rg);
+ ga_instruction_set::region_mim rm(td.mim, td.rg, 0);
gis.whole_instructions[rm].m = td.m;
ga_if_hierarchy if_hierarchy;
ga_compile_node(root, workspace, gis,
@@ -6394,7 +6641,7 @@ namespace getfem {
pga_tree_node root = gis.trees.back().root;
if (root) {
// Compile tree
- ga_instruction_set::region_mim rm(td.mim, td.rg);
+ ga_instruction_set::region_mim rm(td.mim, td.rg, 0);
ga_instruction_set::region_mim_instructions &rmi
= gis.whole_instructions[rm];
rmi.m = td.m;
@@ -6445,7 +6692,10 @@ namespace getfem {
// Compile tree
// cout << "Will compile "; ga_print_node(root, cout); cout <<
endl;
- ga_instruction_set::region_mim rm(td.mim, td.rg);
+ psecondary_domain psd(0);
+ if (added_tree->secondary_domain.size())
+ psd = workspace.secondary_domain(added_tree->secondary_domain);
+ ga_instruction_set::region_mim rm(td.mim, td.rg, psd);
ga_instruction_set::region_mim_instructions &rmi
= gis.whole_instructions[rm];
rmi.m = td.m;
@@ -6488,7 +6738,9 @@ namespace getfem {
if (mf) {
const std::string &intn1 = root->interpolate_name_test1;
const gmm::sub_interval *Ir = 0, *In = 0;
- if (intn1.size() &&
+ bool secondary = intn1.size() &&
+ workspace.secondary_domain_exists(intn1);
+ if (intn1.size() && !secondary &&
workspace.variable_group_exists(root->name_test1)) {
ga_instruction_set::variable_group_info &vgi =
rmi.interpolate_infos[intn1]
@@ -6502,10 +6754,12 @@ namespace getfem {
In = &(workspace.interval_of_variable(root->name_test1));
}
fem_interpolation_context &ctx
- = intn1.size() ? rmi.interpolate_infos[intn1].ctx
- : gis.ctx;
+ = intn1.size() ?
+ (secondary ? rmi.secondary_domain_infos.ctx
+ : rmi.interpolate_infos[intn1].ctx) : gis.ctx;
bool interpolate
- = (!intn1.empty() && intn1.compare("neighbour_elt")!=0);
+ = (!intn1.empty() && intn1.compare("neighbour_elt")!=0 &&
+ !secondary);
pgai = std::make_shared<ga_instruction_fem_vector_assembly>
(root->tensor(), workspace.unreduced_vector(),
workspace.assembled_vector(), ctx, *Ir, *In, mf, mfg,
@@ -6528,15 +6782,23 @@ namespace getfem {
const mesh_fem **mfg1 = 0, **mfg2 = 0;
const std::string &intn1 = root->interpolate_name_test1;
const std::string &intn2 = root->interpolate_name_test2;
- fem_interpolation_context &ctx1
- = intn1.empty() ? gis.ctx
- : rmi.interpolate_infos[intn1].ctx;
- fem_interpolation_context &ctx2
- = intn2.empty() ? gis.ctx
- : rmi.interpolate_infos[intn2].ctx;
- bool interpolate
- = (!intn1.empty() && intn1.compare("neighbour_elt")!=0)
- || (!intn2.empty() && intn2.compare("neighbour_elt")!=0);
+ bool secondary1 = intn1.size() &&
+ workspace.secondary_domain_exists(intn1);
+ bool secondary2 = intn2.size() &&
+ workspace.secondary_domain_exists(intn2);
+ fem_interpolation_context &ctx1
+ = intn1.size() ?
+ (secondary1 ? rmi.secondary_domain_infos.ctx
+ : rmi.interpolate_infos[intn1].ctx) : gis.ctx;
+ fem_interpolation_context &ctx2
+ = intn2.size() ?
+ (secondary2 ? rmi.secondary_domain_infos.ctx
+ : rmi.interpolate_infos[intn2].ctx) : gis.ctx;
+ bool interpolate
+ = (!intn1.empty() && intn1.compare("neighbour_elt")!=0 &&
+ !secondary1) ||
+ (!intn2.empty() && intn2.compare("neighbour_elt")!=0 &&
+ !secondary2);
add_interval_to_gis(workspace, root->name_test1, gis);
add_interval_to_gis(workspace, root->name_test2, gis);
@@ -6544,7 +6806,7 @@ namespace getfem {
const gmm::sub_interval *Ir1 = 0, *In1 = 0, *Ir2 = 0, *In2=0;
const scalar_type *alpha1 = 0, *alpha2 = 0;
- if (!intn1.empty() &&
+ if (!intn1.empty() && !secondary1 &&
workspace.variable_group_exists(root->name_test1)) {
ga_instruction_set::variable_group_info &vgi =
rmi.interpolate_infos[intn1]
@@ -6560,7 +6822,7 @@ namespace getfem {
In1 = &(workspace.interval_of_variable(root->name_test1));
}
- if (!intn2.empty() &&
+ if (!intn2.empty() && !secondary2 &&
workspace.variable_group_exists(root->name_test2)) {
ga_instruction_set::variable_group_info &vgi =
rmi.interpolate_infos[intn2]
@@ -6627,6 +6889,7 @@ namespace getfem {
}
+
//=========================================================================
// Execution of a compiled set of assembly terms
//=========================================================================
@@ -6755,15 +7018,16 @@ namespace getfem {
}
void ga_exec(ga_instruction_set &gis, ga_workspace &workspace) {
- base_matrix G;
+ base_matrix G1, G2;
base_small_vector un;
- scalar_type J(0);
+ scalar_type J1(0), J2(0);
for (const std::string &t : gis.transformations)
workspace.interpolate_transformation(t)->init(workspace);
- for (const auto &instr : gis.whole_instructions) {
+ for (auto &instr : gis.whole_instructions) {
const getfem::mesh_im &mim = *(instr.first.mim());
+ psecondary_domain psd = instr.first.psd();
const getfem::mesh &m = *(instr.second.m);
GMM_ASSERT1(&m == &(mim.linked_mesh()), "Incompatibility of meshes");
const ga_instruction_list &gilb = instr.second.begin_instructions;
@@ -6780,98 +7044,274 @@ namespace getfem {
// for (size_type j = 0; j < gil.size(); ++j)
// cout << typeid(*(gil[j])).name() << endl;
- const mesh_region ®ion = *(instr.first.region());
-
- // iteration on elements (or faces of elements)
- size_type old_cv = size_type(-1);
- bgeot::pgeometric_trans pgt = 0, pgt_old = 0;
- pintegration_method pim = 0;
- papprox_integration pai = 0;
- bgeot::pstored_point_tab pspt = 0, old_pspt = 0;
- bgeot::pgeotrans_precomp pgp = 0;
- bool first_gp = true;
- for (getfem::mr_visitor v(region, m, true); !v.finished(); ++v) {
- if (mim.convex_index().is_in(v.cv())) {
- // cout << "proceed with elt " << v.cv() << " face " << v.f() <<
endl;
- if (v.cv() != old_cv) {
- pgt = m.trans_of_convex(v.cv());
- pim = mim.int_method_of_element(v.cv());
- m.points_of_convex(v.cv(), G);
-
- if (pim->type() == IM_NONE) continue;
- GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods cannot
"
- "be used in high level generic assembly");
- pai = pim->approx_method();
- pspt = pai->pintegration_points();
- if (pspt->size()) {
- if (pgp && gis.pai == pai && pgt_old == pgt) {
- gis.ctx.change(pgp, 0, 0, G, v.cv(), v.f());
- } else {
- if (pai->is_built_on_the_fly()) {
- gis.ctx.change(pgt, 0, (*pspt)[0], G, v.cv(), v.f());
- pgp = 0;
- } else {
- pgp = gis.gp_pool(pgt, pspt);
- gis.ctx.change(pgp, 0, 0, G, v.cv(), v.f());
- }
- pgt_old = pgt; gis.pai = pai;
- }
- if (gis.need_elt_size)
- gis.elt_size = convex_radius_estimate(pgt, G)*scalar_type(2);
- }
- old_cv = v.cv();
- } else {
- if (pim->type() == IM_NONE) continue;
- gis.ctx.set_face_num(v.f());
- }
- if (pspt != old_pspt) { first_gp = true; old_pspt = pspt; }
- if (pspt->size()) {
- // iterations on Gauss points
- gis.nbpt = pai->nb_points_on_convex();
- size_type first_ind = 0;
- if (v.f() != short_type(-1)) {
- gis.nbpt = pai->nb_points_on_face(v.f());
- first_ind = pai->ind_first_point_on_face(v.f());
- }
- for (gis.ipt = 0; gis.ipt < gis.nbpt; ++(gis.ipt)) {
- if (pgp) gis.ctx.set_ii(first_ind+gis.ipt);
- else gis.ctx.set_xref((*pspt)[first_ind+gis.ipt]);
- if (gis.ipt == 0 || !(pgt->is_linear())) {
- J = gis.ctx.J();
- // Computation of unit normal vector in case of a boundary
- if (v.f() != short_type(-1)) {
- gis.Normal.resize(G.nrows());
- un.resize(pgt->dim());
- gmm::copy(pgt->normals()[v.f()], un);
- gmm::mult(gis.ctx.B(), un, gis.Normal);
- scalar_type nup = gmm::vect_norm2(gis.Normal);
- J *= nup;
- gmm::scale(gis.Normal, 1.0/nup);
- gmm::clean(gis.Normal, 1e-13);
- } else gis.Normal.resize(0);
- }
- auto ipt_coeff = pai->coeff(first_ind+gis.ipt);
- gis.coeff = J * ipt_coeff;
- bool enable_ipt = (gmm::abs(ipt_coeff) > 0.0 ||
- workspace.include_empty_int_points());
- if (!enable_ipt) gis.coeff = scalar_type(0);
- if (first_gp) {
- for (size_type j = 0; j < gilb.size(); ++j) j+=gilb[j]->exec();
- first_gp = false;
- }
- if (gis.ipt == 0) {
- for (size_type j = 0; j < gile.size(); ++j) j+=gile[j]->exec();
- }
- if (enable_ipt || gis.ipt == 0 || gis.ipt == gis.nbpt-1) {
- for (size_type j = 0; j < gil.size(); ++j) j+=gil[j]->exec();
- }
- GA_DEBUG_INFO("");
- }
- }
- }
+ if (!psd) { // standard integration on a single domain
+
+ const mesh_region ®ion = *(instr.first.region());
+
+ // iteration on elements (or faces of elements)
+ size_type old_cv = size_type(-1);
+ bgeot::pgeometric_trans pgt = 0, pgt_old = 0;
+ pintegration_method pim = 0;
+ papprox_integration pai = 0;
+ bgeot::pstored_point_tab pspt = 0, old_pspt = 0;
+ bgeot::pgeotrans_precomp pgp = 0;
+ bool first_gp = true;
+ for (getfem::mr_visitor v(region, m, true); !v.finished(); ++v) {
+ if (mim.convex_index().is_in(v.cv())) {
+ // cout << "proceed with elt " << v.cv() << " face " << v.f()<<endl;
+ if (v.cv() != old_cv) {
+ pgt = m.trans_of_convex(v.cv());
+ pim = mim.int_method_of_element(v.cv());
+ m.points_of_convex(v.cv(), G1);
+
+ if (pim->type() == IM_NONE) continue;
+ GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods "
+ "cannot be used in high level generic assembly");
+ pai = pim->approx_method();
+ pspt = pai->pintegration_points();
+ if (pspt->size()) {
+ if (pgp && gis.pai == pai && pgt_old == pgt) {
+ gis.ctx.change(pgp, 0, 0, G1, v.cv(), v.f());
+ } else {
+ if (pai->is_built_on_the_fly()) {
+ gis.ctx.change(pgt, 0, (*pspt)[0], G1, v.cv(), v.f());
+ pgp = 0;
+ } else {
+ pgp = gis.gp_pool(pgt, pspt);
+ gis.ctx.change(pgp, 0, 0, G1, v.cv(), v.f());
+ }
+ pgt_old = pgt; gis.pai = pai;
+ }
+ if (gis.need_elt_size)
+ gis.elt_size = convex_radius_estimate(pgt, G1)*scalar_type(2);
+ }
+ old_cv = v.cv();
+ } else {
+ if (pim->type() == IM_NONE) continue;
+ gis.ctx.set_face_num(v.f());
+ }
+ if (pspt != old_pspt) { first_gp = true; old_pspt = pspt; }
+ if (pspt->size()) {
+ // iterations on Gauss points
+ size_type first_ind = 0;
+ if (v.f() != short_type(-1)) {
+ gis.nbpt = pai->nb_points_on_face(v.f());
+ first_ind = pai->ind_first_point_on_face(v.f());
+ } else {
+ gis.nbpt = pai->nb_points_on_convex();
+ }
+ for (gis.ipt = 0; gis.ipt < gis.nbpt; ++(gis.ipt)) {
+ if (pgp) gis.ctx.set_ii(first_ind+gis.ipt);
+ else gis.ctx.set_xref((*pspt)[first_ind+gis.ipt]);
+ if (gis.ipt == 0 || !(pgt->is_linear())) {
+ J1 = gis.ctx.J();
+ // Computation of unit normal vector in case of a boundary
+ if (v.f() != short_type(-1)) {
+ gis.Normal.resize(G1.nrows());
+ un.resize(pgt->dim());
+ gmm::copy(pgt->normals()[v.f()], un);
+ gmm::mult(gis.ctx.B(), un, gis.Normal);
+ scalar_type nup = gmm::vect_norm2(gis.Normal);
+ J1 *= nup;
+ gmm::scale(gis.Normal, 1.0/nup);
+ gmm::clean(gis.Normal, 1e-13);
+ } else gis.Normal.resize(0);
+ }
+ auto ipt_coeff = pai->coeff(first_ind+gis.ipt);
+ gis.coeff = J1 * ipt_coeff;
+ bool enable_ipt = (gmm::abs(ipt_coeff) > 0.0 ||
+ workspace.include_empty_int_points());
+ if (!enable_ipt) gis.coeff = scalar_type(0);
+ if (first_gp) {
+ for (size_type j=0; j < gilb.size(); ++j) j+=gilb[j]->exec();
+ first_gp = false;
+ }
+ if (gis.ipt == 0) {
+ for (size_type j=0; j < gile.size(); ++j) j+=gile[j]->exec();
+ }
+ if (enable_ipt || gis.ipt == 0 || gis.ipt == gis.nbpt-1) {
+ for (size_type j=0; j < gil.size(); ++j) j+=gil[j]->exec();
+ }
+ GA_DEBUG_INFO("");
+ }
+ }
+ }
+ }
+ GA_DEBUG_INFO("-----------------------------");
+
+ } else { // Integration on the product of two domains (secondary domain)
+
+ auto &sdi = instr.second.secondary_domain_infos;
+ const mesh_region ®ion1 = *(instr.first.region());
+
+ // iteration on elements (or faces of elements)
+ size_type old_cv1=size_type(-1), old_cv2=size_type(-1);
+ size_type nbpt1 = 0, nbpt2 = 0;
+ bgeot::pgeometric_trans pgt1 = 0, pgt1_old = 0, pgt2 = 0, pgt2_old = 0;
+ pintegration_method pim1 = 0, pim2 = 0;
+ papprox_integration pai1 = 0, pai2 = 0;
+ bgeot::pstored_point_tab pspt1=0, old_pspt1=0, pspt2=0, old_pspt2=0;
+ bgeot::pgeotrans_precomp pgp1 = 0, pgp2 = 0;
+ bool first_gp = true;
+ for (getfem::mr_visitor v1(region1, m, true); !v1.finished(); ++v1) {
+ if (mim.convex_index().is_in(v1.cv())) {
+ // cout << "proceed with elt " << v1.cv()<<" face " << v1.f()<<endl;
+ if (v1.cv() != old_cv1) {
+ pgt1 = m.trans_of_convex(v1.cv());
+ pim1 = mim.int_method_of_element(v1.cv());
+ m.points_of_convex(v1.cv(), G1);
+
+ if (pim1->type() == IM_NONE) continue;
+ GMM_ASSERT1(pim1->type() == IM_APPROX, "Sorry, exact methods "
+ "cannot be used in high level generic assembly");
+ pai1 = pim1->approx_method();
+ pspt1 = pai1->pintegration_points();
+ if (pspt1->size()) {
+ if (pgp1 && gis.pai == pai1 && pgt1_old == pgt1) {
+ gis.ctx.change(pgp1, 0, 0, G1, v1.cv(), v1.f());
+ } else {
+ if (pai1->is_built_on_the_fly()) {
+ gis.ctx.change(pgt1, 0, (*pspt1)[0], G1, v1.cv(), v1.f());
+ pgp1 = 0;
+ } else {
+ pgp1 = gis.gp_pool(pgt1, pspt1);
+ gis.ctx.change(pgp1, 0, 0, G1, v1.cv(), v1.f());
+ }
+ pgt1_old = pgt1; gis.pai = pai1;
+ }
+ if (gis.need_elt_size)
+ gis.elt_size = convex_radius_estimate(pgt1,G1)*scalar_type(2);
+ }
+ old_cv1 = v1.cv();
+ } else {
+ if (pim1->type() == IM_NONE) continue;
+ gis.ctx.set_face_num(v1.f());
+ }
+ if (pspt1 != old_pspt1) { first_gp = true; old_pspt1 = pspt1; }
+ if (pspt1->size()) {
+ // iterations on Gauss points
+ size_type first_ind1 = 0;
+ if (v1.f() != short_type(-1)) {
+ nbpt1 = pai1->nb_points_on_face(v1.f());
+ first_ind1 = pai1->ind_first_point_on_face(v1.f());
+ } else {
+ nbpt1 = pai1->nb_points_on_convex();
+ }
+
+ const mesh &m2 = psd->mim().linked_mesh();
+ const mesh_region ®ion2 = psd->give_region(m, v1.cv(), v1.f());
+ for (getfem::mr_visitor v2(region2, m2, true);
+ !v2.finished(); ++v2) {
+ if (v2.cv() != old_cv2) {
+ pgt2 = m2.trans_of_convex(v2.cv());
+ pim2 = psd->mim().int_method_of_element(v2.cv());
+ m2.points_of_convex(v2.cv(), G2);
+
+ if (pim2->type() == IM_NONE) continue;
+ GMM_ASSERT1(pim2->type() == IM_APPROX, "Sorry, exact methods "
+ "cannot be used in high level generic assembly");
+ pai2 = pim2->approx_method();
+ pspt2 = pai2->pintegration_points();
+ if (pspt2->size()) {
+ if (pgp2 && sdi.pai == pai2 && pgt2_old == pgt2) {
+ sdi.ctx.change(pgp2, 0, 0, G2, v2.cv(), v2.f());
+ } else {
+ if (pai2->is_built_on_the_fly()) {
+ sdi.ctx.change(pgt2, 0, (*pspt2)[0], G2,v2.cv(),v2.f());
+ pgp2 = 0;
+ } else {
+ pgp2 = gis.gp_pool(pgt2, pspt2);
+ sdi.ctx.change(pgp2, 0, 0, G2, v2.cv(), v2.f());
+ }
+ pgt2_old = pgt2; sdi.pai = pai2;
+ }
+ }
+ old_cv2 = v2.cv();
+ } else {
+ if (pim2->type() == IM_NONE) continue;
+ sdi.ctx.set_face_num(v2.f());
+ }
+ if (pspt2 != old_pspt2) { first_gp = true; old_pspt2 = pspt2; }
+ if (pspt2->size()) {
+ // iterations on Gauss points
+ size_type first_ind2 = 0;
+ if (v2.f() != short_type(-1)) {
+ nbpt2 = pai2->nb_points_on_face(v2.f());
+ first_ind2 = pai2->ind_first_point_on_face(v2.f());
+ } else {
+ nbpt2 = gis.nbpt = pai2->nb_points_on_convex();
+ }
+ gis.nbpt = nbpt1 * nbpt2;
+ gis.ipt = 0;
+ for (size_type ipt1=0; ipt1 < nbpt1; ++ipt1) {
+ for (size_type ipt2=0; ipt2 < nbpt2; ++ipt2, ++(gis.ipt)) {
+
+ if (pgp1) gis.ctx.set_ii(first_ind1+ipt1);
+ else gis.ctx.set_xref((*pspt1)[first_ind1+ipt1]);
+ if (pgp2) sdi.ctx.set_ii(first_ind2+ipt2);
+ else sdi.ctx.set_xref((*pspt2)[first_ind2+ipt2]);
+
+ if (gis.ipt == 0 || !(pgt1->is_linear())) {
+ J1 = gis.ctx.J();
+ if (v1.f() != short_type(-1)) {
+ gis.Normal.resize(G1.nrows());
+ un.resize(pgt1->dim());
+ gmm::copy(pgt1->normals()[v1.f()], un);
+ gmm::mult(gis.ctx.B(), un, gis.Normal);
+ scalar_type nup = gmm::vect_norm2(gis.Normal);
+ J1 *= nup;
+ gmm::scale(gis.Normal, 1.0/nup);
+ gmm::clean(gis.Normal, 1e-13);
+ } else gis.Normal.resize(0);
+ }
+
+ if (gis.ipt == 0 || !(pgt2->is_linear())) {
+ J2 = sdi.ctx.J();
+ if (v2.f() != short_type(-1)) {
+ sdi.Normal.resize(G2.nrows());
+ un.resize(pgt2->dim());
+ gmm::copy(pgt2->normals()[v2.f()], un);
+ gmm::mult(sdi.ctx.B(), un, sdi.Normal);
+ scalar_type nup = gmm::vect_norm2(sdi.Normal);
+ J2 *= nup;
+ gmm::scale(sdi.Normal, 1.0/nup);
+ gmm::clean(sdi.Normal, 1e-13);
+ } else sdi.Normal.resize(0);
+ }
+
+ auto ipt_coeff = pai1->coeff(first_ind1+ipt1)
+ * pai2->coeff(first_ind2+ipt2);
+ gis.coeff = J1 * J2 * ipt_coeff;
+ bool enable_ipt = (gmm::abs(ipt_coeff) > 0.0 ||
+ workspace.include_empty_int_points());
+ if (!enable_ipt) gis.coeff = scalar_type(0);
+
+ if (first_gp) {
+ for (size_type j=0; j < gilb.size(); ++j)
+ j+=gilb[j]->exec();
+ first_gp = false;
+ }
+ if (gis.ipt == 0) {
+ for (size_type j=0; j < gile.size(); ++j)
+ j+=gile[j]->exec();
+ }
+ if (enable_ipt || gis.ipt == 0 || gis.ipt == gis.nbpt-1) {
+ for (size_type j=0; j < gil.size(); ++j)
+ j+=gil[j]->exec();
+ }
+ GA_DEBUG_INFO("");
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ GA_DEBUG_INFO("-----------------------------");
}
- GA_DEBUG_INFO("-----------------------------");
+
}
+
for (const std::string &t : gis.transformations)
workspace.interpolate_transformation(t)->finalize();
}
diff --git a/src/getfem_generic_assembly_interpolation.cc
b/src/getfem_generic_assembly_interpolation.cc
index cbbee46..765d446 100644
--- a/src/getfem_generic_assembly_interpolation.cc
+++ b/src/getfem_generic_assembly_interpolation.cc
@@ -932,4 +932,39 @@ namespace getfem {
->set_correspondance(elt_corr);
}
+
+ //=========================================================================
+ // Secondary domains
+ //=========================================================================
+
+
+ class standard_secondary_domain : public virtual_secondary_domain {
+
+ public:
+
+ virtual const mesh_region &give_region(const mesh &,
+ size_type, short_type) const
+ { return region; }
+ // virtual void init(const ga_workspace &workspace) const = 0;
+ // virtual void finalize() const = 0;
+
+ standard_secondary_domain(const mesh_im &mim__, const mesh_region ®ion_)
+ : virtual_secondary_domain(mim__, region_) {}
+ };
+
+ void add_standard_secondary_domain
+ (model &md, const std::string &name, const mesh_im &mim,
+ const mesh_region &rg) {
+ psecondary_domain p = std::make_shared<standard_secondary_domain>(mim, rg);
+ md.add_secondary_domain(name, p);
+ }
+
+ void add_standard_secondary_domain
+ (ga_workspace &workspace, const std::string &name, const mesh_im &mim,
+ const mesh_region &rg) {
+ psecondary_domain p = std::make_shared<standard_secondary_domain>(mim, rg);
+ workspace.add_secondary_domain(name, p);
+ }
+
+
} /* end of namespace */
diff --git a/src/getfem_generic_assembly_semantic.cc
b/src/getfem_generic_assembly_semantic.cc
index 30e87ac..8d05147 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -518,21 +518,33 @@ namespace getfem {
}
break;
- case GA_NODE_INTERPOLATE: case GA_NODE_SECONDARY_DOMAIN:
+ case GA_NODE_SECONDARY_DOMAIN:
+ pnode->interpolate_name = tree.secondary_domain;
+ if (tree.secondary_domain.size() == 0)
+ ga_throw_error(pnode->expr, pnode->pos, "Secondary domain used "
+ "in a single domain term.");
+ // continue with what follows
+ case GA_NODE_INTERPOLATE:
if (pnode->name.compare("X") == 0) {
- if (pnode->node_type == GA_NODE_INTERPOLATE)
+ if (pnode->node_type == GA_NODE_INTERPOLATE) {
pnode->node_type = GA_NODE_INTERPOLATE_X;
- else
- pnode->node_type = GA_NODE_SECONDARY_DOMAIN_X;
- pnode->init_vector_tensor(meshdim);
+ pnode->init_vector_tensor(meshdim);
+ } else {
+ auto psd = workspace.secondary_domain(tree.secondary_domain);
+ pnode->node_type = GA_NODE_SECONDARY_DOMAIN_X;
+ pnode->init_vector_tensor(psd->mim().linked_mesh().dim());
+ }
break;
}
if (pnode->name.compare("Normal") == 0) {
- if (pnode->node_type == GA_NODE_INTERPOLATE)
+ if (pnode->node_type == GA_NODE_INTERPOLATE) {
pnode->node_type = GA_NODE_INTERPOLATE_NORMAL;
- else
+ pnode->init_vector_tensor(meshdim);
+ } else {
+ auto psd = workspace.secondary_domain(tree.secondary_domain);
pnode->node_type = GA_NODE_SECONDARY_DOMAIN_NORMAL;
- pnode->init_vector_tensor(meshdim);
+ pnode->init_vector_tensor(psd->mim().linked_mesh().dim());
+ }
break;
}
// else continue with what follows
@@ -2563,6 +2575,7 @@ namespace getfem {
size_type ref_elt_dim,
bool eval_fixed_size,
bool ignore_X, int option) {
+ // cout << "Begin semantic anaylsis" << endl;
GMM_ASSERT1(predef_operators_nonlinear_elasticity_initialized &&
predef_operators_plasticity_initialized &&
predef_operators_contact_initialized, "Internal error");
@@ -2583,6 +2596,7 @@ namespace getfem {
tree.clear();
}
ga_valid_operand(tree.root);
+ // cout << "end of semantic anaylsis" << endl;
}
diff --git a/src/getfem_generic_assembly_tree.cc
b/src/getfem_generic_assembly_tree.cc
index 7aaa926..caffa4f 100644
--- a/src/getfem_generic_assembly_tree.cc
+++ b/src/getfem_generic_assembly_tree.cc
@@ -1165,7 +1165,7 @@ namespace getfem {
else if (is_elementary)
str << "," << pnode->elementary_name << ")";
else if (is_secondary)
- str << "," << pnode->interpolate_name << ")";
+ str << ")";
else if (is_xfem_plus || is_xfem_minus)
str << ")";
@@ -1666,18 +1666,7 @@ namespace getfem {
"variable, test function, X or Normal.");
tree.current_node->name = std::string(&((*expr)[token_pos]),
token_length);
-
- t_type = ga_get_token(*expr, pos, token_pos, token_length);
- if (t_type != GA_COMMA)
- ga_throw_error(expr, pos, "Bad format for Secondary_domain "
- "arguments.");
- t_type = ga_get_token(*expr, pos, token_pos, token_length);
- if (t_type != GA_NAME)
- ga_throw_error(expr, pos,
- "Second argument of Secondary_domain should be a "
- "transformation name.");
- tree.current_node->interpolate_name
- = std::string(&((*expr)[token_pos]), token_length);
+ tree.current_node->interpolate_name = tree.secondary_domain;
t_type = ga_get_token(*expr, pos, token_pos, token_length);
if (t_type != GA_RPAR)
ga_throw_error(expr, pos-1, "Missing a parenthesis after "
diff --git a/src/getfem_generic_assembly_workspace.cc
b/src/getfem_generic_assembly_workspace.cc
index 193afc0..03ba618 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -398,7 +398,6 @@ namespace getfem {
bool function_expr, size_type for_interpolation,
const std::string varname_interpolation) {
if (tree.root) {
-
// Eliminate the term if it corresponds to disabled variables
if ((tree.root->test_function_type >= 1 &&
is_disabled_variable(tree.root->name_test1)) ||
@@ -428,6 +427,7 @@ namespace getfem {
bool found = false;
for (size_type i = 0; i < trees.size(); ++i) {
if (trees[i].mim == &mim && trees[i].m == &m &&
+ trees[i].secondary_domain.compare(tree.secondary_domain) == 0 &&
trees[i].order == order &&
trees[i].name_test1.compare(tree.root->name_test1) == 0 &&
trees[i].interpolate_name_test1.compare
@@ -455,6 +455,7 @@ namespace getfem {
trees.push_back(tree_description());
trees.back().mim = &mim; trees.back().m = &m;
trees.back().rg = &rg;
+ trees.back().secondary_domain = tree.secondary_domain;
trees.back().ptree = new ga_tree;
trees.back().ptree->swap(tree);
pga_tree_node root = trees.back().ptree->root;
@@ -506,13 +507,19 @@ namespace getfem {
size_type ga_workspace::add_expression(const std::string &expr,
const mesh_im &mim,
const mesh_region &rg_,
- size_type add_derivative_order) {
+ size_type add_derivative_order,
+ const std::string &secondary_dom) {
const mesh_region &rg = register_region(mim.linked_mesh(), rg_);
// cout << "adding expression " << expr << endl;
GA_TIC;
size_type max_order = 0;
std::vector<ga_tree> ltrees(1);
ga_read_string(expr, ltrees[0], macro_dictionnary());
+ if (secondary_dom.size()) {
+ GMM_ASSERT1(secondary_domain_exists(secondary_dom),
+ "Unknow secondary domain " << secondary_dom);
+ ltrees[0].secondary_domain = secondary_dom;
+ }
// cout << "read : " << ga_tree_to_string(ltrees[0]) << endl;
ga_semantic_analysis(ltrees[0], *this, mim.linked_mesh(),
ref_elt_dim_of_mesh(mim.linked_mesh()),
diff --git a/src/getfem_models.cc b/src/getfem_models.cc
index f7eaf97..afa4f71 100644
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@ -525,7 +525,8 @@ namespace getfem {
{
ga_workspace workspace(*this);
for (const auto &ge : generic_expressions)
- workspace.add_expression(ge.expr, ge.mim, ge.region);
+ workspace.add_expression(ge.expr, ge.mim, ge.region,
+ 2, ge.secondary_domain);
workspace.set_assembled_matrix(distro_rTM);
workspace.assembly(2);
});
@@ -2362,7 +2363,8 @@ namespace getfem {
(*this, ib, brick.vlist, brick.dlist);
ga_workspace workspace(*this);
- size_type order = workspace.add_expression(expr, dummy_mim, region);
+ size_type order = workspace.add_expression
+ (expr, dummy_mim, region);
GMM_ASSERT1(order <= 1, "Wrong order for a Neumann term");
expr = workspace.extract_Neumann_term(varname);
if (expr.size()) {
@@ -2706,7 +2708,8 @@ namespace getfem {
(ad.varname, ad.expr, ad.region, ad.order, ad.before);
for (const auto &ge : generic_expressions)
- workspace.add_expression(ge.expr, ge.mim, ge.region);
+ workspace.add_expression(ge.expr, ge.mim, ge.region,
+ 2, ge.secondary_domain);
if (version & BUILD_RHS) {
if (is_complex()) {
@@ -3275,6 +3278,7 @@ model_complex_plain_vector &
std::string expr, directvarname, directdataname;
model::varnamelist vl_test1;
+ std::string secondary_domain;
void asm_real_tangent_terms(const model &md, size_type /* ib */,
const model::varnamelist &,
@@ -3298,7 +3302,7 @@ model_complex_plain_vector &
size_type nbgdof = md.nb_dof();
ga_workspace workspace(md, true);
GMM_TRACE2(name << ": generic source term assembly");
- workspace.add_expression(expr, *(mims[0]), region);
+ workspace.add_expression(expr, *(mims[0]), region, 1,
secondary_domain);
model::varnamelist vlmd; md.variable_list(vlmd);
for (size_type i = 0; i < vlmd.size(); ++i)
if (md.is_disabled_variable(vlmd[i]))
@@ -3343,8 +3347,9 @@ model_complex_plain_vector &
std::string brickname,
const model::varnamelist &vl_test1_,
const std::string &directvarname_,
- const std::string &directdataname_)
- : vl_test1(vl_test1_) {
+ const std::string &directdataname_,
+ const std::string &secdom)
+ : vl_test1(vl_test1_), secondary_domain(secdom) {
if (brickname.size() == 0)
brickname = "Generic source term assembly brick";
expr = expr_;
@@ -3357,7 +3362,8 @@ model_complex_plain_vector &
};
- static bool check_compatibility_vl_test(model &md,const model::varnamelist
vl_test) {
+ static bool check_compatibility_vl_test(model &md,
+ const model::varnamelist vl_test) {
model::varnamelist org;
for (size_type i = 0; i < vl_test.size(); ++i) {
if (md.is_affine_dependent_variable(vl_test[i]))
@@ -3369,13 +3375,15 @@ model_complex_plain_vector &
return true;
}
- size_type add_source_term
+ size_type add_source_term_
(model &md, const mesh_im &mim, const std::string &expr, size_type region,
std::string brickname, std::string directvarname,
- const std::string &directdataname, bool return_if_nonlin) {
+ const std::string &directdataname, bool return_if_nonlin,
+ const std::string &secondary_domain) {
ga_workspace workspace(md);
- size_type order = workspace.add_expression(expr, mim, region);
+ size_type order = workspace.add_expression(expr, mim, region, 1,
+ secondary_domain);
GMM_ASSERT1(order <= 1, "Wrong order for a source term");
model::varnamelist vl, vl_test1, vl_test2, dl;
bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 1);
@@ -3392,7 +3400,8 @@ model_complex_plain_vector &
} else directvarname = "";
pbrick pbr = std::make_shared<gen_source_term_assembly_brick>
- (expr, brickname, vl_test1, directvarname, directdataname);
+ (expr, brickname, vl_test1, directvarname, directdataname,
+ secondary_domain);
model::termlist tl;
for (size_type i = 0; i < vl_test1.size(); ++i)
@@ -3402,7 +3411,24 @@ model_complex_plain_vector &
return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
}
+
+ size_type add_source_term
+ (model &md, const mesh_im &mim, const std::string &expr, size_type region,
+ std::string brickname, std::string directvarname,
+ const std::string &directdataname, bool return_if_nonlin) {
+ return add_source_term_(md, mim, expr, region, brickname, directvarname,
+ directdataname, return_if_nonlin, "");
+ }
+ size_type add_twodomain_source_term
+ (model &md, const mesh_im &mim, const std::string &expr, size_type region,
+ const std::string &secondary_domain,
+ std::string brickname, std::string directvarname,
+ const std::string &directdataname, bool return_if_nonlin) {
+ return add_source_term_(md, mim, expr, region, brickname, directvarname,
+ directdataname, return_if_nonlin, secondary_domain);
+ }
+
// ----------------------------------------------------------------------
//
// Linear generic assembly brick
@@ -3414,6 +3440,7 @@ model_complex_plain_vector &
std::string expr;
bool is_lower_dim;
model::varnamelist vl_test1, vl_test2;
+ std::string secondary_domain;
virtual void asm_real_tangent_terms(const model &md, size_type ib,
const model::varnamelist &/* vl */,
@@ -3437,7 +3464,7 @@ model_complex_plain_vector &
if (recompute_matrix) {
size_type nbgdof = md.nb_dof();
ga_workspace workspace(md, true);
- workspace.add_expression(expr, *(mims[0]), region);
+ workspace.add_expression(expr, *(mims[0]), region, 2,
secondary_domain);
model::varnamelist vlmd; md.variable_list(vlmd);
for (size_type i = 0; i < vlmd.size(); ++i)
if (md.is_disabled_variable(vlmd[i]))
@@ -3469,8 +3496,9 @@ model_complex_plain_vector &
bool is_sym,
bool is_coer, std::string brickname,
const model::varnamelist &vl_test1_,
- const model::varnamelist &vl_test2_)
- : vl_test1(vl_test1_), vl_test2(vl_test2_) {
+ const model::varnamelist &vl_test2_,
+ const std::string &secdom)
+ : vl_test1(vl_test1_), vl_test2(vl_test2_), secondary_domain(secdom) {
if (brickname.size() == 0) brickname = "Generic linear assembly brick";
expr = expr_;
is_lower_dim = mim.is_lower_dimensional();
@@ -3506,13 +3534,16 @@ model_complex_plain_vector &
return true;
}
- size_type add_linear_term
+
+
+ size_type add_linear_term_
(model &md, const mesh_im &mim, const std::string &expr, size_type region,
bool is_sym, bool is_coercive, std::string brickname,
- bool return_if_nonlin) {
+ bool return_if_nonlin, const std::string &secondary_domain) {
ga_workspace workspace(md, true);
- size_type order = workspace.add_expression(expr, mim, region);
+ size_type order = workspace.add_expression(expr, mim, region,
+ 2, secondary_domain);
model::varnamelist vl, vl_test1, vl_test2, dl;
bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
@@ -3522,8 +3553,9 @@ model_complex_plain_vector &
std::string const_expr= workspace.extract_constant_term(mim.linked_mesh());
if (const_expr.size()) {
- add_source_term_generic_assembly_brick
- (md, mim, const_expr, region, brickname+" (source term)");
+ add_source_term_
+ (md, mim, const_expr, region, brickname+" (source term)",
+ "", "", false, secondary_domain);
}
// GMM_ASSERT1(order <= 1,
@@ -3535,7 +3567,8 @@ model_complex_plain_vector &
if (vl_test1.size()) {
pbrick pbr = std::make_shared<gen_linear_assembly_brick>
- (expr, mim, is_sym, is_coercive, brickname, vl_test1, vl_test2);
+ (expr, mim, is_sym, is_coercive, brickname, vl_test1, vl_test2,
+ secondary_domain);
model::termlist tl;
for (size_type i = 0; i < vl_test1.size(); ++i)
tl.push_back(model::term_description(vl_test1[i], vl_test2[i], false));
@@ -3545,6 +3578,22 @@ model_complex_plain_vector &
return size_type(-1);
}
+ size_type add_linear_term
+ (model &md, const mesh_im &mim, const std::string &expr, size_type region,
+ bool is_sym, bool is_coercive, std::string brickname,
+ bool return_if_nonlin) {
+ return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
+ brickname, return_if_nonlin, "");
+ }
+
+ size_type add_linear_twodomain_term
+ (model &md, const mesh_im &mim, const std::string &expr, size_type region,
+ const std::string &secondary_domain, bool is_sym, bool is_coercive,
+ std::string brickname, bool return_if_nonlin) {
+ return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
+ brickname, return_if_nonlin, secondary_domain);
+ }
+
// ----------------------------------------------------------------------
//
@@ -3556,6 +3605,7 @@ model_complex_plain_vector &
std::string expr;
bool is_lower_dim;
+ std::string secondary_domain;
virtual void real_post_assembly_in_serial(const model &md, size_type ,
const model::varnamelist &,
@@ -3569,7 +3619,7 @@ model_complex_plain_vector &
GMM_ASSERT1(mims.size() == 1,
"Generic linear assembly brick needs one and only one "
"mesh_im");
- md.add_generic_expression(expr, *(mims[0]), region);
+ md.add_generic_expression(expr, *(mims[0]), region, secondary_domain);
}
virtual std::string declare_volume_assembly_string
@@ -3581,9 +3631,12 @@ model_complex_plain_vector &
gen_nonlinear_assembly_brick(const std::string &expr_, const mesh_im &mim,
bool is_sym,
- bool is_coer, std::string brickname = "") {
+ bool is_coer,
+ std::string brickname,
+ const std::string &secdom) {
if (brickname.size() == 0) brickname = "Generic linear assembly brick";
expr = expr_;
+ secondary_domain = secdom;
is_lower_dim = mim.is_lower_dimensional();
set_flags(brickname, false /* is linear*/,
is_sym /* is symmetric */, is_coer /* is coercive */,
@@ -3592,12 +3645,14 @@ model_complex_plain_vector &
};
- size_type add_nonlinear_term
+ size_type add_nonlinear_term_
(model &md, const mesh_im &mim, const std::string &expr, size_type region,
- bool is_sym, bool is_coercive, std::string brickname) {
+ bool is_sym, bool is_coercive, const std::string &brickname,
+ const std::string &secondary_domain) {
ga_workspace workspace(md);
- size_type order = workspace.add_expression(expr, mim, region);
+ size_type order = workspace.add_expression(expr, mim, region, 2,
+ secondary_domain);
GMM_ASSERT1(order < 2, "Order two test functions (Test2) are not allowed"
" in assembly string for nonlinear terms");
model::varnamelist vl, vl_test1, vl_test2, ddl, dl;
@@ -3607,13 +3662,30 @@ model_complex_plain_vector &
else vl.push_back(ddl[i]);
if (order == 0) { is_coercive = is_sym = true; }
pbrick pbr = std::make_shared<gen_nonlinear_assembly_brick>
- (expr, mim, is_sym, is_coercive, brickname);
+ (expr, mim, is_sym, is_coercive, brickname, secondary_domain);
model::termlist tl; // No term
// tl.push_back(model::term_description(true, is_sym));
// TODO to be changed.
return md.add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
}
+
+ size_type add_nonlinear_term
+ (model &md, const mesh_im &mim, const std::string &expr, size_type region,
+ bool is_sym, bool is_coercive, const std::string &brickname) {
+ return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
+ brickname, "");
+ }
+
+ size_type add_nonlinear_twodomain_term
+ (model &md, const mesh_im &mim, const std::string &expr, size_type region,
+ const std::string &secondary_domain, bool is_sym, bool is_coercive,
+ const std::string &brickname) {
+ return add_nonlinear_term_(md, mim, expr, region, is_sym, is_coercive,
+ brickname, secondary_domain);
+ }
+
+
// ----------------------------------------------------------------------
//
// Generic elliptic brick