Effect of discretization error and adaptive mesh generation in diffuse optical absorption imaging: II

In diffuse optical tomography (DOT), the discretization error in the numerical solutions of the forward and inverse problems results in error in the reconstructed optical images. In this first part of our work, we analyse the error in the reconstructed optical absorption images, resulting from the discretization of the forward and inverse problems. Our analysis identifies several factors which influence the extent to which the discretization impacts on the accuracy of the reconstructed images. For example, the mutual dependence of the forward and inverse problems, the number of sources and detectors, their configuration and their orientation with respect to optical absorptive heterogeneities, and the formulation of the inverse problem. As a result, our error analysis shows that the discretization of one problem cannot be considered independent of the other problem. While our analysis focuses specifically on the discretization error in DOT, the approach can be extended to quantify other error sources in DOT and other inverse parameter estimation problems.

We consider the inverse electrocardiographic problem of computing epicardial potentials from a body-surface potential map. We study how to improve numerical approximation of the inverse problem when the finite-element method is used. Being ill-posed, the inverse problem requires different discretization strategies from its corresponding forward problem. We propose refinement guidelines that specifically address the ill-posedness of the problem. The resulting guidelines necessitate the use of hybrid finite elements composed of tetrahedra and prism elements. Also, in order to maintain consistent numerical quality when the inverse problem is discretized into different scales, we propose a new family of regularizers using the variational principle underlying finite-element methods. These variational-formed regularizers serve as an alternative to the traditional Tikhonov regularizers, but preserves the L(2) norm and thereby achieves consistent regularization in multiscale simulations. The variational formulation also enables a simple construction of the discrete gradient operator over irregular meshes, which is difficult to define in traditional discretization schemes. We validated our hybrid element technique and the variational regularizers by simulations on a realistic 3-D torso/heart model with empirical heart data. Results show that discretization based on our proposed strategies mitigates the ill-conditioning and improves the inverse solution, and that the variational formulation may benefit a broader range of potential-based bioelectric problems.

In the last 25 years, mathematical and numerical modeling of cardiovascular problems has experienced a spectacular progress. Moving progressively from the stage of idealized models supporting traditional in vitro and in vivo approaches – basically with proofs of concept – Cardiovascular Mathematics has been transforming into a unique tool for patient-specific analysis. This has been made possible by the strong integration with imaging devices. Nowadays, the work pipeline starting from a medical image in the Digital Imaging and Communications in Medicine (DICOM) format to retrieve the patient-specific geometry of an artery or the left ventricle to be used as computational domain in a numerical simulation is well established. This has initiated a paradigm shift in the interactions between mathematicians, bioengineers, and medical doctors. The latter are progressively getting used to in silico analysis as an important component of their research [3]. The introduction of numerical procedures as a part of an established clinical routine and more in general of a consolidated support to the decision-making process of physicians still requires some steps both in terms of infrastructures (to bring computational tools to the operating room or the bedside) and methods. In particular, the quality of the numerical results needs to be assessed and certified. The reliability of simulations calls for an accurate quantification and possibly reduction of uncertainty. In this scenario and in view of terrific advancements of measuring techniques, an important research line – quite established in other fields – is data assimilation [4,6,9], that is, the integration of numerical simulations and measurements. We may say that numerical models provide a background knowledge (based on physical principles and constitutive laws, not patient-specific), whereas measures give a foreground (individual) information; accuracy of in silico procedures relies on the correct integration of these two levels of knowledge of the problem [9,10]. As a matter of fact, numerical models depend on parameters and boundary/initial conditions. Their knowledge is in practice usually incomplete, inaccurate, and noisy. In some cases, the sensitivity itself of the final results to these parameters is unclear and needs to be quantified. Such an uncertainty may be significantly reduced by the availability of patient-specific measures; on the other hand, the quality of measures can be strongly enhanced by comparison with mathematical models (after all, noise is in general not educated, it does not know PDEs . . . ) [4,5,7,8,9,11]. We think therefore that, in the future, a crucial step for the clinical application of Cardiovascular Mathematics will be the passage from numerical simulation to data assimilation. To be concrete, borrowing a motto by Michel Fortin supporting mesh adaptivity in finite element simulations, stating that a ‘grid should not be considered a data, but an unknown of the problem’, we think that many parameters currently tuned on the basis of nonpatient specific criteria (even in patient-specific geometries) should be considered part of the unknowns to be found by assimilation procedures.

