Source code for custEM.fem.frequency_domain_approaches
# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
import dolfin as df
import numpy as np
from custEM.fem import ApproachBaseFD
from custEM.fem import check_source
from custEM.fem import add_sigma_vals_to_list
from custEM.fem import check_sigma_vals
from custEM.misc import logger_print as lp
"""
In this file, the different EM modeling approaches "A_V_nodal",
"A_V_mixed", "F_U_nodal", "E_vector", and "H_vector" are implemented.
The formulations are derived by Rochlitz (2020) (dissertation) and the
references therein.
"""
[docs]
class A_V_nodal(ApproachBaseFD):
"""
FE implementation of nodal potential approach, see (Badea, 2001) or
(Puzyref, 2013).
"""
def __init__(self, FS, MP):
"""
Initialize FE class for *An_t* or *An_s* approach.
Required arguments
------------------
- FS, type class
FunctionSpaces instance
- MP, type class
ModelParameters instance
"""
# inherit default attributes and common methods from ApproachBaseFD
super(self.__class__, self).__init__()
# initialize trial and test functions and pipe FE parameters
(self.v_1r, self.v_2r, self.v_3r, self.v_4r,
self.v_1i, self.v_2i, self.v_3i, self.v_4i) = df.TestFunctions(FS.M)
(self.u_1r, self.u_2r, self.u_3r, self.u_4r, self.u_1i,
self.u_2i, self.u_3i, self.u_4i) = df.TrialFunctions(FS.M)
self.FS = FS
self.MP = MP
self.p = self.FS.p
# mute print of the following FE class attributes
self.mute = ['v_1r', 'v_1i', 'v_2r', 'v_2i', 'v_3r', 'v_3i', 'DOM',
'v_4r', 'v_4i', 'u_1r', 'u_1i', 'u_2r', 'u_2i',
'u_3r', 'u_3i', 'u_4r', 'u_4i', 'dx', 'FS', 'MP', 'mute']
[docs]
def build_var_form(self, **fem_kwargs):
"""
Print FE parameters and initialize variational formulations for
assembly. Valid keyword arguments are described in the
**ApproachBaseFD** and **PrimaryFields** classes.
"""
self.update_fem_parameters(fem_kwargs)
self.bc = 'ZD' # force zero Dirichlet conditions for A-V-nodal
check_source(self.s_type, self.tx, self.pf_EH_flag)
self.lhs_form()
if '_t' in self.MP.approach:
self.rhs_form_total()
self.FS.add_primary_fields(self, self.MP, calc=False, **fem_kwargs)
else:
self.FS.anom_flag = check_sigma_vals(self.MP)
self.FS.add_primary_fields(self, self.MP, **fem_kwargs)
self.rhs_form_secondary()
# Notice! using old way of implementing multiple domains with different
# sigma values for A-V nodal approach, which is slower. The alternative
# way based on DG-paramter functions can be implemented as well if
# required by users. For all other appraoches, it is already the
# standard.
[docs]
def lhs_form(self):
"""
Left-hand side formulation as the basis for the system matrix assembly.
"""
C = df.Constant((self.MP.mu * self.MP.omega))
c1 = add_sigma_vals_to_list(self.MP.sigma)
dx = self.FS.DOM.dx
self.L = ((1./C) * df.inner(df.grad(self.u_1r),
df.grad(self.v_1r)) * dx(0) -
df.inner(c1[0] * self.u_1i, self.v_1r) * dx(0) -
df.inner(self.u_4i.dx(0), c1[0] * self.v_1r) * dx(0) -
df.inner(c1[0] * self.u_1r, self.v_1i) * dx(0) -
(1./C) * df.inner(df.grad(self.u_1i),
df.grad(self.v_1i)) * dx(0) -
df.inner(self.u_4r.dx(0), c1[0] * self.v_1i) * dx(0) +
(1./C) * df.inner(df.grad(self.u_2r),
df.grad(self.v_2r)) * dx(0) -
df.inner(c1[0] * self.u_2i, self.v_2r) * dx(0) -
df.inner(self.u_4i.dx(1), c1[0] * self.v_2r) * dx(0) -
df.inner(c1[0] * self.u_2r, self.v_2i) * dx(0) -
(1./C) * df.inner(df.grad(self.u_2i),
df.grad(self.v_2i)) * dx(0) -
df.inner(self.u_4r.dx(1), c1[0] * self.v_2i) * dx(0) +
(1./C) * df.inner(df.grad(self.u_3r),
df.grad(self.v_3r)) * dx(0) -
df.inner(c1[0] * self.u_3i, self.v_3r) * dx(0) -
df.inner(self.u_4i.dx(2), c1[0] * self.v_3r) * dx(0) -
df.inner(c1[0] * self.u_3r, self.v_3i) * dx(0) -
(1./C) * df.inner(df.grad(self.u_3i),
df.grad(self.v_3i)) * dx(0) -
df.inner(self.u_4r.dx(2), c1[0] * self.v_3i) * dx(0) -
df.inner(c1[0] * self.u_1i, self.v_4r.dx(0)) * dx(0) -
df.inner(c1[0] * self.u_2i, self.v_4r.dx(1)) * dx(0) -
df.inner(c1[0] * self.u_3i, self.v_4r.dx(2)) * dx(0) -
df.inner(c1[0] * df.grad(self.u_4i),
df.grad(self.v_4r)) * dx(0) -
df.inner(c1[0] * self.u_1r, self.v_4i.dx(0)) * dx(0) -
df.inner(c1[0] * self.u_2r, self.v_4i.dx(1)) * dx(0) -
df.inner(c1[0] * self.u_3r, self.v_4i.dx(2)) * dx(0) -
df.inner(c1[0] * df.grad(self.u_4r),
df.grad(self.v_4i)) * dx(0))
for j in range(1, self.FS.DOM.n_domains):
if type(self.MP.sigma[j]) is list and \
type(self.MP.sigma[j][0]) is float:
self.L += ((1./C) * df.inner(df.grad(self.u_1r),
df.grad(self.v_1r)) * dx(j) -
df.inner(c1[j] * self.u_1i, self.v_1r) * dx(j) -
df.inner(self.u_4i.dx(0), c1[j] * self.v_1r)*dx(j) -
df.inner(c1[j] * self.u_1r, self.v_1i) * dx(j) -
(1./C) * df.inner(df.grad(self.u_1i),
df.grad(self.v_1i)) * dx(j) -
df.inner(self.u_4r.dx(0), c1[j] * self.v_1i)*dx(j) +
(1./C) * df.inner(df.grad(self.u_2r),
df.grad(self.v_2r)) * dx(j) -
df.inner(c1[j] * self.u_2i, self.v_2r) * dx(j) -
df.inner(self.u_4i.dx(1), c1[j] * self.v_2r)*dx(j) -
df.inner(c1[j] * self.u_2r, self.v_2i) * dx(j) -
(1./C) * df.inner(df.grad(self.u_2i),
df.grad(self.v_2i)) * dx(j) -
df.inner(self.u_4r.dx(1), c1[j] * self.v_2i)*dx(j) +
(1./C) * df.inner(df.grad(self.u_3r),
df.grad(self.v_3r)) * dx(j) -
df.inner(c1[j] * self.u_3i, self.v_3r) * dx(j) -
df.inner(self.u_4i.dx(2), c1[j] * self.v_3r)*dx(j) -
df.inner(c1[j] * self.u_3r, self.v_3i) * dx(j) -
(1./C) * df.inner(df.grad(self.u_3i),
df.grad(self.v_3i)) * dx(j) -
df.inner(self.u_4r.dx(2), c1[j] * self.v_3i)*dx(j) -
df.inner(c1[j] * self.u_1i, self.v_4r.dx(0))*dx(j) -
df.inner(c1[j] * self.u_2i, self.v_4r.dx(1))*dx(j) -
df.inner(c1[j] * self.u_3i, self.v_4r.dx(2))*dx(j) -
df.inner(c1[j] * df.grad(self.u_4i),
df.grad(self.v_4r)) * dx(j) -
df.inner(c1[j] * self.u_1r, self.v_4i.dx(0))*dx(j) -
df.inner(c1[j] * self.u_2r, self.v_4i.dx(1))*dx(j) -
df.inner(c1[j] * self.u_3r, self.v_4i.dx(2))*dx(j) -
df.inner(c1[j] * df.grad(self.u_4r),
df.grad(self.v_4i)) * dx(j))
else:
self.L += (+(1.0/C) * df.inner(df.grad(self.u_1r),
df.grad(self.v_1r)) * dx(j) -
df.inner(c1[j][0][0] * self.u_1i, self.v_1r)*dx(j) -
df.inner(self.u_4i.dx(0),
c1[j][0][0] * self.v_1r) * dx(j) -
df.inner(c1[j][0][0] * self.u_1r, self.v_1i)*dx(j) -
(1.0/C) * df.inner(df.grad(self.u_1i),
df.grad(self.v_1i)) * dx(j) -
df.inner(self.u_4r.dx(0),
c1[j][0][0] * self.v_1i) * dx(j) +
(1.0/C) * df.inner(df.grad(self.u_2r),
df.grad(self.v_2r)) * dx(j) -
df.inner(c1[j][1][1] * self.u_2i, self.v_2r)*dx(j) -
df.inner(self.u_4i.dx(1),
c1[j][1][1] * self.v_2r) * dx(j) -
df.inner(c1[j][1][1] * self.u_2r, self.v_2i)*dx(j) -
(1.0/C) * df.inner(df.grad(self.u_2i),
df.grad(self.v_2i)) * dx(j) -
df.inner(self.u_4r.dx(1),
c1[j][1][1] * self.v_2i) * dx(j) +
(1.0/C) * df.inner(df.grad(self.u_3r),
df.grad(self.v_3r)) * dx(j) -
df.inner(c1[j][2][2] * self.u_3i, self.v_3r)*dx(j) -
df.inner(self.u_4i.dx(2),
c1[j][2][2] * self.v_3r) * dx(j) -
df.inner(c1[j][2][2] * self.u_3r, self.v_3i)*dx(j) -
(1.0/C) * df.inner(df.grad(self.u_3i),
df.grad(self.v_3i)) * dx(j) -
df.inner(self.u_4r.dx(2),
c1[j][2][2] * self.v_3i) * dx(j) -
df.inner(c1[j][0][0] * self.u_1i, self.v_4r.dx(0)) *
dx(j) - df.inner(c1[j][1][1] * self.u_2i,
self.v_4r.dx(1)) *
dx(j) - df.inner(c1[j][2][2] * self.u_3i,
self.v_4r.dx(2)) *
dx(j) - df.inner(c1[j] * df.grad(self.u_4i),
df.grad(self.v_4r)) * dx(j) -
df.inner(c1[j][0][0] * self.u_1r, self.v_4i.dx(0)) *
dx(j) - df.inner(c1[j][1][1] * self.u_2r,
self.v_4i.dx(1)) *
dx(j) - df.inner(c1[j][2][2] * self.u_3r,
self.v_4i.dx(2)) *
dx(j) - df.inner(c1[j] * df.grad(self.u_4r),
df.grad(self.v_4i)) * dx(j))
self.L = [self.L]
[docs]
def rhs_form_total(self):
"""
Right-hand side formulation for total field approach. Right-hand side
vector will be filled later by the **TotalFieldAssembler** instance.
"""
dx_0 = self.FS.DOM.dx_0
self.R = [(df.inner(df.Constant(("0.0")), self.v_1r) * dx_0 +
df.inner(df.Constant(("0.0")), self.v_1i) * dx_0 +
df.inner(df.Constant(("0.0")), self.v_2r) * dx_0 +
df.inner(df.Constant(("0.0")), self.v_2i) * dx_0 +
df.inner(df.Constant(("0.0")), self.v_3r) * dx_0 +
df.inner(df.Constant(("0.0")), self.v_3i) * dx_0 +
df.inner(df.Constant(("0.0")), self.v_4r) * dx_0 +
df.inner(df.Constant(("0.0")), self.v_4i) * dx_0)]
[docs]
def rhs_form_secondary(self, fi=0):
"""
Right-hand side formulation for secondary field approach.
"""
f_0 = df.Constant(('0'))
dx = self.FS.DOM.dx
PF = self.FS.PF
om1 = df.Constant(1./self.MP.omegas[fi])
c2 = add_sigma_vals_to_list(self.MP.delta_sigma)
self.R = []
for ti in range(self.n_tx):
R = (df.inner(f_0, self.v_1r) * dx(0) +
df.inner(f_0, self.v_1i) * dx(0) +
df.inner(f_0, self.v_2r) * dx(0) +
df.inner(f_0, self.v_2i) * dx(0) +
df.inner(f_0, self.v_3r) * dx(0) +
df.inner(f_0, self.v_3i) * dx(0) +
df.inner(f_0, self.v_4r) * dx(0) +
df.inner(f_0, self.v_4i) * dx(0))
if self.pf_EH_flag == 'H':
lp(self.MP.logger, 20,
'Notice! Using alternative primary magnetic fields.\n'
'For *A-Phi-nodal*, this requires transformation from '
'vector to scalar fields. Continuing ...', pre_dash=False)
import custEM as ce
Solver = ce.core.Solver(self.FS, self)
M = df.PETScMatrix()
w1, w2 = df.PETScVector(), df.PETScVector()
k_r = df.TrialFunction(self.FS.V)
k_i = df.TrialFunction(self.FS.V)
l_r = df.TestFunction(self.FS.V)
l_i = df.TestFunction(self.FS.V)
br = df.inner(k_r, l_r) * dx(0)
bi = df.inner(k_i, l_i) * dx(0)
H_0_r = df.interpolate(PF.H_0_r_cg[ti], self.FS.V)
H_0_i = df.interpolate(PF.H_0_i_cg[ti], self.FS.V)
Kr = df.inner(df.inv(self.MP.sigma_air[0]) * df.curl(H_0_r),
l_r) * dx(0)
Ki = df.inner(df.inv(self.MP.sigma_air[0]) * df.curl(H_0_i),
l_i) * dx(0)
sig_gi = add_sigma_vals_to_list(self.MP.sigma_0, inverse=True)
for j in range(1, len(self.MP.sigma)):
br += df.inner(k_r, l_r) * dx(j)
bi += df.inner(k_i, l_i) * dx(j)
Kr += df.inner(sig_gi[j-1] * df.curl(H_0_r), l_r) * dx(j)
Ki += df.inner(sig_gi[j-1] * df.curl(H_0_i), l_i) * dx(j)
H1, b1 = df.assemble_system(br, Kr, A_tensor=M, b_tensor=w1)
H2, b2 = df.assemble_system(bi, Ki, A_tensor=M, b_tensor=w2)
E_0_r = Solver.solve_system_mumps(H1, [b1], self.FS.V,
sym=True)
E_0_r_cg = df.interpolate(E_0_r[ti], self.FS.V_cg)
E_0_i = Solver.solve_system_mumps(H2, [b2], self.FS.V,
sym=True, first_call=False)
E_0_i_cg = df.interpolate(E_0_i[ti], self.FS.V_cg)
# singular hack for increased PF accuracy
# df.File('E_r.xml') << E_0_r_cg
# df.File('E_i.xml') << E_0_i_cg
# raise SystemExit
# mesh = df.Mesh('meshes/_xml/example_3_mesh_p2.xml')
# FS = df.VectorFunctionSpace(mesh, 'CG', 1)
# Er = df.Function(self.FS.V_cg, 'E_r.xml')
# Ei = df.Function(self.FS.V_cg, 'E_i.xml')
# Er2 = df.interpolate(Er, FS)
# Ei2 = df.interpolate(Ei, FS)
# df.File('E_r2.xml') << Er2
# df.File('E_i2.xml') << Ei2
# V_cg = df.VectorFunctionSpace(self.FS.mesh, 'CG', 1)
# E_0_r_cg = df.Function(V_cg, 'E_r2.xml')
# E_0_i_cg = df.Function(V_cg, 'E_i2.xml')
elif self.pf_EH_flag == 'E':
E_0_r_cg = PF.E_0_r_cg[ti]
E_0_i_cg = PF.E_0_i_cg[ti]
for j in range(1, self.FS.DOM.n_domains):
if type(self.MP.delta_sigma[j]) is list and \
type(self.MP.delta_sigma[j][0]) is float and \
abs(self.MP.delta_sigma[j][0]) < 1e-8:
if j != 1:
lp(self.MP.logger, 20,
'Notice! "delta_sigma" < 1e-8 is assumed to be 0.',
pre_dash=False)
R += (df.inner(f_0, self.v_1r) * dx(j) +
df.inner(f_0, self.v_1i) * dx(j) +
df.inner(f_0, self.v_2r) * dx(j) +
df.inner(f_0, self.v_2i) * dx(j) +
df.inner(f_0, self.v_3r) * dx(j) +
df.inner(f_0, self.v_3i) * dx(j) +
df.inner(f_0, self.v_4r) * dx(j) +
df.inner(f_0, self.v_4i) * dx(j))
else:
if self.pf_EH_flag in ['E', 'H']:
if type(self.MP.delta_sigma[j]) is list and \
type(self.MP.delta_sigma[j][0]) is float:
R += (-df.inner(om1 * c2[j] * E_0_r_cg.sub(0),
self.v_1r) * dx(j) +
df.inner(om1 * c2[j] * E_0_i_cg.sub(0),
self.v_1i) * dx(j) -
df.inner(om1 * c2[j] * E_0_r_cg.sub(1),
self.v_2r) * dx(j) +
df.inner(om1 * c2[j] * E_0_i_cg.sub(1),
self.v_2i) * dx(j) -
df.inner(om1 * c2[j] * E_0_r_cg.sub(2),
self.v_3r) * dx(j) +
df.inner(om1 * c2[j] * E_0_i_cg.sub(2),
self.v_3i) * dx(j) -
df.inner(om1 * c2[j] * E_0_r_cg,
df.grad(self.v_4r)) * dx(j) +
df.inner(om1 * c2[j] * E_0_i_cg,
df.grad(self.v_4i)) * dx(j))
else:
# manual implementation of anisotropic delta sigma
c2 = add_sigma_vals_to_list(self.MP.sigma)
c1 = add_sigma_vals_to_list(self.MP.sigma_0)
R += (-df.inner(c2[j][0][0] * E_0_r_cg.sub(0),
om1 * self.v_1r) * dx(j) +
df.inner(c2[j][0][0] * E_0_i_cg.sub(0),
om1 * self.v_1i) * dx(j) -
df.inner(c2[j][1][1] * E_0_r_cg.sub(1),
om1 * self.v_2r) * dx(j) +
df.inner(c2[j][1][1] * E_0_i_cg.sub(1),
om1 * self.v_2i) * dx(j) -
df.inner(c2[j][2][2] * E_0_r_cg.sub(2),
om1 * self.v_3r) * dx(j) +
df.inner(c2[j][2][2] * E_0_i_cg.sub(2),
om1 * self.v_3i) * dx(j) -
df.inner(om1 * c2[j] * E_0_r_cg,
df.grad(self.v_4r)) * dx(j) +
df.inner(om1 * c2[j] * E_0_i_cg,
df.grad(self.v_4i)) * dx(j))
R += -(-df.inner(c1[j-1][0][0] * E_0_r_cg.sub(0),
om1 * self.v_1r) * dx(j) +
df.inner(c1[j-1][0][0] * E_0_i_cg.sub(0),
om1 * self.v_1i) * dx(j) -
df.inner(c1[j-1][1][1] * E_0_r_cg.sub(1),
om1 * self.v_2r) * dx(j) +
df.inner(c1[j-1][1][1] * E_0_i_cg.sub(1),
om1 * self.v_2i) * dx(j) -
df.inner(c1[j-1][2][2] * E_0_r_cg.sub(2),
om1 * self.v_3r) * dx(j) +
df.inner(c1[j-1][2][2] * E_0_i_cg.sub(2),
om1 * self.v_3i) * dx(j) -
df.inner(om1 * c1[j-1] * E_0_r_cg,
df.grad(self.v_4r)) * dx(j) +
df.inner(om1 * c1[j-1] * E_0_i_cg,
df.grad(self.v_4i)) * dx(j))
self.R.append(R)
[docs]
class A_V_mixed(ApproachBaseFD):
"""
FE implementation of mixed potential approach, see (Ansari, 2014) or
a modified formulation by (Schwarzbach, 2009).
"""
def __init__(self, FS, MP):
"""
Initialize FE class for *Am_t* or *Am_s* approach.
Required arguments
------------------
- FS, type class
FunctionSpaces instance
- MP, type class
ModelParameters instance
"""
# inherit default attributes and common methods from ApproachBaseFD
super(self.__class__, self).__init__()
# initialize trial and test functions and pipe FE parameters
(self.v_1r, self.v_1i, self.v_2r, self.v_2i) = df.TestFunctions(FS.M)
(self.u_1r, self.u_1i, self.u_2r, self.u_2i) = df.TrialFunctions(FS.M)
self.FS = FS
self.MP = MP
self.p = self.FS.p
# mute print of the following FE class attributes
self.mute = ['v_1r', 'v_1i', 'v_2r', 'v_2i', 'u_1r', 'u_1i', 'DOM',
'u_2r', 'u_2i', 'dx', 'FS', 'MP', 'mute']
[docs]
def build_var_form(self, **fem_kwargs):
"""
Print FE parameters and initialize variational formulations for
assembly. Valid keyword arguments are described in the
**ApproachBaseFD** and **PrimaryFields** classes.
"""
self.update_fem_parameters(fem_kwargs)
check_source(self.s_type, self.tx, self.pf_EH_flag)
self.lhs_forms()
if '_t' in self.MP.approach:
self.FS.add_primary_fields(self, self.MP, calc=False, **fem_kwargs)
self.rhs_form_total()
else:
self.FS.anom_flag = check_sigma_vals(self.MP)
self.FS.add_primary_fields(self, self.MP, **fem_kwargs)
self.rhs_form_secondary()
[docs]
def lhs_forms(self):
"""
Left-hand side formulation as the basis for the system matrix assembly.
"""
self.L = []
for fi, omega in enumerate(self.MP.omegas):
om1 = df.Constant(1. / omega)
self.L.append(
om1 * df.inner(df.curl(self.u_1r), self.MP.mu_inv_func *
df.curl(self.v_1r)) * self.FS.DOM.dx -
df.inner(self.u_1i,
self.MP.sigma_func * self.v_1r) * self.FS.DOM.dx -
df.inner(df.grad(self.u_2i),
self.MP.sigma_func * self.v_1r) * self.FS.DOM.dx -
df.inner(self.u_1r,
self.MP.sigma_func * self.v_1i) * self.FS.DOM.dx -
om1 * df.inner(df.curl(self.u_1i), self.MP.mu_inv_func *
df.curl(self.v_1i)) * self.FS.DOM.dx -
df.inner(df.grad(self.u_2r),
self.MP.sigma_func * self.v_1i) * self.FS.DOM.dx -
df.inner(self.u_1i, self.MP.sigma_func *
df.grad(self.v_2r)) * self.FS.DOM.dx +
df.inner(df.grad(self.u_2i), self.MP.sigma_func *
df.grad(self.v_2r)) * self.FS.DOM.dx -
df.inner(self.u_1r, self.MP.sigma_func *
df.grad(self.v_2i)) * self.FS.DOM.dx +
df.inner(df.grad(self.u_2r), self.MP.sigma_func *
df.grad(self.v_2i)) * self.FS.DOM.dx)
[docs]
def rhs_form_total(self):
"""
Right-hand side formulation for total field approach. Right-hand side
vector will be filled later by the **TotalFieldAssembler** instance.
"""
self.R = [(df.inner(df.Constant(("0.0", "0.0", "0.0")),
self.v_1r) * self.FS.DOM.dx_0 +
df.inner(df.Constant(("0.0", "0.0", "0.0")),
self.v_1i) * self.FS.DOM.dx_0 +
df.inner(df.Constant(("0.0")),
self.v_2i) * self.FS.DOM.dx_0 +
df.inner(df.Constant(("0.0")),
self.v_2r) * self.FS.DOM.dx_0)]
[docs]
def rhs_form_secondary(self, fi=0):
"""
Right-hand side formulation for secondary potential field approach.
"""
om1 = df.Constant( 1./ self.MP.omegas[fi])
self.R = []
if self.pf_EH_flag == 'E':
for ti in range(self.n_tx):
self.R.append(
-om1 * df.inner(
self.MP.dsigma_func * self.FS.PF.E_0_r_cg[ti],
self.v_1r) * self.FS.DOM.dx +
om1 * df.inner(
self.MP.dsigma_func * self.FS.PF.E_0_i_cg[ti],
self.v_1i) * self.FS.DOM.dx -
om1 * df.inner(
self.MP.dsigma_func * self.FS.PF.E_0_r_cg[ti],
df.grad(self.v_2r)) * self.FS.DOM.dx +
om1 * df.inner(
self.MP.dsigma_func * self.FS.PF.E_0_i_cg[ti],
df.grad(self.v_2i)) * self.FS.DOM.dx)
elif self.pf_EH_flag == 'H':
for ti in range(self.n_tx):
self.R.append(
-om1 * df.inner(
self.MP.dsigma_func * self.MP.sigma0_inv_func *
df.curl(self.FS.PF.H_0_r_cg[ti]),
self.v_1r) * self.FS.DOM.dx +
om1 * df.inner(
self.MP.dsigma_func * self.MP.sigma0_inv_func *
df.curl(self.FS.PF.H_0_i_cg[ti]),
self.v_1i) * self.FS.DOM.dx -
om1 * df.inner(
self.MP.dsigma_func * self.MP.sigma0_inv_func *
df.curl(self.FS.PF.H_0_r_cg[ti]),
df.grad(self.v_2r)) * self.FS.DOM.dx +
om1 * df.inner(
self.MP.dsigma_func * self.MP.sigma0_inv_func *
df.curl(self.FS.PF.H_0_i_cg[ti]),
df.grad(self.v_2i)) * self.FS.DOM.dx)
[docs]
class F_U_mixed(ApproachBaseFD):
"""
FE implementation of Tau-Omega approach on mixed elements
(Mitsuhata, 2004), called F-U potential approach in custEM.
"""
def __init__(self, FS, MP):
"""
Initialize FE class for *Fm_s* approach.
Required arguments
------------------
- FS, type class
FunctionSpaces instance
- MP, type class
ModelParameters instance
"""
# inherit default attributes and common methods from ApproachBaseFD
super(self.__class__, self).__init__()
# initialize trial and test functions and pipe FE parameters
(self.v_1r, self.v_1i, self.v_2r, self.v_2i) = df.TestFunctions(FS.M)
(self.u_1r, self.u_1i, self.u_2r, self.u_2i) = df.TrialFunctions(FS.M)
self.FS = FS
self.MP = MP
self.p = self.FS.p
# mute print of the following FE class attributes
self.mute = ['v_1r', 'v_1i', 'v_2r', 'v_2i', 'u_1r', 'u_1i', 'DOM',
'u_2r', 'u_2i', 'dx', 'FS', 'MP', 'mute']
[docs]
def build_var_form(self, **fem_kwargs):
"""
Print FE parameters and initialize variational formulations for
assembly. Valid keyword arguments are described in the
**ApproachBaseFD** and **PrimaryFields** classes.
"""
self.update_fem_parameters(fem_kwargs)
check_source(self.s_type, self.tx, self.pf_EH_flag)
self.lhs_forms()
if '_t' in self.MP.approach:
self.FS.add_primary_fields(self, self.MP, calc=False, **fem_kwargs)
self.rhs_form_total()
else:
self.FS.anom_flag = check_sigma_vals(self.MP)
self.FS.add_primary_fields(self, self.MP, **fem_kwargs)
self.rhs_form_secondary()
[docs]
def lhs_forms(self):
"""
Left-hand side formulation as the basis for the system matrix assembly.
"""
dx = self.FS.DOM.dx
c1inv = add_sigma_vals_to_list(self.MP.sigma, inverse=True)
mu = df.Constant(self.MP.mu)
om1 = df.Constant(1./self.MP.omega)
self.L = []
for fi, omega in enumerate(self.MP.omegas):
om1 = df.Constant(1. / omega)
self.L.append(
om1 * df.inner(df.curl(self.u_1r), self.MP.sigma_inv_func *
df.curl(self.v_1r)) * self.FS.DOM.dx -
self.MP.mu_func * df.inner(self.u_1i,
self.v_1r) * self.FS.DOM.dx -
self.MP.mu_func * df.inner(df.grad(self.u_2i),
self.v_1r) * self.FS.DOM.dx -
self.MP.mu_func * df.inner(self.u_1r,
self.v_1i) * self.FS.DOM.dx -
om1 * df.inner(df.curl(self.u_1i), self.MP.sigma_inv_func *
df.curl(self.v_1i)) * self.FS.DOM.dx -
self.MP.mu_func * df.inner(df.grad(self.u_2r),
self.v_1i) * self.FS.DOM.dx -
df.inner(self.u_1i, self.MP.mu_func *
df.grad(self.v_2r)) * self.FS.DOM.dx +
df.inner(df.grad(self.u_2i), self.MP.mu_func *
df.grad(self.v_2r)) * self.FS.DOM.dx -
df.inner(self.u_1r, self.MP.mu_func *
df.grad(self.v_2i)) * self.FS.DOM.dx +
df.inner(df.grad(self.u_2r), self.MP.mu_func *
df.grad(self.v_2i)) * self.FS.DOM.dx)
[docs]
def rhs_form_total(self):
"""
Right-hand side formulation for total field approach. Right-hand side
vector will be filled later by the **TotalFieldAssembler** instance.
"""
self.R = (df.inner(df.Constant(("0.0", "0.0", "0.0")),
self.v_1r) * self.FS.DOM.dx +
df.inner(df.Constant(("0.0", "0.0", "0.0")),
self.v_1i) * self.FS.DOM.dx +
df.inner(df.Constant(("0.0")), self.v_2i) * self.FS.DOM.dx +
df.inner(df.Constant(("0.0")), self.v_2r) * self.FS.DOM.dx)
[docs]
def rhs_form_secondary(self, fi=0):
"""
Right-hand side formulation for secondary potential field approach.
"""
om1 = df.Constant(1. / self.MP.omegas[fi])
self.R = []
if self.pf_EH_flag == 'E':
for ti in range(self.n_tx):
self.R.append(
om1 * df.inner(
self.MP.sigma_inv_func * self.MP.dsigma_func *
self.FS.PF.E_0_r_cg[ti],
df.curl(self.v_1r)) * self.FS.DOM.dx -
om1 * df.inner(
self.MP.sigma_inv_func * self.MP.dsigma_func *
self.FS.PF.E_0_i_cg[ti],
df.curl(self.v_1i)) * self.FS.DOM.dx +
df.inner(df.Constant((0.)), self.v_2r) * self.FS.DOM.dx +
df.inner(df.Constant((0.)), self.v_2i) * self.FS.DOM.dx)
elif self.pf_EH_flag == 'H':
for ti in range(self.n_tx):
self.R.append(
om1 * df.inner(
(self.MP.sigma0_inv_func - self.MP.sigma_inv_func) *
df.curl(self.FS.PF.H_0_r_cg[ti]),
df.curl(self.v_1r)) * self.FS.DOM.dx -
om1 * df.inner(
(self.MP.sigma0_inv_func - self.MP.sigma_inv_func) *
df.curl(self.FS.PF.H_0_i_cg[ti]),
df.curl(self.v_1i)) * self.FS.DOM.dx +
df.inner(df.Constant((0.)), self.v_2r) * self.FS.DOM.dx +
df.inner(df.Constant((0.)), self.v_2i) * self.FS.DOM.dx)
[docs]
class F_U_nodal(ApproachBaseFD):
"""
FE implementation of F-U (Tau-Omega) approach on nodal elements. This
approach was not implemented yet as it requires a homogeneous conductivity
distribution (at least in our derivation of the equations) in the complete
computational domain. Hence, it is irrelevant for geo-EM applications.
"""
def __init__(self, FS, MP):
"""
Initialize FE class for *Fn_s* approach.
Required arguments
------------------
- FS, type class
FunctionSpaces instance
- MP, type class
ModelParameters instance
"""
# dummy, not implemented because unsuited for geo-EM modeling
pass
[docs]
class E_vector(ApproachBaseFD):
"""
FE implementation of E-field approach (Grayver, 2014, Schwarzbach 2009).
"""
def __init__(self, FS, MP):
"""
Initialize FE class for *E_t* or *E_s* approach.
Required arguments
------------------
- FS, type class
FunctionSpaces instance
- MP, type class
ModelParameters instance
"""
# inherit default attributes and common methods from ApproachBaseFD
super(self.__class__, self).__init__()
# initialize trial and test functions and pipe FE parameters
(self.v_r, self.v_i) = df.TestFunctions(FS.M)
(self.u_r, self.u_i) = df.TrialFunctions(FS.M)
self.vh=df.TestFunction(FS.M1)
self.uh=df.TrialFunction(FS.M1)
self.FS = FS
self.MP = MP
self.p = self.FS.p
# mute print of the following FE class attributes
self.mute = ['v_r', 'v_i', 'u_r', 'u_i','vh','uh', 'DOM',
'dx', 'FS', 'MP', 'mute']
[docs]
def build_var_form(self,**fem_kwargs):
"""
Print FE parameters and initialize variational formulations for
assembly. Valid keyword arguments are described in the
**ApproachBaseFD** and **PrimaryFields** classes.
"""
self.update_fem_parameters(fem_kwargs)
check_source(self.s_type, self.tx, self.pf_EH_flag,
approach=self.MP.approach)
if self.quasi_static and not self.ip:
self.lhs_forms(self.MP.system_it)
else:
self.lhs_forms_full()
if '_t' in self.MP.approach or self.MP.approach == 'MT':
self.FS.add_primary_fields(self, self.MP, calc=False, **fem_kwargs)
self.rhs_form_total(self.MP.system_it)
else:
self.FS.anom_flag = check_sigma_vals(self.MP)
self.FS.add_primary_fields(self, self.MP, **fem_kwargs)
self.rhs_form_secondary()
[docs]
def lhs_forms(self,system_it):
"""
Left-hand side formulation as the basis for the system matrix assembly.
"""
self.L = []
if self.MP.system_it==False:
for fi, omega in enumerate(self.MP.omegas):
om = df.Constant(omega)
self.L.append(
df.inner(self.MP.mu_inv_func * df.curl(self.u_r),
df.curl(self.v_r)) * self.FS.DOM.dx_0 -
om * df.inner(self.MP.sigma_func * self.u_i,
self.v_r) * self.FS.DOM.dx_0 -
om * df.inner(self.MP.sigma_func * self.u_r,
self.v_i) * self.FS.DOM.dx_0 -
df.inner(self.MP.mu_inv_func * df.curl(self.u_i),
df.curl(self.v_i)) * self.FS.DOM.dx_0)
else:
self.H = []
self.B = []
self.A = []
for fi, omega in enumerate(self.MP.omegas):
om = df.Constant(omega)
self.L.append(
om * df.inner(self.MP.sigma_func * self.u_r,
self.v_r) * self.FS.DOM.dx_0 +
df.inner(self.MP.mu_inv_func * df.curl(self.u_i),
df.curl(self.v_r)) * self.FS.DOM.dx_0 -
df.inner(self.MP.mu_inv_func * df.curl(self.u_r),
df.curl(self.v_i)) * self.FS.DOM.dx_0 +
om * df.inner(self.MP.sigma_func * self.u_i,
self.v_i) * self.FS.DOM.dx_0)
self.B.append(df.inner(self.MP.mu_inv_func * df.curl(self.uh),
df.curl(self.vh)) * self.FS.DOM.dx_0)
self.H.append(
df.inner(self.MP.mu_inv_func * df.curl(self.uh),
df.curl(self.vh)) * self.FS.DOM.dx_0 +
om * df.inner(self.MP.sigma_func * self.uh,
self.vh) * self.FS.DOM.dx_0)
self.A.append(
om * df.inner(self.MP.sigma_func * self.uh,
self.vh) * self.FS.DOM.dx_0)
[docs]
def lhs_forms_full(self):
"""
Left-hand side formulation as the basis for the system matrix assembly.
Here, electric permittivities and induced polarization parameters are
considered in addition to the conductivity values.
"""
self.L = []
if self.MP.system_it==False:
for fi, omega in enumerate(self.MP.omegas):
om = df.Constant(omega)
if self.ip:
if fi == 0:
lp(self.MP.logger, 20,
' - modeling with induced polarization effects '
'enabled - ', post_dash=False)
if self.MP.ip_m == 0.:
lp(self.MP.logger, 50,
'Error! Specify at least the model parameter '
'**ip_m** if induced polarization effects '
'should be taken into account by setting '
'**ip=True**\n.'
' - default value for ip_c = 1. - \n'
' - default value for ip_tau = 1. (s) - \n'
'Aborting ...')
raise SystemExit
sigma_real_func, sigma_imag_func = \
self.FS.build_complex_sigma_func(
omega, self.MP.sigma,
self.MP.ip_tau, self.MP.ip_m, self.MP.ip_c)
else:
if fi == 0:
sigma_real_func = self.MP.sigma_func
sigma_imag_func = df.Constant((0.))
if fi == 0:
if self.quasi_static:
self.eps_func = df.Constant((0.))
else:
self.eps_func = self.MP.eps_func
self.L.append(
df.inner(self.MP.mu_inv_func * df.curl(self.u_r),
df.curl(self.v_r)) * self.FS.DOM.dx_0 -
(om * om * self.eps_func + om * sigma_imag_func) *
df.inner(self.u_r, self.v_r) * self.FS.DOM.dx_0 -
om * df.inner(sigma_real_func * self.u_i,
self.v_r) * self.FS.DOM.dx_0 -
om * df.inner(sigma_real_func * self.u_r,
self.v_i) * self.FS.DOM.dx_0 -
df.inner(self.MP.mu_inv_func * df.curl(self.u_i),
df.curl(self.v_i)) * self.FS.DOM.dx_0 +
(om * om * self.eps_func + om * sigma_imag_func) *
df.inner(self.u_i, self.v_i) * self.FS.DOM.dx_0)
else:
self.H = []
self.B = []
self.A = []
for fi, omega in enumerate(self.MP.omegas):
om = df.Constant(omega)
if self.ip:
if fi == 0:
lp(self.MP.logger, 20,
' - modeling with induced polarization effects '
'enabled - ', post_dash=False)
if self.MP.ip_m == 0.:
lp(self.MP.logger, 50,
'Error! Specify at least the model parameter '
'**ip_m** if induced polarization effects '
'should be taken into account by setting '
'**ip=True**\n.'
' - default value for ip_c = 1. - \n'
' - default value for ip_tau = 1. (s) - \n'
'Aborting ...')
raise SystemExit
sigma_real_func, sigma_imag_func = \
self.FS.build_complex_sigma_func(
omega, self.MP.sigma,
self.MP.ip_tau, self.MP.ip_m, self.MP.ip_c)
else:
if fi == 0:
sigma_real_func = self.MP.sigma_func
sigma_imag_func = df.Constant((0.))
if fi == 0:
if self.quasi_static:
self.eps_func = df.Constant((0.))
else:
self.eps_func = self.MP.eps_func
self.L.append(
om * df.inner(sigma_real_func * self.u_r,
self.v_r) * self.FS.DOM.dx_0 +
df.inner(self.MP.mu_inv_func * df.curl(self.u_i),
df.curl(self.v_r)) * self.FS.DOM.dx_0 -
(om * om * self.eps_func + om * sigma_imag_func) *
df.inner(self.u_i, self.v_r) * self.FS.DOM.dx_0 -
df.inner(self.MP.mu_inv_func * df.curl(self.u_r),
df.curl(self.v_i)) * self.FS.DOM.dx_0 +
(om * om * self.eps_func + om * sigma_imag_func) *
df.inner(self.u_r, self.v_i) * self.FS.DOM.dx_0 +
om * df.inner(sigma_real_func * self.u_i,
self.v_i) * self.FS.DOM.dx_0)
self.B.append(
df.inner(self.MP.mu_inv_func * df.curl(self.uh),
df.curl(self.vh)) * self.FS.DOM.dx_0 -
(om * om * self.eps_func + om * sigma_imag_func) *
df.inner(self.uh, self.vh) * self.FS.DOM.dx_0)
self.H.append(
df.inner(self.MP.mu_inv_func * df.curl(self.uh),
df.curl(self.vh)) * self.FS.DOM.dx_0 -
(om * om * self.eps_func + om * sigma_imag_func) *
df.inner(self.uh, self.vh) * self.FS.DOM.dx_0 +
om * df.inner(sigma_real_func * self.uh,
self.vh) * self.FS.DOM.dx_0)
self.A.append(
om * df.inner(sigma_real_func * self.uh,
self.vh) * self.FS.DOM.dx_0)
# old code
# eps = [df.Constant(elem) for elem in self.MP.eps]
# c1 = add_sigma_vals_to_list(self.MP.sigma)
# L = ((1./mu) * df.inner(df.curl(self.u_r),
# df.curl(self.v_r)) * dx(0) -
# (eps[0] * om * om) * df.inner(self.u_r,
# self.v_r) * dx(0) -
# om * df.inner(c1[0] * self.u_i, self.v_r) * dx(0) -
# om * df.inner(c1[0] * self.u_r, self.v_i) * dx(0) -
# (1./mu) * df.inner(df.curl(self.u_i),
# df.curl(self.v_i)) * dx(0) +
# (eps[0] * om * om) * df.inner(self.u_i,
# self.v_i) * dx(0))
# for j in range(1, self.FS.DOM.n_domains):
# L += (
# (1./mu) * df.inner(
# df.curl(self.u_r), df.curl(self.v_r)) * dx(j) -
# (eps[j] * om * om) * df.inner(self.u_r,
# self.v_r) * dx(j) -
# om * df.inner(c1[j] * self.u_i, self.v_r) * dx(j) -
# om * df.inner(c1[j] * self.u_r, self.v_i) * dx(j) -
# (1./mu) * df.inner(
# df.curl(self.u_i), df.curl(self.v_i)) * dx(j) +
# (eps[j] * om * om) * df.inner(self.u_i,
# self.v_i) * dx(j))
# self.L.append(L)
[docs]
def rhs_form_total(self,system_it):
"""
Right-hand side formulation for total field approach. Right-hand side
vector will be filled later by the **TotalFieldAssembler** instance.
"""
dx_0 = self.FS.DOM.dx_0
if self.MP.system_it==False:
self.R = [(df.inner(df.Constant(("0.0", "0.0", "0.0")),
self.v_r) * self.FS.DOM.dx_0 +
df.inner(df.Constant(("0.0", "0.0", "0.0")),
self.v_i) * self.FS.DOM.dx_0)]
else:
self.R = [(df.inner(df.Constant(("0.0", "0.0", "0.0")),
self.v_r) * self.FS.DOM.dx_0 +
df.inner(df.Constant(("0.0", "0.0", "0.0")),
self.v_i) * self.FS.DOM.dx_0)]
[docs]
def rhs_form_secondary(self, fi=0):
"""
Right-hand side formulation for secondary electric field approach.
"""
om = df.Constant(self.MP.omegas[fi])
PF = self.FS.PF
self.R = []
if self.pf_EH_flag == 'E':
for ti in range(self.n_tx):
if self.FS.PF.nedelec_pf:
self.R.append(
om * df.inner(self.MP.dsigma_func * PF.E_0_i[ti],
self.v_r) * self.FS.DOM.dx_0 +
om * df.inner(self.MP.dsigma_func * PF.E_0_r[ti],
self.v_i) * self.FS.DOM.dx_0)
else:
self.R.append(
om * df.inner(self.MP.dsigma_func * PF.E_0_i_cg[ti],
self.v_r) * self.FS.DOM.dx_0 +
om * df.inner(self.MP.dsigma_func * PF.E_0_r_cg[ti],
self.v_i) * self.FS.DOM.dx_0)
elif self.pf_EH_flag == 'H':
for ti in range(self.n_tx):
if self.FS.PF.nedelec_pf:
self.R.append(
om * df.inner(
self.MP.dsigma_func * self.MP.sigma0_inv_func *
df.curl(PF.H_0_i[ti]),
self.v_r) * self.FS.DOM.dx_0 +
om * df.inner(
self.MP.dsigma_func * self.MP.sigma0_inv_func *
df.curl(PF.H_0_r[ti]), self.v_i) *
self.FS.DOM.dx_0)
else:
self.R.append(
om * df.inner(
self.MP.dsigma_func * self.MP.sigma0_inv_func *
df.curl(PF.H_0_i_cg[ti]),
self.v_r) * self.FS.DOM.dx_0 +
om * df.inner(
self.MP.dsigma_func * self.MP.sigma0_inv_func *
df.curl(PF.H_0_r_cg[ti]), self.v_i) *
self.FS.DOM.dx_0)
[docs]
class H_vector(ApproachBaseFD):
"""
FE implementation of H-fild approach (Rochlitz, 2019).
"""
def __init__(self, FS, MP):
"""
Initialize FE class for *An_t* or *An_s* approach.
Required arguments
------------------
- FS, type class
FunctionSpaces instance
- MP, type class
ModelParameters instance
"""
# inherit default attributes and common methods from ApproachBaseFD
super(self.__class__, self).__init__()
# initialize trial and test functions and pipe FE parameters
(self.v_r, self.v_i) = df.TestFunctions(FS.M)
(self.u_r, self.u_i) = df.TrialFunctions(FS.M)
self.vh=df.TestFunction(FS.M1)
self.uh=df.TrialFunction(FS.M1)
self.FS = FS
self.MP = MP
self.p = self.FS.p
# mute print of the following FE class attributes
self.mute = ['v_r', 'v_i', 'u_r', 'u_i', 'vh', 'uh', 'DOM',
'dx', 'FS', 'MP', 'mute']
[docs]
def build_var_form(self, add_primary=True, **fem_kwargs):
"""
Print FE parameters and initialize variational formulations for
assembly. Valid keyword arguments are described in the
**ApproachBaseFD** and **PrimaryFields** classes.
"""
self.update_fem_parameters(fem_kwargs)
check_source(self.s_type, self.tx, self.pf_EH_flag)
if self.quasi_static:
self.lhs_forms(self.MP.system_it)
else:
self.lhs_form_full()
if '_s' in self.MP.approach:
lp(self.MP.logger, 50,
'*full* formulation not implemented on RHS yet for '
'secondary magenticc approach. Aborting ...')
raise SystemExit
if '_t' in self.MP.approach:
lp(self.MP.logger, 50,
'Error! Total H-field approach not working yet!\n'
'Aborting ...')
raise SystemExit
self.rhs_form_total()
self.FS.add_primary_fields(self, self.MP, calc=False, **fem_kwargs)
else:
self.FS.anom_flag = check_sigma_vals(self.MP)
if add_primary:
self.FS.add_primary_fields(self, self.MP, **fem_kwargs)
self.rhs_form_secondary()
[docs]
def lhs_forms(self,system_it):
"""
Left-hand side formulation as the basis for the system matrix assembly.
"""
self.L = []
if self.MP.system_it==False:
for omega in self.MP.omegas:
om = df.Constant(omega)
self.L.append(
df.inner(self.MP.sigma_inv_func * df.curl(self.u_r),
df.curl(self.v_r)) * self.FS.DOM.dx_0 -
om * df.inner(self.MP.mu_func * self.u_i,
self.v_r) * self.FS.DOM.dx_0 -
om * df.inner(self.MP.mu_func * self.u_r,
self.v_i) * self.FS.DOM.dx_0 -
df.inner(self.MP.sigma_inv_func * df.curl(self.u_i),
df.curl(self.v_i)) * self.FS.DOM.dx_0)
else:
self.H = []
self.B = []
self.A = []
for omega in self.MP.omegas:
om = df.Constant(omega)
self.L.append(
om * df.inner(self.MP.mu_func * self.u_r,
self.v_r) * self.FS.DOM.dx_0 +
df.inner(self.MP.sigma_inv_func * df.curl(self.u_i),
df.curl(self.v_r)) * self.FS.DOM.dx_0 -
df.inner(self.MP.sigma_inv_func * df.curl(self.u_r),
df.curl(self.v_i)) * self.FS.DOM.dx_0 +
om * df.inner(self.MP.mu_func * self.u_i,
self.v_i) * self.FS.DOM.dx_0)
self.B.append(df.inner(self.MP.sigma_inv_func * df.curl(self.uh),
df.curl(self.vh)) * self.FS.DOM.dx_0)
self.H.append(
df.inner(self.MP.sigma_inv_func * df.curl(self.uh),
df.curl(self.vh)) * self.FS.DOM.dx_0 +
om * df.inner(self.MP.mu_func * self.uh,
self.vh) * self.FS.DOM.dx_0)
self.A.append(
om * df.inner(self.MP.mu_func * self.uh,
self.vh) * self.FS.DOM.dx_0)
[docs]
def lhs_form_full(self):
"""
Left-hand side contributions for full (not quasi-static) H-field
approach, in development, not fully tested yet.
"""
lp(self.MP.logger, 50,
'Error! Need to correct permittivity support for H field approach. '
'Aborting ...')
raise SystemExit
om = df.Constant(self.MP.omega)
mu = df.Constant(self.MP.mu)
eps = df.Constant(self.MP.eps)
s_func = build_dg_func(self.FS.S_dg,
self.FS.DOM.domain_func,
self.MP.sigma)
self.s_func = s_func
fact1 = (om * eps) / (om * om * eps * eps + s_func * s_func)
fact2 = s_func / (om * om * eps * eps + s_func * s_func)
self.fact1 = fact1
self.fact2 = fact2
lp(self.MP.logger, 30,
'Warning! Full magenetic field approach is experimental.\n'
'Only supported for single frequency modeling. Continuing ...',
post_dash=True)
dx_0 = self.FS.DOM.dx_0
inv_s_func = build_dg_func(self.FS.S_dg,
self.FS.DOM.domain_func,
self.MP.sigma, inverse=True)
self.inv_s_func = inv_s_func
self.L = (df.inner(df.curl(self.u_r),
df.curl(fact2 * self.v_r)) * dx_0 +
df.inner(df.curl(self.u_i), df.curl(fact1 * self.v_r))
* dx_0 -
mu * om * df.inner(self.u_i, self.v_r) * dx_0 +
df.inner(df.curl(self.u_r), df.curl(fact1 * self.v_i))
* dx_0 -
mu * om * df.inner(self.u_r, self.v_i) * dx_0 -
df.inner(df.curl(self.u_i),
df.curl(fact2 * self.v_i)) * dx_0)
[docs]
def rhs_form_total(self):
"""
Right-hand side formulation for total field approach. Right-hand side
vector will be filled later by the **TotalFieldAssembler** instance.
"""
# class Delta(df.UserExpression):
# def __init__(self, eps, x0, **kwargs):
# self.eps = eps
# self.x0 = x0
# df.UserExpression.__init__(self, **kwargs)
# def eval(self, values, x):
# eps = self.eps
# values[0] = 0.
# values[1] = 0.
# values[2] = eps/np.pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
# def value_shape(self): return (3, )
# delta = Delta(eps=1E-4, x0=np.array([0., 0., 1.]), degree=5)
dx_0 = self.FS.DOM.dx_0
self.R = [(df.inner(df.Constant(("0.0", "0.0", "0.0")),
self.v_r) * dx_0 +
df.inner(df.Constant(("0.0", "0.0", "0.0")),
self.v_i) * dx_0)]
[docs]
def rhs_form_secondary(self, fi=0):
"""
Right-hand side formulation for secondary magentic field approach.
"""
self.R = []
if self.pf_EH_flag == 'H':
for ti in range(self.n_tx):
if self.FS.PF.nedelec_pf:
self.R.append(
df.inner(
(self.MP.sigma0_inv_func - self.MP.sigma_inv_func) *
df.curl(self.FS.PF.H_0_r[ti]),
df.curl(self.v_r)) * self.FS.DOM.dx_0 -
df.inner(
(self.MP.sigma0_inv_func - self.MP.sigma_inv_func) *
df.curl(self.FS.PF.H_0_i[ti]),
df.curl(self.v_i)) * self.FS.DOM.dx_0)
else:
self.R.append(
df.inner(
(self.MP.sigma0_inv_func - self.MP.sigma_inv_func) *
df.curl(self.FS.PF.H_0_r_cg[ti]),
df.curl(self.v_r)) * self.FS.DOM.dx_0 -
df.inner(
(self.MP.sigma0_inv_func - self.MP.sigma_inv_func) *
df.curl(self.FS.PF.H_0_i_cg[ti]),
df.curl(self.v_i)) * self.FS.DOM.dx_0)
elif self.pf_EH_flag == 'E':
for ti in range(self.n_tx):
if self.FS.PF.nedelec_pf:
self.R.append(
df.inner(
(self.MP.sigma_inv_func * self.MP.dsigma_func) *
self.FS.PF.E_0_r[ti],
df.curl(self.v_r)) * self.FS.DOM.dx_0 -
df.inner(
(self.MP.sigma_inv_func * self.MP.dsigma_func) *
self.FS.PF.E_0_i[ti],
df.curl(self.v_i)) * self.FS.DOM.dx_0)
else:
self.R.append(
df.inner(
(self.MP.sigma_inv_func * self.MP.dsigma_func) *
self.FS.PF.E_0_r_cg[ti],
df.curl(self.v_r)) * self.FS.DOM.dx_0 -
df.inner(
(self.MP.sigma_inv_func * self.MP.dsigma_func) *
self.FS.PF.E_0_i_cg[ti],
df.curl(self.v_i)) * self.FS.DOM.dx_0)