# diffusion/laplace_coupling_lcbcs.pyΒΆ

Description

Two Laplace equations with multiple linear combination constraints.

The two equations are coupled by a periodic-like boundary condition constraint with a shift, given as a non-homogeneous linear combination boundary condition.

Find such that:

where , are given functions (shifts), is the periodic coordinate mapping and , and are unknown constant values - the unknown DOFs in , and are replaced by the integral mean values.

View the results using:

$./postproc.py square_quad.vtk -b --wireframe -d'u1,plot_warp_scalar,rel_scaling=1:u2,plot_warp_scalar,rel_scaling=1' source code r""" Two Laplace equations with multiple linear combination constraints. The two equations are coupled by a periodic-like boundary condition constraint with a shift, given as a non-homogeneous linear combination boundary condition. Find :math:u such that: .. math:: \int_{\Omega_1} \nabla v_1 \cdot \nabla u_1 = 0 \;, \quad \forall v_1 \;, \int_{\Omega_2} \nabla v_2 \cdot \nabla u_2 = 0 \;, \quad \forall v_2 \;, u_1 = 0 \mbox{ on } \Gamma_{bottom} \;, u_2 = 1 \mbox{ on } \Gamma_{top} \;, u_1(\ul{x}) = u_2(\ul{x}) + a(\ul{x}) \mbox{ for } \ul{x} \in \Gamma = \bar\Omega_1 \cap \bar\Omega_2 u_1(\ul{x}) = u_1(\ul{y}) + b(\ul{y}) \mbox{ for } \ul{x} \in \Gamma_{left}, \ul{y} \in \Gamma_{right}, \ul{y} = P(\ul{x}) \;, u_1 = c_{11} \mbox{ in } \Omega_{m11} \subset \Omega_1 \;, u_1 = c_{12} \mbox{ in } \Omega_{m12} \subset \Omega_1 \;, u_2 = c_2 \mbox{ in } \Omega_{m2} \subset \Omega_2 \;, where :math:a(\ul{x}), :math:b(\ul{y}) are given functions (shifts), :math:P is the periodic coordinate mapping and :math:c_{11}, :math:c_{12} and :math:c_2 are unknown constant values - the unknown DOFs in :math:\Omega_{m11}, :math:\Omega_{m12} and :math:\Omega_{m2} are replaced by the integral mean values. View the results using::$ ./postproc.py square_quad.vtk -b --wireframe -d'u1,plot_warp_scalar,rel_scaling=1:u2,plot_warp_scalar,rel_scaling=1'
"""
from __future__ import absolute_import
import numpy as nm

import sfepy.discrete.fem.periodic as per
from sfepy import data_dir

options = {
'nls' : 'newton',
'ls' : 'ls',
}

def get_shift1(ts, coors, region):
val = 0.1 * coors[:, 0]

return val

def get_shift2(ts, coors, region):
val = nm.empty_like(coors[:, 1])
val.fill(0.3)

return val

functions = {
'get_shift1' : (get_shift1,),
'get_shift2' : (get_shift2,),
'match_y_line' : (per.match_y_line,),
'match_x_line' : (per.match_x_line,),
}

fields = {
'scalar1': ('real', 1, 'Omega1', 1),
'scalar2': ('real', 1, 'Omega2', 1),
}

materials = {
}

variables = {
'u1' : ('unknown field', 'scalar1', 0),
'v1' : ('test field', 'scalar1', 'u1'),
'u2' : ('unknown field', 'scalar2', 1),
'v2' : ('test field', 'scalar2', 'u2'),
}

regions = {
'Omega1' : 'cells of group 1',
'Omega2' : 'cells of group 2',
'Omega_m1' : 'r.Omega1 -v (r.Gamma +s vertices of surface)',
'Omega_m11' : 'r.Omega_m1 *v vertices in (x < 0)',
'Omega_m12' : 'r.Omega_m1 *v vertices in (x > 0)',
'Omega_m2' : 'r.Omega2 -v (r.Gamma +s vertices of surface)',
'Left' : ('vertices in (x < -0.499)', 'facet'),
'Right' : ('vertices in (x > 0.499)', 'facet'),
'Bottom' : ('vertices in ((y < -0.499))', 'facet'),
'Top' : ('vertices in ((y > 0.499))', 'facet'),
'Gamma' : ('r.Omega1 *v r.Omega2 -v vertices of surface', 'facet'),
'Gamma1' : ('copy r.Gamma', 'facet', 'Omega1'),
'Gamma2' : ('copy r.Gamma', 'facet', 'Omega2'),
}

ebcs = {
'fix1' : ('Top', {'u2.all' : 1.0}),
'fix2' : ('Bottom', {'u1.all' : 0.0}),
}

lcbcs = {
'shifted1' : (('Gamma1', 'Gamma2'),
{'u1.all' : 'u2.all'},
'match_x_line', 'shifted_periodic',
'get_shift1'),
'shifted2' : (('Left', 'Right'),
{'u1.all' : 'u1.all'},
'match_y_line', 'shifted_periodic',
'get_shift2'),
'mean11' : ('Omega_m11', {'u1.all' : None}, None, 'integral_mean_value'),
'mean12' : ('Omega_m12', {'u1.all' : None}, None, 'integral_mean_value'),
'mean2' : ('Omega_m2', {'u2.all' : None}, None, 'integral_mean_value'),
}

integrals = {
'i1' : 2,
}

equations = {
'eq1' : """
dw_laplace.i1.Omega1(v1, u1) = 0
""",
'eq2' : """
dw_laplace.i1.Omega2(v2, u2) = 0
""",
}

solvers = {
'ls' : ('ls.scipy_direct', {}),
'newton' : ('nls.newton', {
'i_max'      : 1,
'eps_a'      : 1e-10,
}),
}