Re: [deal.II] How to compute derivatives of the composite function with the reference cell mapping
Thank you Wolfgang, this does look relevant! /Anton
On Sat, Aug 5, 2023 at 12:37 PM Wolfgang Bangerth
wrote:
>
> On 8/5/23 04:30, Anton Evgrafov wrote:
> > I have a problem involving a
> > field living in R^{dim^2} and a field in R^{dim}.
>
> This may or may not be relevant to you, but it might be worth taking a look at
> hyper.deal:
>https://github.com/hyperdeal/hyperdeal
>
> Best
> W.
>
> --
>
> Wolfgang Bangerth email: bange...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to a topic in the Google
> Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/dealii/FONV6StZZus/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/0f072f81-1297-fa96-4744-d4f22908b19b%40colostate.edu.
--
Anton Evgrafov
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/CAFBVMaqjwpaKYfenaQ2kSe7BSEcZ1TsLuiV6wigx3K2ERgubfA%40mail.gmail.com.
Re: [deal.II] How to compute derivatives of the composite function with the reference cell mapping
Simon,
> That sounds doable and cheaper. If you are going for this, FEFieldFunction
> might be useful:
> https://www.dealii.org/developer/doxygen/deal.II/classFunctions_1_1FEFieldFunction.html
I have now compared the implementation using UnitToReal function
wrapper + QuadratureGenerator vs using an auxiliary triangulation with
one reference cell onto which I interpolate my "level set function",
so to speak, + DiscreteQuadratureGenerator. The latter is waaay
faster.
> Just out of curiosity (and if it is not a research secret), what is your use
> case? Why do you need to generate many quadratures for the same cell?
Well, I am certainly misusing deal.II. I have a problem involving a
field living in R^{dim^2} and a field in R^{dim}. The weak form
involves products between the two, and integrals are over
dim^2-dimensional domains, which may be decomposed into two iterated
dim-dimensional integrals. The inner integral then goes over a ball
which cuts elements. This is where non-matching quadrature comes
along, but I have to compute them for each outer quadrature point, so
to speak. Also the integrand is not polynomial all, so some
quadratures are of higher degree.
The proper way of doing this would be to utilize a FEM library that
supports mesh/elements/quadratures in dim^2.Can one do this with
deal.ii? Most of the code seems to be dimension-independent, but then
there are plenty of places where e.g. Point constructors are
specialized for dimensions 1..3.
Right now I am doing a "quick and dirty" test implementation, where I
am basically manually working on a tensor product mesh with tensor
product FEM etc, which is highly inelegant. If you are
doing/interested in research, I would be happy to discuss the problem
in detail/collaborate.
/Anton
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/CAFBVMarTRWobVTMcf7NRO2rCqiDca9eTKQM7QA%3DXxcWP%3DKggAQ%40mail.gmail.com.
Re: [deal.II] How to compute derivatives of the composite function with the reference cell mapping
Dear Simon,
Thank you very much for your reply. The wrapper code looks
embarrassingly simple!
> If you are working with a perfectly Cartesian background I would expect that
> it would actually be faster to generate the quadrature in real space and then
> transform the quadrature to reference space before passing it to FEValues. If
> you want to do this the functions cell->bounding_box() and
> BoundingBox::real_to_unit(..) are useful.
That I did not think about for some reason :)
If my element is not aligned with the axes, won't I risk getting
quadrature points outside of the element though? Or is this what you
meant by the "perfect Cartesian" situation? Also I am not 100% clear
at the moment about how to change the quadrature weights after the
physical->unit transformation, as I would need access to Jacobian
determinants.
Another option I was contemplating is to apply the
NonMatching::DiscreteQuadratureGenerator on a fake triangulation
with just one element. Then I simply would need to interpolate
(evaluate) my real-space function at this one element and let the
class' implementation do the magic for me. This is probably cheaper
and about equally accurate?
/Anton
On Fri, Aug 4, 2023 at 12:07 PM Simon Sticko wrote:
>
> Hi.
>
> I have done this in the past. I dug up the wrapper function I had for it, in
> case we decide that we should add it to the library:
>
> https://github.com/dealii/dealii/pull/15838
>
> The problem with this approach is that evaluating the reference space level
> set function gets very expensive. QuadratureGenerator uses many function
> calls and we have to use the Jacobian and its derivative in the
> transformation. This makes this approach significantly slower than using an
> interpolated level set function.
>
> If you are working with a perfectly Cartesian background I would expect that
> it would actually be faster to generate the quadrature in real space and then
> transform the quadrature to reference space before passing it to FEValues. If
> you want to do this the functions cell->bounding_box() and
> BoundingBox::real_to_unit(..) are useful.
>
> Best,
> Simon
>
> On 03/08/2023 23:11, Anton wrote:
> > Hello,
> >
> > I would like to use the NonMatching::QuadratureGenerator class
> > directly, as I need to generate multiple (as in many!, unfortunately)
> > nonmatching quadratures for the same cell. Therefore I would like to avoid
> > the route of interpolating the level set function onto the whole mesh etc,
> > as described in CutFEM tutorial.
> >
> > Therefore I am thinking about calling the "generate" function, which only
> > needs a level set function and a bounding box as inputs. According to the
> > documentation the level set function should provide the gradient and the
> > Hessian. I would like to do this in the reference cell coordinates, so
> > that I can simply use this quadrature later via FEView as one normally does.
> > * The bounding box I presume is the hypercube [0,1]^{dim} for the meshes
> > I am interested in. (Can one request them for the reference cell and not
> > the triangulation cell?)
> >
> > * The level set function should also be reasonably straightforward. I
> > know the function in the physical coordinates (say, f(x)) and all its
> > derivatives. I can also get the map from reference to physical coordinates
> > via FEView, let us call this map x=M(y). The issue is how to extract the
> > derivatives of this map so that the chain rule can be applied? I am mostly
> > working with Q1 maps, so of course I could figure out the derivatives "by
> > hand". I am just wondering if there is a more intelligent (code-wise) way
> > of doing this?
> >
> > --
> > The deal.II project is located at http://www.dealii.org/
> > <http://www.dealii.org/>
> > For mailing list/forum options, see
> > https://groups.google.com/d/forum/dealii?hl=en
> > <https://groups.google.com/d/forum/dealii?hl=en>
> > ---
> > You received this message because you are subscribed to the Google Groups
> > "deal.II User Group" group.
> > To unsubscribe from this group and stop receiving emails from it, send an
> > email to dealii+unsubscr...@googlegroups.com
> > <mailto:dealii+unsubscr...@googlegroups.com>.
> > To view this discussion on the web visit
> > https://groups.google.com/d/msgid/dealii/4ac38563-e467-4f1d-8cb4-ee35379515a6n%40googlegroups.com
> >
> > <https://groups.google.com/d/msgid/dealii/4ac38563-e467-4f1d-8cb4-ee35379515a6n%40googlegroups.com?utm_medium=email_source=footer>.
>
> --
> The deal.II project is
[deal.II] How to compute derivatives of the composite function with the reference cell mapping
Hello,
I would like to use the NonMatching::QuadratureGenerator class
directly, as I need to generate multiple (as in many!, unfortunately)
nonmatching quadratures for the same cell. Therefore I would like to avoid
the route of interpolating the level set function onto the whole mesh etc,
as described in CutFEM tutorial.
Therefore I am thinking about calling the "generate" function, which only
needs a level set function and a bounding box as inputs. According to the
documentation the level set function should provide the gradient and the
Hessian. I would like to do this in the reference cell coordinates, so
that I can simply use this quadrature later via FEView as one normally does.
* The bounding box I presume is the hypercube [0,1]^{dim} for the meshes I
am interested in. (Can one request them for the reference cell and not the
triangulation cell?)
* The level set function should also be reasonably straightforward. I
know the function in the physical coordinates (say, f(x)) and all its
derivatives. I can also get the map from reference to physical coordinates
via FEView, let us call this map x=M(y). The issue is how to extract the
derivatives of this map so that the chain rule can be applied? I am mostly
working with Q1 maps, so of course I could figure out the derivatives "by
hand". I am just wondering if there is a more intelligent (code-wise) way
of doing this?
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/4ac38563-e467-4f1d-8cb4-ee35379515a6n%40googlegroups.com.
Re: [deal.II] Re: Eigenproblem and creating a preconditioner out of a linear operator, Eigensolver Selection
Dear Daniel, I am working on a similar problem (planetary acoustic oscillations). I am interested in looking at your tutorial for an electromagnetic cavity, but it seems that the pull request was deleted. I wonder if it can be revived. Thank you, Anton. On Monday, July 8, 2019 at 3:51:08 AM UTC-7 Daniel Garcia-Sanchez wrote: > Hi Andreas, > > > On Sunday, July 7, 2019 at 4:26:39 PM UTC+2, Andreas Hegendörfer wrote: >> >> I also tried a spectral transformation with the Arpack solver with the >> same result as without spectral transformation. I am interested in the >> smallest real eigenvalues. I know from previous calculations with the >> Krylov Schur solver form SLEPc that using or not using a spectral >> transformation makes a very big difference here. >> > > Yes, using a spectral transformation makes a huge difference if your > eigenvalues are from the interior of the spectrum, although if you are > interested in the smallest real eigenvalue, the difference might not be > that important. > > The drawback of an spectral transformation is that you have to use a > direct solver. > > You can use an iterative solver with an spectral transformation for very > large problems. However, using an iterative linear solver with an spectral > transformation makes the overall solution process less robust. > > Although the direct solver approach may seem too costly, the factorization > is only carried out at the beginning of the eigenvalue calculation and this > cost is amortized in each subsequent application of the operator. > > I think that the drawback of an iterative solver is that it requires lot > of memory and does not scale very well. This figure can give you an idea of > the scalability of MUMPS (direct parallel solver). > > https://www.researchgate.net/figure/Strong-scalability-of-an-iterative-solver-with-different-preconditioners-versus-MUMPS_fig2_282172435 > > >> I do not want to use SLEPc here because I think handling the PETSc >> Matrices and vectors is too uncomfortable for my application. Am I right at >> this point? What do you think about using SLEPc here? >> > > I'm writing a tutorial about how to calculate the eigenmodes of an > electromagnetic cavity. So far, I've written the code, I'm writing now the > documentation. You can take a look. You can find the code in this pull > request: > > https://github.com/dealii/dealii/pull/8345 > > The tutorial uses MPI, and SLEPc with std::complex. > > The eigenvalues in this tutorial are from the interior of the spectrum, > therefore I have to use an spectral transformation with a direct solver. > > Note that PETSc has to be compiled with MUMPS if you want to run that > tutorial. > > Best, > Daniel > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/5f867abb-b809-4230-994f-e5a920e95336n%40googlegroups.com.
Re: [deal.II] Vector eigenvalue problem with SLEPc
>>> *That might not easily be possible with our PETSc interfaces.* *PETSc has the concept of MatShell, which is in essence what deal.II'sLinearOperator is: A class that implements a matrix-vector operation,without giving access to individual matrix elements. You could trywhether you can implement something derived from MatrixBase that is aMatShell operation in the same way as the existing classes derived fromMatrixBase model other types of PETSc matrices.* *I don't know how difficult that would be, though.* Ok. I see. Seems like using ARPACK might be easier then. >>> You can do that, but I believe that it's also possible to work with the original system where one of the two matrices is only semidefinite. I have redone it with the original system, but it seems that PETSc is complaining about missing diagonal entry. *[0]PETSC ERROR: Object is in wrong state* *[0]PETSC ERROR: Matrix is missing diagonal entry 3* I presume this happens due to a zero block in one (or both) of the original matrices. Hmmm, should I then just explicitly see all zero diagonal values to 0.0. Anton. On Wednesday, August 25, 2021 at 12:20:38 PM UTC-6 Wolfgang Bangerth wrote: > On 8/25/21 11:53 AM, Anton Ermakov wrote: > > I I¯ A 0 ¯I*eigenvalue + I¯0 B ¯ I I * eigenvector =0 > > I I_ 0 0 _I I_ C D _I I > > ¯ ¯ > > where eigenvector consists of a vector displacement and a scalar > > pressure perturbation: > > > > eigenvector = transpose([U P]) > > > > It seems from the literature (e.g., Wang & Bathe 1997) that it is > > standard to rewrite this system to eliminate the pressure variable: > > > > eigenvalue * A * U - B * D^-1 * C * U = 0 > > > > or > > > > eigenvalue * A * U + M * U = 0 > > > > with M = - B * D^-1 * C and P = -D^-1 * C * U > > You can do that, but I believe that it's also possible to work with the > original system where one of the two matrices is only semidefinite. > > > What I am having trouble with is how to efficiently represent matrix > > products and inverses that go into matrix M. All the matrices (A, B, C, > > D) are sparse, but due to the inverse, it seems that M will be a full > > matrix and it might be unaffordable to store it in the memory. There was > > a discussion on a similar topic here: > > > > https://groups.google.com/g/dealii/c/5P_fv0zg7jg/m/CPYcYrbsDQAJ > > > > It seems that LinearOperator was used in that discussion to manipulate > > the matrix products and inverses, which seemed like an approach to > > follow if you use ARPACK eigensolvers. However, I am using SLEPc and it > > seems the inputs to SLEPc solvers would have to be objects of class > > PETScWrappers::MatrixBase > > < > https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassPETScWrappers_1_1MatrixBase.html=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Ce376b299e56c44690ada08d967f13673%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637655108014800275%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000=QQLBoH27O1sZ6GMXDatmucrxAye8edSeN3XNqAybONI%3D=0>. > > > > > > > > So, my question is how to prepare such a matrix M involving matrix > > products and matrix inverse efficiently and represent it as > > PETScWrappers::MatrixBase > > That might not easily be possible with our PETSc interfaces. > > PETSc has the concept of MatShell, which is in essence what deal.II's > LinearOperator is: A class that implements a matrix-vector operation, > without giving access to individual matrix elements. You could try > whether you can implement something derived from MatrixBase that is a > MatShell operation in the same way as the existing classes derived from > MatrixBase model other types of PETSc matrices. > > I don't know how difficult that would be, though. > > Best > W. > > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/88a0aebf-86dd-436f-8f56-8aef2edae4d8n%40googlegroups.com.
[deal.II] Vector eigenvalue problem with SLEPc
Dear Deal.II community, I am working on implementing a vector eigenproblem in deal.ii, starting from the step-36 example. The question I am studying is the normal oscillation of a fluid planet. The eigenproblem is of this type: _ _ I I¯ A 0 ¯I*eigenvalue + I¯0 B ¯ I I * eigenvector =0 I I_ 0 0 _I I_ C D _I I ¯ ¯ where eigenvector consists of a vector displacement and a scalar pressure perturbation: eigenvector = transpose([U P]) It seems from the literature (e.g., Wang & Bathe 1997) that it is standard to rewrite this system to eliminate the pressure variable: eigenvalue * A * U - B * D^-1 * C * U = 0 or eigenvalue * A * U + M * U = 0 with M = - B * D^-1 * C and P = -D^-1 * C * U What I am having trouble with is how to efficiently represent matrix products and inverses that go into matrix M. All the matrices (A, B, C, D) are sparse, but due to the inverse, it seems that M will be a full matrix and it might be unaffordable to store it in the memory. There was a discussion on a similar topic here: https://groups.google.com/g/dealii/c/5P_fv0zg7jg/m/CPYcYrbsDQAJ It seems that LinearOperator was used in that discussion to manipulate the matrix products and inverses, which seemed like an approach to follow if you use ARPACK eigensolvers. However, I am using SLEPc and it seems the inputs to SLEPc solvers would have to be objects of class PETScWrappers::MatrixBase <https://www.dealii.org/current/doxygen/deal.II/classPETScWrappers_1_1MatrixBase.html> . So, my question is how to prepare such a matrix M involving matrix products and matrix inverse efficiently and represent it as PETScWrappers::MatrixBase <https://www.dealii.org/current/doxygen/deal.II/classPETScWrappers_1_1MatrixBase.html> suitable for a SLEPc eigensolver. In addition, I am interested in the interior eigenvalues. It seems that shift-invert spectral transformation is available in SLEPc - this is what I plan to use. However, SLEPc also has a polynomial filtering spectral transformation. It does not seem that deal.ii provides a wrapper to that one though. So, I wonder if anyone has had expertise with making an interface in deal.ii to use SLEPc's polynomial filtering transformation. Thanks a lot, Anton. -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/426e3d15-a42f-4e90-a3f7-117c66e71c92n%40googlegroups.com.
[deal.II] Re: shift-invert transformation problem with slepc
Ok, I see - it is that EPS_TARGET_MAGNITUDE flag that was missing. It works now! Thanks a lot, Simon. Anton. On Tuesday, August 24, 2021 at 12:14:20 AM UTC-6 simon...@gmail.com wrote: > Hi, > > I think you need to set both which eigen pair and the target value. For > example, if you look for the eigenvalue closest in magnitude to 0, you need > to set the following > > const double target_eigenvalue = 0; > eigensolver.set_which_eigenpairs(EPS_TARGET_MAGNITUDE); > eigensolver.set_target_eigenvalue(target_eigenvalue); > > Depending on slepc version it might be that you can only choose > EPS_TARGET_MAGNITUDE. See question 8 here: > > https://slepc.upv.es/documentation/faq.htm > > Best, > Simon > > On Tuesday, August 24, 2021 at 3:39:29 AM UTC+2 anton.i...@gmail.com > wrote: > >> Dear All, >> >> I am trying to use a shift-invert transformation in step-36. I did >> minimal changes in the solve function adding the following code that >> describes what spectral transformation I need to do: >> >> SolverControl solver_control (dof_handler.n_dofs(), 1e-9); >> SLEPcWrappers::SolverKrylovSchur eigensolver(solver_control); >> eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL); >> eigensolver.set_problem_type (EPS_GHEP); >> >> double eigen_shift = 1.0e-4; >> SLEPcWrappers::TransformationShiftInvert::AdditionalData >> additional_data(eigen_shift); >> SLEPcWrappers::TransformationShiftInvert shift (MPI_COMM_WORLD, >> additional_data); >> eigensolver.set_transformation (shift); >> >> eigensolver.solve(stiffness_matrix, >> mass_matrix, >> eigenvalues, >> eigenfunctions, >> eigenfunctions.size()); >> >> and it gives me the following error: >> >> [0]PETSC ERROR: Shift-and-invert requires a target 'which' (see >> EPSSetWhichEigenpairs), for instance -st_type sinvert -eps_target 0 >> -eps_target_magnitude >> >> Any help to resolve this error is greatly appreciated. I also attached >> the full code of modified step-36. >> >> Thank you, >> Anton. >> > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/879b9d55-f2ae-456f-828b-91bb26b6ba20n%40googlegroups.com.
[deal.II] shift-invert transformation problem with slepc
Dear All,
I am trying to use a shift-invert transformation in step-36. I did minimal
changes in the solve function adding the following code that describes what
spectral transformation I need to do:
SolverControl solver_control (dof_handler.n_dofs(), 1e-9);
SLEPcWrappers::SolverKrylovSchur eigensolver(solver_control);
eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL);
eigensolver.set_problem_type (EPS_GHEP);
double eigen_shift = 1.0e-4;
SLEPcWrappers::TransformationShiftInvert::AdditionalData
additional_data(eigen_shift);
SLEPcWrappers::TransformationShiftInvert shift (MPI_COMM_WORLD,
additional_data);
eigensolver.set_transformation (shift);
eigensolver.solve(stiffness_matrix,
mass_matrix,
eigenvalues,
eigenfunctions,
eigenfunctions.size());
and it gives me the following error:
[0]PETSC ERROR: Shift-and-invert requires a target 'which' (see
EPSSetWhichEigenpairs), for instance -st_type sinvert -eps_target 0
-eps_target_magnitude
Any help to resolve this error is greatly appreciated. I also attached the
full code of modified step-36.
Thank you,
Anton.
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/59b3e662-f3c8-4ae5-b48d-f6dfad601153n%40googlegroups.com.
/* -
*
* Copyright (C) 2009 - 2019 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, 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.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* -
*
* Authors: Toby D. Young, Polish Academy of Sciences,
* Wolfgang Bangerth, Texas A University
*/
// @sect3{Include files}
// As mentioned in the introduction, this program is essentially only a
// slightly revised version of step-4. As a consequence, most of the following
// include files are as used there, or at least as used already in previous
// tutorial programs:
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
// IndexSet is used to set the size of each PETScWrappers::MPI::Vector:
#include
// PETSc appears here because SLEPc depends on this library:
#include
#include
// And then we need to actually import the interfaces for solvers that SLEPc
// provides:
#include
#include
// We also need some standard C++:
#include
#include
// Finally, as in previous programs, we import all the deal.II class and
// function names into the namespace into which everything in this program
// will go:
namespace Step36
{
using namespace dealii;
// @sect3{The EigenvalueProblem class template}
// Following is the class declaration for the main class template. It looks
// pretty much exactly like what has already been shown in step-4:
template
class EigenvalueProblem
{
public:
EigenvalueProblem(const std::string _file);
void run();
private:
void make_grid_and_dofs();
void assemble_system();
unsigned int solve();
void output_results() const;
Triangulation triangulation;
FE_Q fe;
DoFHandlerdof_handler;
// With these exceptions: For our eigenvalue problem, we need both a
// stiffness matrix for the left hand side as well as a mass matrix for
// the right hand side. We also need not just one solution function, but a
// whole set of these for the eigenfunctions we want to compute, along
// with the corresponding eigenvalues:
PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix;
std::vector eigenfunctions;
std::vector eigenvalues;
// And then we need an object that will store several run-time parameters
// that we will specify in an input file:
ParameterHandler parameters;
// Finally, we will have an object that contains "constraints" on our
// degrees of freedom. This could include hanging node constraints if we
// had adaptively refined meshes (which we don't have in the current
// program). Here, we will store the constraint
Re: [deal.II] treatment of handing nodes in the heat conduction problem
Yes, that part now works. Thank you, Wolfgang. Anton. On Wednesday, March 17, 2021 at 3:13:17 PM UTC-7 Wolfgang Bangerth wrote: > On 3/17/21 10:41 PM, Anton Ermakov wrote: > > Thanks for the reply. I think I got it. Basically, I had to manually > place the > > local matrices into the global matrix, without taking into account the > hanging > > node constraints, which are later taken care of by > > > > /constraints.condense(system_matrix, system_rhs);/ > > > > In the case of the mass matrix, I had to use this: > > > > /for (unsigned int i = 0; i < dofs_per_cell; ++i)/ > > > > /for (unsigned int j = 0; j < dofs_per_cell; ++j)/ > > > > > /mass_matrix.add(local_dof_indices[i], local_dof_indices[j], cell_matrix(i, > j));/ > > > > Instead of > > > > /constraints.distribute_local_to_global(cell_matrix, local_dof_indices, > > mass_matrix);/ > > Ah yes, that makes sense. So do I understand right that whatever you have > now > actually works? > > Best > W. > > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/fdc7b151-885e-48d9-aace-8e25e8669597n%40googlegroups.com.
Re: [deal.II] treatment of handing nodes in the heat conduction problem
Thanks for the reply. I think I got it. Basically, I had to manually place the local matrices into the global matrix, without taking into account the hanging node constraints, which are later taken care of by *constraints.condense(system_matrix, system_rhs);* In the case of the mass matrix, I had to use this: * for (unsigned int i = 0; i < dofs_per_cell; ++i)* * for (unsigned int j = 0; j < dofs_per_cell; ++j)* * mass_matrix.add(local_dof_indices[i], local_dof_indices[j], cell_matrix(i, j));* Instead of * constraints.distribute_local_to_global(cell_matrix, local_dof_indices, mass_matrix);* Anton. On Tuesday, March 16, 2021 at 3:22:54 AM UTC-7 Wolfgang Bangerth wrote: > On 3/16/21 6:12 AM, Anton Ermakov wrote: > > > > I attach the minimally changed Step-26 code to illustrate the problem. > It is > > different from the original step-26 in three important places: > > > > 1) in the run function, the mesh is refined in a corner > > 2) in the setup_system the mass and Laplace matrices are manually > created. > > 3) Initial temperature field is set to 100 K and the boundary condition > is > > temperature fixed at 50 K. No source term. > > I don't immediately see a problem. What happens if you use your matrices > with > the exact set up of step-26? Simplify your problem: Make one change at a > time, > not three. > > If you do get different results, compare the matrix you compute with the > one > computed by the MatrixTools functions. How are they different? And why? > > Best > W. > > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/665c688f-fca3-4c21-b56b-7193a09cdb47n%40googlegroups.com.
Re: [deal.II] project_boundary_values for components
Wolfgang, Thank you for the help. I think I have it working in the non-distributed case, but project_boundary_values segfaults on me without any explanations (compiled in debug mode) when run with a distributed triangulation. Is this "normal"? /Anton On Monday, April 8, 2019 at 12:11:19 AM UTC+2, Wolfgang Bangerth wrote: > > On 4/6/19 12:29 PM, Anton wrote: > > > > Basically if I interpolate a step function (which can be exactly > presented in > > for example DGQ1 space if the discontinuity is at the node) I get a > 1-element > > wide transition. In ASCII pictures, I expect ___ but I get > ___/--- > > instead :) > > Ah, yes -- ASCII art to the fore! Yes, the problem is that the DGQ(k) > elements > have support (interpolation) points at the vertices, and so you'll get the > same value for both cells at these vertices. That's unavoidable. > Projection is > the solution. > > > > So if I understand you correctly, one would need to (1) obtain > > ConstraintMatrix from project_boundary_values; (2) go over each row of > the > > constraint matrix and check which component it is; (3) if it is the > "right" > > component then copy the row. Could I ask how does one figure out the > > component number from the row number in the ConstraintMatrix - I cannot > > immediately figure this out from DoFHandler interface? > > You can't do it this way, but there are functions in DoFTools that allow > you > to extract *all* DoF indices for a particular component from a DoFHandler, > and > then you can cross-check the two objects. > > Best > W. > > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > > www: http://www.math.colostate.edu/~bangerth/ > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] project_boundary_values for components
Thank you Wolfgang. Basically if I interpolate a step function (which can be exactly presented in for example DGQ1 space if the discontinuity is at the node) I get a 1-element wide transition. In ASCII pictures, I expect ___ but I get ___/--- instead :) So if I understand you correctly, one would need to (1) obtain ConstraintMatrix from project_boundary_values; (2) go over each row of the constraint matrix and check which component it is; (3) if it is the "right" component then copy the row. Could I ask how does one figure out the component number from the row number in the ConstraintMatrix - I cannot immediately figure this out from DoFHandler interface? /Anton On Saturday, April 6, 2019 at 5:11:09 PM UTC+2, Wolfgang Bangerth wrote: > > On 4/5/19 11:25 AM, Anton wrote: > > I need to apply a boundary condition (discontinuous function) to a few > > components of a dof_handler only. I am currently using > > > > | > > voidVectorTools::interpolate_boundary_values > > < > https://dealii.org/current/doxygen/deal.II/group__constraints.html#ga2376b9a282a2b62be5c55ab62a11a132>(constDoFHandlerType > > > >,consttypes::boundary_id > > < > https://dealii.org/current/doxygen/deal.II/namespacetypes.html#aaf4eb6ec214fa642dfd956f11a9cd2d7>boundary_component,constFunction > > > > <https://dealii.org/current/doxygen/deal.II/classFunction.html> > > >_function,ConstraintMatrix > > <https://dealii.org/current/doxygen/deal.II/classConstraintMatrix.html>,constComponentMask > > > > > <https://dealii.org/current/doxygen/deal.II/classComponentMask.html>_mask=ComponentMask > > > > > <https://dealii.org/current/doxygen/deal.II/classComponentMask.html>()) > > | > > > > for this purpose, but because of the discontinuity I get unwanted > effects near > > the discontinuity - that is even when the discontinuity is aligned with > the > > mesh, and even though I am using discontinuous elements (FE_FaceQ), I > get > > unwanted effects because the interpolation point falls exactly between > the > > elements. Is there a way of applying project_boundary_values instead to > fix > > this issue? Basically I cannot seem to find support for component > selection > > in connection with project_boundary_values. > > Can you explain what these unwanted effects are, for example in a picture? > > You're right that the project_b_v() function (currently) doesn't support > selecting individual components. But you can just project onto all > components > into a separate object, and then you only transfer those elements that > correspond to the components you care about. > > Best > Wolfgang > > -- > > Wolfgang Bangerth email: bang...@colostate.edu > > www: http://www.math.colostate.edu/~bangerth/ > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] project_boundary_values for components
Hello, I need to apply a boundary condition (discontinuous function) to a few components of a dof_handler only. I am currently using void VectorTools::interpolate_boundary_values <https://dealii.org/current/doxygen/deal.II/group__constraints.html#ga2376b9a282a2b62be5c55ab62a11a132> (const DoFHandlerType< dim, spacedim > , const types::boundary_id <https://dealii.org/current/doxygen/deal.II/namespacetypes.html#aaf4eb6ec214fa642dfd956f11a9cd2d7> boundary_component, const Function <https://dealii.org/current/doxygen/deal.II/classFunction.html>< spacedim, number > _function, ConstraintMatrix <https://dealii.org/current/doxygen/deal.II/classConstraintMatrix.html> & constraints, const ComponentMask <https://dealii.org/current/doxygen/deal.II/classComponentMask.html> & component_mask=ComponentMask <https://dealii.org/current/doxygen/deal.II/classComponentMask.html>()) for this purpose, but because of the discontinuity I get unwanted effects near the discontinuity - that is even when the discontinuity is aligned with the mesh, and even though I am using discontinuous elements (FE_FaceQ), I get unwanted effects because the interpolation point falls exactly between the elements. Is there a way of applying project_boundary_values instead to fix this issue? Basically I cannot seem to find support for component selection in connection with project_boundary_values. Sincerely, /Anton -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] DataPostprocessor for two DoFHandler s
Hello, I would like to compute a scalar field for outputting purposes. The field is a result of postprocessing of two solution vectors belonging to two different DoFHandler s, but living on the same Triangulation. Is there functionality for doing something like this in the library, or does one have to combine the two vectors into one by creating an FESystem with an associated DoFHandler first? Sincerely, /Anton -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