Numerical experiments were performed to assess the effect of numerical discretization error on the convergence rate of polynomial chaos (PC) approximations for a transonic axial compressor stage. A random variable with a uniform distribution and expected value of one was introduced into the expression for turbulent viscosity of the k-ω SST turbulence model. Model uncertainty was quantified from the expected value and standard deviation estimates obtained via univariate non-intrusive polynomial chaos. Spectral projection and point collocation were both used and their results were compared. The effect of discretization error on convergence of the PC approximation was investigated using a grid refinement study with four grids. The PC expansion was computed for each grid while maintaining the same boundary conditions, basis functions, model evaluations, random variable distribution, and polynomial order. The quantities of interest (QOIs) were total–to–total pressure ratio, total–to–total temperature, and adiabatic efficiency. The grid resolution was found to have an influence on resulting surrogate models and the estimates of expected value and standard deviation for all QOIs. However, the estimates converged towards final values as the mesh was refined. Point collocation provided different estimates from spectral projection and the difference was also found to depend on the mesh size.

Inverse problem solution is an integral part of data interpretation for well testing in petroleum reservoirs. In case of two-phase well tests with water injection, forward problem is based on the multiphase flow model in porous media and solved numerically. The inverse problem is based on a misfit or likelihood objective function. Adjoint methods have proved robust and efficient for gradient calculation of the objective function in this type of problems. However, if time-lapse electrical resistivity measurements during the well test are included in the objective function, both the forward and inverse problems become multiphysical, and straightforward application of the adjoint method is problematic. In this paper we present a novel adjoint algorithm for the inverse problems considered. It takes into account the structure of cross dependencies between flow and electrical equations and variables, as well as specifics of the equations (mixed parabolic-hyperbolic for flow and elliptic for electricity), numerical discretizations and grids, and measurements in the inverse problem. Decomposition is proposed for the adjoint problem which makes possible step-wise solution of the electric adjoint equations, like in the forward problem, after which a cross-term is computed and added to the right-hand side of the flow adjoint equations at this timestep. The overall procedure provides accurate gradient calculation for the multiphysical objective function while preserving robustness and efficiency of the adjoint methods. Example cases of the adjoint gradient calculation are presented and compared to straightforward difference-based gradient calculation in terms of accuracy and efficiency.

BackgroundThe aim of electroencephalogram (EEG) source localization is to find the brain areas responsible for EEG waves of interest. It consists of solving forward and inverse problems. The forward problem is solved by starting from a given electrical source and calculating the potentials at the electrodes. These evaluations are necessary to solve the inverse problem which is defined as finding brain sources which are responsible for the measured potentials at the EEG electrodes.MethodsWhile other reviews give an extensive summary of the both forward and inverse problem, this review article focuses on different aspects of solving the forward problem and it is intended for newcomers in this research field.ResultsIt starts with focusing on the generators of the EEG: the post-synaptic potentials in the apical dendrites of pyramidal neurons. These cells generate an extracellular current which can be modeled by Poisson's differential equation, and Neumann and Dirichlet boundary conditions. The compartments in which these currents flow can be anisotropic (e.g. skull and white matter). In a three-shell spherical head model an analytical expression exists to solve the forward problem. During the last two decades researchers have tried to solve Poisson's equation in a realistically shaped head model obtained from 3D medical images, which requires numerical methods. The following methods are compared with each other: the boundary element method (BEM), the finite element method (FEM) and the finite difference method (FDM). In the last two methods anisotropic conducting compartments can conveniently be introduced. Then the focus will be set on the use of reciprocity in EEG source localization. It is introduced to speed up the forward calculations which are here performed for each electrode position rather than for each dipole position. Solving Poisson's equation utilizing FEM and FDM corresponds to solving a large sparse linear system. Iterative methods are required to solve these sparse linear systems. The following iterative methods are discussed: successive over-relaxation, conjugate gradients method and algebraic multigrid method.ConclusionSolving the forward problem has been well documented in the past decades. In the past simplified spherical head models are used, whereas nowadays a combination of imaging modalities are used to accurately describe the geometry of the head model. Efforts have been done on realistically describing the shape of the head model, as well as the heterogenity of the tissue types and realistically determining the conductivity. However, the determination and validation of the in vivo conductivity values is still an important topic in this field. In addition, more studies have to be done on the influence of all the parameters of the head model and of the numerical techniques on the solution of the forward problem.

