# Copyright 2023, 2024 Johannes Reding, Julius Zimmermann
# SPDX-License-Identifier: GPL-3.0-or-later
import logging
import ngsolve
from ossdbs.fem.solver import Solver
from ossdbs.fem.volume_conductor.volume_conductor_model import VolumeConductor
from ossdbs.model_geometry import ModelGeometry
from .conductivity import ConductivityCF
_logger = logging.getLogger(__name__)
[docs]class VolumeConductorNonFloating(VolumeConductor):
"""Model for representing a volume conductor which evaluates the potential."""
def __init__(
self,
geometry: ModelGeometry,
conductivity: ConductivityCF,
solver: Solver,
order: int,
meshing_parameters: dict,
output_path: str = "Results",
) -> None:
super().__init__(
geometry, conductivity, solver, order, meshing_parameters, output_path
)
# dirichlet boundaries (no surface impedance)
boundaries = [
contact.name
for contact in self.contacts.active
if contact.surface_impedance_model is None
]
self._space = self.h1_space(boundaries=boundaries, is_complex=self.is_complex)
self.update_space()
[docs] def update_space(self):
"""Update space (e.g., if mesh changes)."""
self._space.Update()
self._potential = ngsolve.GridFunction(space=self._space)
[docs] def compute_solution(self, frequency: float) -> ngsolve.comp.GridFunction:
"""Evaluate electrical potential of volume conductor.
Parameters
----------
frequency : float
Frequency [Hz] of the input signal.
Returns
-------
Potential
Data object representing the potential of volume conductor and
floating values of floating contacts.
"""
self._frequency = frequency
_logger.debug("Get conductivity at frequency")
self._sigma = self.conductivity_cf(self.mesh, frequency)
_logger.debug(f"Sigma: {self._sigma}")
# update boundary condition values
_logger.debug("Assign potential values")
boundary_values = self.contacts.voltages
coefficient = self.mesh.boundary_coefficients(boundary_values)
self._potential.Set(coefficient=coefficient, VOL_or_BND=ngsolve.BND)
_logger.debug("Prepare weak form")
u = self._space.TrialFunction()
v = self._space.TestFunction()
_logger.debug("Bilinear form")
# TODO symmetric even if anisotropy?
bilinear_form = ngsolve.BilinearForm(space=self._space)
_logger.debug("Bilinear form, formulation")
bilinear_form += self._sigma * ngsolve.grad(u) * ngsolve.grad(v) * ngsolve.dx
_logger.debug("Linear form")
linear_form = ngsolve.LinearForm(space=self._space)
# add surface impedance as Robin BC
self._surface_impedances = self.contacts.get_surface_impedances(
frequency, is_complex=self.is_complex
)
for contact in self.contacts.active:
if contact.surface_impedance_model is None:
continue
ys = 1.0 / self._surface_impedances[contact.name]
_logger.debug(
f"Contact: {contact.name}, Ys: {ys},"
f" Boundary Value: {boundary_values[contact.name]}"
)
bilinear_form += ngsolve.CF(ys) * u * v * ngsolve.ds(contact.name)
linear_form += (
ngsolve.CF(ys)
* ngsolve.CF(boundary_values[contact.name])
* v
* ngsolve.ds(contact.name)
)
_logger.debug("Solve BVP")
self.solver.bvp(bilinear_form, linear_form, self._potential)