Researchers have investigated the dynamics of oscillation of a single-link assistive manipulator and carried out vibration analysis of the same due to its great application in real life. Thus, studying the dynamics of flexible assistive robots and developing the model is indeed an interesting research topic. It is possible that the exact solution to the equations associated with such a single-link manipulator does not always exist. We used a numerical method, viz., Galerkin’s method to solve the governing differential equations of the problem, and the corresponding frequency of oscillation was calculated. This is known as the forward problem. On the other hand, the inverse problem is to compute values of the unknowns when the frequency is known by experiment, but the values of the remaining parameters are not known. With this backdrop, the purpose of this chapter is to investigate forward and inverse problems related to single-link robotic manipulators using Euler’s beam theory to define the system dynamics. Firstly, the frequencies of the first four modes of vibration have been calculated. Furthermore, the frequency is assumed to be known (experimental outcome), and the inverse case is studied using Young’s modulus and density of the material of manufacturing of the manipulator as unknown parameters. Besides, the system of single-link manipulators has been studied in uncertain environments using fuzzy theory for forward and inverse problems.

International Journal of Computational Engineering ScienceVol. 04, No. 01, pp. 19-65 (2003) No AccessFULLY NONLINEAR MODELING AND ANALYSIS OF PRECISION MEMBRANESP. FRANK PAI and LEYLAND G. YOUNGP. FRANK PAIDepartment of Mechanical and Aerospace Engineering, University of Missouri-Columbia, Columbia, MO 65211, USA Search for more papers by this author and LEYLAND G. YOUNGDepartment of Mechanical and Aerospace Engineering, University of Missouri-Columbia, Columbia, MO 65211, USA Search for more papers by this author https://doi.org/10.1142/S1465876303000788Cited by:4 PreviousNext AboutSectionsPDF/EPUB ToolsAdd to favoritesDownload CitationsTrack CitationsRecommend to Library ShareShare onFacebookTwitterLinked InRedditEmail AbstractHigh precision membranes are used in many current space applications. This paper presents a fully nonlinear membrane theory with forward and inverse analyses of high precision membrane structures. The fully nonlinear membrane theory is derived from Jaumann strains and stresses, exact coordinate transformations, the concept of local relative displacements, and orthogonal virtual rotations. In this theory, energy and Newtonian formulations are fully correlated, and every structural term can be interpreted in terms of vectors. Fully nonlinear ordinary differential equations (ODEs) governing the large static deformations of known axisymmetric membranes under known axisymmetric loading (i.e., forward problems) are presented as first-order ODEs, and a method for obtaining numerically exact solutions using the multiple shooting procedure is shown. A method for obtaining the undeformed geometry of any axisymmetric membrane with a known inflated geometry and a known internal pressure (i.e., inverse problems) is also derived. Numerical results from forward analysis are verified using results in the literature, and results from inverse analysis are verified using known exact solutions and solutions from the forward analysis. Results show that the membrane theory and the proposed numerical methods for solving nonlinear forward and inverse membrane problems are accurate.Keywords:Nonlinear membrane theorynumerically exact analysis References M. Salama , M. Lou and H. Fang , Proceedings of the AIAA Structures, Structural Dynamics and Materials Conference ( 2000 ) . Google Scholar M. S. Smith , R. S. Schallenkamp and L. G. Seely , AIAA International Balloon Technology Conference . Google Scholar D. Gorinevsky , T. Hyde and C. Cabuz , Proceedings of the AIAA Structures, Structural Dynamics and Materials Conference ( 2001 ) . Google Scholar V. Firt , Statics, Formfinding and Dynamics of Air-Supported Membrane Structures ( Martinus Nijhoff Publishers , Boston , 1983 ) . Google Scholar A. L. Palisoc, Inflatable Reflector Development Program, Task 3 Report, L'Garde Technical Report, LTR-94-AP-008, 1994 . Google Scholar S. V. Damle et al. , AIAA International Balloon Technology Conference . Google ScholarP. F. Pai, A. N. Palazotto and J. M. Greer Jr., Comput. Struct. 66, 823 (1998), DOI: 10.1016/S0045-7949(98)00004-2. Crossref, Google Scholar H. Furuya and Y. Kawasaki , Proceedings of the AIAA Structures, Structural Dynamics and Materials Conference . Google Scholar G. Greschik , M. M. Mikulas and A. Palisoc , Proceedings of the AIAA Structures, Structural Dynamics and Materials Conference ( 1998 ) . Google ScholarJ. Ruze, Proceedings of the IEEE 54, 633 (1966), DOI: 10.1109/PROC.1966.4784. Crossref, Google Scholar A. Palisoc et al. , Fifth SPIE Conference on Space Telescopes and Instruments . Google ScholarP. F. Pai and A. H. Nayfeh, Nonlinear Dynamics 6, 459 (1994). Crossref, Google ScholarP. F. Pai and A. N. Palazotto, J. Eng. Mech. 121, 568 (1995), DOI: 10.1061/(ASCE)0733-9399(1995)121:4(568). Crossref, Google ScholarP. F. Pai and A. H. Nayfeh, Nonlinear Dynamics 2, 445 (1991), DOI: 10.1007/BF00045438. Crossref, Google Scholar K. Washizu , Variational Methods in Elasticity & Plasticity , 3rd edn. ( Pergamon Press Inc , New York , 1982 ) . Google ScholarP. F. Pai and A. N. Palazotto, Int. J. Solids Struct. 32, 3047 (1995), DOI: 10.1016/0020-7683(94)00273-Y. Crossref, Google Scholar H. Golden , T. W. Strganac and R. A. Schapery , AIAA International Balloon Technology Conference ( 1997 ) . Google Scholar J. Stoer and R. Bulirsch , Introduction to Numerical Analysis ( Springer-Verlag , New York , 1980 ) . Crossref, Google Scholar FiguresReferencesRelatedDetailsCited By 4Analytical Investigation of Dynamics of Inflatable Parabolic Membrane ReflectorFushou Liu and Dongping Jin1 Jan 2015 | Journal of Spacecraft and Rockets, Vol. 52, No. 1A point collocation method for geometrically nonlinear membranesKyle F. Kolsti and Donald L. Kunz1 Jan 2013 | International Journal of Solids and Structures, Vol. 50, No. 2Numerical and experimental dynamic characteristics of thin-film membranesLeyland G. Young, Suresh Ramanathan, Jiazhu Hu and P. Frank Pai1 May 2005 | International Journal of Solids and Structures, Vol. 42, No. 9-10Modeling, Analysis and Testing of some Deployable/Inflatable StructuresPerngjin Pai and Leyland Young26 June 2012 Recommended Vol. 04, No. 01 Metrics History Received 26 August 2002 Accepted 4 November 2002 KeywordsNonlinear membrane theorynumerically exact analysisPDF download

'Acoustic probes have shown promise to be quite effective in determining the solid content in solid-liquid suspensions. However, the presence of small amounts of gas in the waste slurries stored in tanks across the DOE complex prevents straightforward application for characterization of these slurries. The proposed research will develop an acoustic probe for monitoring particle size and volume fraction in slurries in the absence and the presence of gas bubbles. Theoretical Analysis Accomplished: Attenuation of sound waves depends on the size distribution of the solids and the volume fraction of solids. These can in principle be calculated from attenuation measured over a range of frequencies. However, small amounts of bubbles distort the measured attenuation. A typical result from theoretical analysis for the attenuation of solid- gas-liquid systems is given in Figure 1. The total attenuation of a sound wave v(o) equals the sum of contributions by a large number of ''bins'' of particle sizes. This notion yields the following equation for the (hitherto) unknown number density of solid particles as a function of particle radius N(a): j k(o,a)N(a)da = v(o), where the kernel k(o,a) is obtained from analysis. If N(a) is given, the above equation is used to calculate the attenuation v(o). This is referred to as solving the ''forward problem''. Solving for N(a) with v(o) given is the ''inverse problem''. A complication that one faces when trying to solve the inverse problem is that the stated problem is mathematically ill-posed, i.e., small fluctuations in v(o) cause large fluctuations in the result for the number density. Therefore the problem needs to be ''regularized'', i.e., the stated problem needs to be changed slightly such as to make it well-posed. This has been done by others for gas-liquid systems in the past. This approach is currently being applied in the present project to solid-liquid systems. As is shown in Figure 2, it successfully recovers the number density that has been used in the forward problem to generate attenuation data. Having this solution technique giving reliable results for the inverse problems of both gas-liquid and solid-liquid systems, the authors shall apply this method in the near future to solid-gas-liquid systems.'

Many materials (e.g., rubber or biologic tissues) are "nearly" incompressible and often assumed to be incompressible in their constitutive equations. This assumption hinders realistic analyses of wave motion including acoustoelasticity. In this study, this constraint is relaxed and the reflected waves from nearly incompressible, hyper-elastic materials are examined. Specifically, reflection coefficients are considered from the interface of water and uni-axially prestretched rubber. Both forward and inverse problems are experimentally and analytically studied with the incident wave perpendicular to the interface. In the forward problem, the wave reflection coefficient at the interface is evaluated with strain energy functions for nearly incompressible materials in order to compute applied strain. For the general inverse problem, mathematical relations are derived that identify both uni-axial strains and normalized material constants from reflected wave data. The validity of this method of analysis is demonstrated via an experiment with stretched rubber. Results demonstrate that applied strains and normalized material coefficients can be simultaneously determined from the reflected wave data alone if they are collected at several different (but unknown) levels of strain. This study therefore indicates that acoustoelasticity, with an appropriate constitutive formulation, can determine strain and material properties in hyper-elastic, nearly incompressible materials.

Recent work on agent-based models of wealth distribution has yielded nonlinear, non-local Fokker-Planck equations whose steady-state solutions describe empirical wealth distributions with remarkable accuracy using only a few free parameters. Because these equations are often used to solve the 'inverse problem' of determining the free parameters given empirical wealth data, there is much impetus to find fast and accurate methods of solving the 'forward problem' of finding the steady state corresponding to given parameters. In this work, we derive and calibrate a lattice Boltzmann equation for this purpose. This article is part of the theme issue 'Fluid dynamics, soft matter and complex systems: recent results and new methods'.

PyTorch-FEA: Autograd-enabled finite element analysis methods with applications for biomechanical analysis of human aorta

Diffuse optical tomography (DOT) in the near infrared involves reconstruction of spatially varying optical properties of turbid medium from boundary measurements based on a forward model of photon propagation. Due to highly non-linear nature of the DOT, high quality image reconstruction is a computationally demanding problem that requires repeated solutions of both the forward and the inverse problems. Therefore, it is highly desirable to develop methods and algorithms that are computationally efficient. In this paper, we propose a domain decomposition approach to address the computational complexity of the DOT problem. We propose a two-level multiplicative overlapping domain decomposition method for the forward problem and a two-level space decomposition method for the inverse problem. We showed the convergence of the inverse solver and derived the computational complexity of each method. We demonstrate the performance of the proposed approach in numerical simulations.

We compute the solution to the photon diffusion equation using a state-space system. The state-space system under consideration is obtained by applying a FEM technique to the photon diffusion equation. Computing the photon density in the medium is commonly called in computational tomography, the forward problem. We investigate and present some preliminary results on the inverse problem. The inverse problem entails the identification of anomalies (benign or malignant) inside the medium. We investigate the inverse problem using an optimization setup. Specifically, we use a baseline homogeneous experiment that represents the healthy tissue to compare with the actual data. We then form an error signal between the heterogeneous and homogeneous case, which we try to minimize over the absorption coefficient of the true medium. We minimize it using the two-norm of the error signal.

The inverse problem for conducting defective lattices