Re: [deal.II] Sparsity pattern error in the context of diffuse domain with hp-capability and CIP

[Simon Sticko]

Hi.


1. Is there a way to deal make a correct sparsity pattern accounting for 
hanging nodes, FENothing elements and coupling of faces that have ghost penalty 
or is this really a feature that isn't implemented in deal.ii ?


If you do not care so much about performance at this stage, you can just avoid 
passing the face_has_flux_coupling lambda to make_flux_sparsity_pattern (the 
version that takes cell and face couplings). You will then get a flux coupling 
over all faces. This gives a sparsity pattern with more entries than you need, 
but you will not get an exception.

The alternative is to choose the face_has_flux_coupling lambda more carefully. 
The problem is that you are trying to assign values to entries in the matrix 
that does not exist in the sparsity pattern.



2. Should I assemble the matrix without taking care of the constraints and 
apply them afterwards ?


No, see 1.



3. Why do the sparsity patterns differ from one call to another ?


The implementation of the two functions are different. The simpler function has 
a bug in it. There is an open issue for this:

https://github.com/dealii/dealii/issues/12037

but since there is a workaround, it has not been prioritized.


Best,
Simon


On 01/05/2024 23:39, 'Matthias Nicola Henssler' via deal.II User Group wrote:

Dear all,

*Background :* I am currently trying to implement a diffuse domain formulation 
to solve the poisson equation whilst penalizing the jumps in gradient in very 
similar way to what step-85 showcases. The nuance relative to the mentioned 
tutorial stems the need for AMR close to the boundary.
Having worked through this problem without neglecting the outside of the domain 
(contrary to what is shown in step-85 with the use of FENothing in the 
FECollection object for all the elements for which the level set is positive on 
all vertices), I found myself wanting to implement this feature again.

*Issue :* Taking good care to assigne FENothing only to the cells that extend 
beyond an arbitrary threshold to avoid computing the face values on a cell 
who's neighbor is neglected, I got the following error when computing the flux 
sparsity pattern. Note that hanging nodes are taken into account. Note that 
this issue arise only if cells are flagged for ghost penalty and if elements 
(outside the domain of interest) are attributed the FENothing finite element.



     DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());

     const unsigned int           n_components = fe_collection.n_components();
     Table<2, DoFTools::Coupling> cell_coupling(n_components, n_components);
     Table<2, DoFTools::Coupling> face_coupling(n_components, n_components);
     cell_coupling[0][0] = DoFTools::always;
     face_coupling[0][0] = DoFTools::always;


     const bool                      keep_constrained_dofs = true;
     constraints.clear();
     DoFTools::make_hanging_node_constraints(dof_handler, constraints);
     constraints.close();
     DoFTools::make_flux_sparsity_pattern(dof_handler,
                                           dsp,
                                           constraints,
                                           keep_constrained_dofs,
                                           cell_coupling,
                                           face_coupling,
                                           numbers::invalid_subdomain_id);


-
The violated condition was:
     matrix_values->column() == column
Additional information:
     You are trying to access the matrix entry with index <14,142>, but
     this entry does not exist in the sparsity pattern of this matrix.
-

Having tried the following sister make_flux_sparsity_pattern() function call 
without explicitly passing the cell and face couplings, I get this error:
-
An error occurred in line <4467> of file 
 
in function
     void dealii::AffineConstraints::add_entries_local_to_global(const std::vector&, const 
dealii::AffineConstraints&, const std::vector&, dealii::SparsityPatternBase&, bool, 
const dealii::Table<2, bool>&) const [with number = double]
The violated condition was:
     false
Additional information:
     You are trying to use functionality in deal.II that is currently not
     implemented. In many cases, this indicates that there simply didn't
     appear much of a need for it, or that the author of the original code
     did not have the time to implement a particular case. If you hit this
     exception, it is therefore worth the time to look into the code to
     find out whether you may be able to implement the missing
     functionality. If you do, please consider providing a patch to the
     deal.II development sources (see the deal.II website on how to
     contribute).


Implying that the both 

Re: [deal.II] Sparsity pattern error in the context of diffuse domain with hp-capability and CIP

[Simon Sticko]

Hi.


1. Is there a way to deal make a correct sparsity pattern accounting for 
hanging nodes, FENothing elements and coupling of faces that have ghost penalty 
or is this really a feature that isn't implemented in deal.ii ?


If you do not care so much about performance at this stage, you can just avoid 
passing the face_has_flux_coupling lambda to make_flux_sparsity_pattern (the 
version that takes cell and face couplings). You will then get a flux coupling 
over all faces. This gives a sparsity pattern with more entries than you need, 
but you will not get an exception.

The alternative is to choose the face_has_flux_coupling lambda more carefully. 
The problem is that you are trying to assign values to entries in the matrix 
that does not exist in the sparsity pattern.



2. Should I assemble the matrix without taking care of the constraints and 
apply them afterwards ?


No, see 1.



3. Why do the sparsity patterns differ from one call to another ?


The implementation of the two functions are different. The simpler function has 
a bug in it. There is an open issue for this:

https://github.com/dealii/dealii/issues/12037

but since there is a workaround, it has not been prioritized.


Best,
Simon


On 01/05/2024 23:39, 'Matthias Nicola Henssler' via deal.II User Group wrote:

Dear all,

*Background :* I am currently trying to implement a diffuse domain formulation 
to solve the poisson equation whilst penalizing the jumps in gradient in very 
similar way to what step-85 showcases. The nuance relative to the mentioned 
tutorial stems the need for AMR close to the boundary.
Having worked through this problem without neglecting the outside of the domain 
(contrary to what is shown in step-85 with the use of FENothing in the 
FECollection object for all the elements for which the level set is positive on 
all vertices), I found myself wanting to implement this feature again.

*Issue :* Taking good care to assigne FENothing only to the cells that extend 
beyond an arbitrary threshold to avoid computing the face values on a cell 
who's neighbor is neglected, I got the following error when computing the flux 
sparsity pattern. Note that hanging nodes are taken into account. Note that 
this issue arise only if cells are flagged for ghost penalty and if elements 
(outside the domain of interest) are attributed the FENothing finite element.



     DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());

     const unsigned int           n_components = fe_collection.n_components();
     Table<2, DoFTools::Coupling> cell_coupling(n_components, n_components);
     Table<2, DoFTools::Coupling> face_coupling(n_components, n_components);
     cell_coupling[0][0] = DoFTools::always;
     face_coupling[0][0] = DoFTools::always;


     const bool                      keep_constrained_dofs = true;
     constraints.clear();
     DoFTools::make_hanging_node_constraints(dof_handler, constraints);
     constraints.close();
     DoFTools::make_flux_sparsity_pattern(dof_handler,
                                           dsp,
                                           constraints,
                                           keep_constrained_dofs,
                                           cell_coupling,
                                           face_coupling,
                                           numbers::invalid_subdomain_id);


-
The violated condition was:
     matrix_values->column() == column
Additional information:
     You are trying to access the matrix entry with index <14,142>, but
     this entry does not exist in the sparsity pattern of this matrix.
-

Having tried the following sister make_flux_sparsity_pattern() function call 
without explicitly passing the cell and face couplings, I get this error:
-
An error occurred in line <4467> of file 
 
in function
     void dealii::AffineConstraints::add_entries_local_to_global(const std::vector&, const 
dealii::AffineConstraints&, const std::vector&, dealii::SparsityPatternBase&, bool, 
const dealii::Table<2, bool>&) const [with number = double]
The violated condition was:
     false
Additional information:
     You are trying to use functionality in deal.II that is currently not
     implemented. In many cases, this indicates that there simply didn't
     appear much of a need for it, or that the author of the original code
     did not have the time to implement a particular case. If you hit this
     exception, it is therefore worth the time to look into the code to
     find out whether you may be able to implement the missing
     functionality. If you do, please consider providing a patch to the
     deal.II development sources (see the deal.II website on how to
     contribute).


Implying that the both 

[deal.II] A question on licenses

[Simon Sticko]

Hi.

The function

GridIn::read_msh

uses the Gmsh-api internally to read a mesh from file.  However, Gmsh is
distributed under GPL and not LGPL.

Does this mean that one can not use this function if one wants to use
deal.II under LGPL?

Best,
Simon

-- 
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/CAMhmiLjUDyd0d2MoUc%2BLPk6-LhJ1R9oXuNM9n6VHfTV9CK9Bhw%40mail.gmail.com.

Re: [deal.II] How to compute derivatives of the composite function with the reference cell mapping

[Simon Sticko]

Hi.


If my element is not aligned with the axes, [...]
is this what meant by the "perfect Cartesian" situation?


Yes, exactly. If this is not what you have, I'm afraid the idea of generating 
the quadratures in real space is not an option.


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?


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

Regarding accuracy, you typically get the same order of accuracy. That is, say 
you that you use Lagrange elements of order p in the finite element space for 
you pde solution, and compare using the exact level set function with one that 
has been interpolated (also onto a space of elements of order p), then you 
typically see that

Exact level set function=> error ~ C_1 h^q
Interpolated level set function => error ~ C_2 h^q

(error in some global norm)

However, I know one expert who have claimed that the difference between the 
constants C_1 and C_2 can be surprisingly large for some problems. But I am 
unaware of any publication that has actually compared this.

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?

Best,
Simon

On 04/08/2023 12:39, Anton Evgrafov wrote:

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). 

Re: [deal.II] How to compute derivatives of the composite function with the reference cell mapping

[Simon Sticko]

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/ 

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
 
.


--
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/f1642630-6453-2028-1345-3db743d814f2%40gmail.com.

Re: [deal.II] Query regarding solution of linear system in deal.ii using linear solver

[Simon Sticko]

Hi.

From your attached files it looks like you have saved the matrix and rhs with 
quite few digits. Do you get the same result if you increase the number of 
digits you write to your std::ofstream object?

Compare:
https://en.cppreference.com/w/cpp/io/manip/setprecision

Best,
Simon

On 29/07/2023 16:36, Krish wrote:

Hello,

I am solving a problem that has 30 degrees of freedom. The newton update 
solution obtained from deal.ii (using both Direct and CG solver) seems to be 
less accurate (differs after 4 significant digits)  from the newton update 
solution obtained from Matlab using the same tangent_matrix and system_rhs. I 
am getting the two following solutions from the same tangent_matrix and 
system_rhs (data files containing the tangent_matrix and system_rhs attached).

     Matlab solution deal.ii solution

                  0                   0
   -0.003033422503555  -0.00303356000
                                      0                   0
   -0.002439900952733  -0.00243999000
                       0                   0
   -0.002737307578376  -0.00273743000
    0.001007226333711   0.00100724000
   -0.001069740617723  -0.00106978000
    0.001663965890031   0.00166397000
   -0.002605547185476  -0.00260561000
    0.001113032138953   0.00111304000
   -0.001469379092482  -0.00146943000
    0.000916401048149   0.00091638000
                       0                   0
    0.001674445763014   0.00167444000
                                       0                   0
    0.001041269520982   0.00104126000
   -0.000  -0.000
                   0                   0
    0.003033422503555   0.00303356000
                      0                   0
    0.002439900952733   0.00243999000
                                      0                   0
    0.002737307578376   0.00273743000
    0.001007226333711   0.00100724000
    0.001069740617723   0.00106978000
    0.001663965890031   0.00166397000
    0.002605547185476   0.00260561000
    0.001113032138953   0.00111304000
    0.001469379092482   0.00146943000

I am trying to figure out what could be the reason behind it and how to get a 
more accurate solution.

Thank you in advance.

Krish

--
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/ac75c294-7b25-4351-82e5-dd7ad3fd410bn%40googlegroups.com
 
.


--
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/7fb5ded9-2382-bd5b-a9a2-39df063df956%40gmail.com.

Re: [deal.II] point_value for a vector-valued Problem / FESystem

[Simon Sticko]

Hi.

There are several point_value functions. Right now, you are calling one that 
only works for scalar elements:

https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#a7be5c7eed52308898dfaad91c4cff204

You need to call the more general function:

https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#acd358e9b110ccbf4a7f76796d206b9c7

by passing a vector with each components value as the last argument to the 
function:

Vector value(n_solution_components);
VectorTools::point_value(dof_handler,
 solution,
 point,
 value);

Best,
Simon

On 20/07/2023 12:54, 'Jost Arndt' via deal.II User Group wrote:

Dear everyone,

I have trouble understanding in how to evaluate a solution at a given point 
correctly.
Described in Step-3 for example the function point_value is used to evaluate 
the solution at (0.3, 0.3), and I want to do something quite similar.
However, my solution is vector valued (i.e. I use FESystem ), and I get the 
following error message

 >>An error occurred in line <181> of file 
<./include/deal.II/numerics/vector_tools_point_value.templates.h> in function
     typename VectorType::value_type dealii::VectorTools::point_value(const dealii::Mapping&, const dealii::DoFHandler&, const VectorType&, const 
dealii::Point&) [with int dim =
2; VectorType = dealii::Vector; int spacedim = 2; typename 
VectorType::value_type = double]
The violated condition was:
     dof.get_fe(0).n_components() == 1
Additional information:
     Finite element is not scalar as is necessary for this function


However in the documentation it says "Evaluate a possibly vector-valued finite element 
function defined by the given DoFHandler 
 and nodal vector 
fe_function at the given point point". So I am curious how to evaluate my solution at a 
few fixed points? Did I miss something?

Performance-wise it does not have to be perfect at all, as I want to evaluate 
only ~400 points on a few hundered solutions (time-steps).

Best,

Jost

--
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/7e2a3e13-aef1-4a0e-a357-8dafc3bdf958n%40googlegroups.com
 
.


--
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/229072fe-6f09-1470-35b6-ee34ddc308fd%40gmail.com.

Re: [deal.II] Re: interface between dealii and matlab for data transfer (using mex files?)

[Simon Sticko]

Hi,

I have done the mex-coupling with dealii before. As Bruno said, it is easiest 
to mex with cmake. Example in attached zip. Hope it helps.


Best,
Simon Sticko

On 15/08/2022 15:05, Bruno Turcksin wrote:

Simon

Le lun. 15 août 2022 à 08:52, Simon Wiesheier mailto:simon.wieshe...@gmail.com>> a écrit :


I also thought about transferring the data by reading from/writing to a 
file.
But I was looking for a more general approach without having to create 
files.

Writing files is no less general than coupling directly Matlab and deal.II The 
main advantage of coupling the two codes directly is that it is faster but it 
should only matter if you have to transfer a very large amount of data or if 
you do it very often.

I was also pondering if it is possible to compile the mex file like a 
normal dealii program (by adding a few Matlab macros to the files generated by 
cmake).
But I can not say how much is involved to do this due to my modest 
knowledge about cmake.
If it were possible and someone wanted to guide me through that process, I 
would be willed to try that out; otherwise I go with the easier approach and 
work with files.

The way you probably want to do that is to start from a CMakeLists.txt similar to the one that 
you are using to compile your deal.II code and use 
https://cmake.org/cmake/help/latest/module/FindMatlab.html 
<https://cmake.org/cmake/help/latest/module/FindMatlab.html> to compile the mex file. I 
have never done it so I don't know if it "just works" or if it is more involved.

Best,

Bruno

--
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/CAGVt9eMM_kwCXbSHZONL2p8icdeiwDF%3DQrvQ-qbyQFYxqEPL4g%40mail.gmail.com
 
<https://groups.google.com/d/msgid/dealii/CAGVt9eMM_kwCXbSHZONL2p8icdeiwDF%3DQrvQ-qbyQFYxqEPL4g%40mail.gmail.com?utm_medium=email_source=footer>.


--
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/ad0f6fad-39f8-db0d-d5f2-bd849504dcbb%40gmail.com.
<>

Re: [deal.II] FEInterfaceValues and hp::FECollection

[Simon Sticko]

Hi,

Unfortunately no-one has implemented FEInterfaceValues for hp yet (I think it 
would be great if someone did). The only workaround I know is to NOT use 
FEInterfaceValues and instead:

1. Use 2 hp::FEFaceValues objects, one initialized with cell 0 and one 
initialized with cell 1.

2. Split the jump term in 4 terms: [u][v] = u_0*v_0 - u0*v_1 - u_1*v_0 + u_1*v_1

3. Assemble each term in a separate local matrix (A_00, A_01, A_10, A_11), where

   A_01 - Rows correspond to local dofs on cell 0, columns to local dofs on 
cell 1.

4. Use global_dof_indices of cell 0 and cell 1 to map each local matrix to the 
global matrix.

Best,
Simon



On 09/08/2022 00:58, Marco Feder wrote:

Dear all,

I'm using an hp::FECollection object to describe different FE spaces in my 
grid. In particular, I'm in the following scenario:

x- - - - - -x
|             |
|     0      |
|             |
x - - F - -x
|             |
|     1      |
|             |
x- - - - - -x

On cell number 0 I have a FESystem with FE_Q(1) and FE_Q(1), which corresponds 
to fe_collection[0].
On cell number 1, I have again a FESystem with FE_Q(1) and FE_Nothing, 
corresponding to fe_collection[1].

I need to compute the jump [d_n u_h]_F, being F the face shared by the two 
cells. Unfortunately, I can't find a proper way to initialize FEInterfaceValues 
so that it understands that I need the jump of the first components only. 
Indeed, if I use fe_collection[0] (or fe_collection[1]) as first argument in 
the constructor, I have a runtime error in reinit(cell) saying (as expected):

/    The FiniteElement you provided to FEValues and the FiniteElement that
     belongs to the DoFHandler that provided the cell iterator do not
     match./
/
/
since we have different spaces on cell0 and cell1. Is there a possible 
workaround to this?

Best,
Marco

--
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/2d71b802-132c-4f8a-a20d-03542939036an%40googlegroups.com
 
.


--
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/23e81ebd-cf4e-60ad-ea66-ad825581b533%40gmail.com.

Re: [deal.II] mesh_generator

[Simon Sticko]

Hi,

Did you call

triangulation.execute_coarsening_and_refinement();

after setting the refinement flag?

Best,
Simon

On 23/06/2022 17:22, LY XXXiao wrote:

Hi there,

Has anyone experienced this issue in mesh refinement:

I have already a basic mesh with 120 cells, and I only want to refine one 
element among them, so I set a refine_flag for this element. But the refinement 
does not happen ?? Is there anything else that I missed? I learned step 1 and 
step 6 for meshing problem.

Best regards,
Longying


--
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/0cfddd51-4380-4cd3-81af-9259046b2c8en%40googlegroups.com
 
.


--
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/49efb643-4ca8-f2a6-c9a7-c7c20fa12ee8%40gmail.com.

Re: [deal.II] How to add degrees of freedom to a specified node

[Simon Sticko]

Hi,

Since several XFEM papers formulate this as a Heaviside enrichment around the 
crack, you might find this class useful:

https://www.dealii.org/current/doxygen/deal.II/classFE__Enriched.html

Relating to Wolfgangs answer, note in particular the section "Using enriched and 
non-enriched FEs together".

Best,
Simon


On 18/06/2022 07:43, Wolfgang Bangerth wrote:


Wang Yuan,


  I am using extended finite element method to do crack expansion. I have a 
question, can I increase the degree of freedom for the nodes designated around 
the crack in Deal?


Not easily. deal.II enumerates degrees of freedom based on the finite element 
used on a cell, so if you want to have extra degrees of freedom, you need to 
describe that by using a different element on that cell (using the 
hp-framework, in the same way as step-46 for example).


Meanwhile, another question is can I assign different numbers of Gaussian 
integral points to different elements? Thank you for your help!


Yes! You can use a hp::QCollection together with hp::FEValues. Or you can do it 
by hand: You create multiple quadrature rules and corresponding FEValues 
objects, and then you choose which one is appropriate when you are on a 
concrete cell.

Best
  W.



--
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/56f2ee5e-9bba-fd46-7499-7b925aa3353d%40gmail.com.

Re: [deal.II] Integral over part of domain.

[Simon Sticko]

Hi,

This type of integral can be computed with the tools presented in the Step-85 
tutorial:

https://www.dealii.org/developer/doxygen/deal.II/step_85.html

which you can use if you install the master branch of deal.II (or wait for the 
soon-released 9.4 version).

In particular, take a look at the assemble_system() function in Step-85 and the 
NonMatching::FEValues class:

https://www.dealii.org/developer/doxygen/deal.II/classNonMatching_1_1FEValues.html

Best,
Simon


On 13/06/2022 15:17, Oleg Kmechak wrote:

Hello there,

Facing such a problem.
Want to evaluate integral over circle (see picture below). Integral has such a 
form:
I = Integral[ solution_value(x, y) * my_func(x, y) * dx dy,  over circle].

Currently I am using the simplest way to evaluate such integral - rectangular 
approximation. And most time cost part - is using 
Functions::FEFieldFunction to get solution_value(x, y). Solving FEM is 
faster than evaluation of such integrals.

Also, my_func(x, y) has hidden parameter 'order'. Basically, I am using 
integral to expand field around tip (in the center of circle) into series. And 
also with higher order, accuracy of series parameter (result of integration) is 
poorer.

So maybe once again. How to evaluate fast and accurately integral over such 
circle (not whole domain)?

Also, maybe I can introduce a function which is equal to 1 inside of circle and 
0 - outside. And then integrate over whole domain.
Also, maybe I can adopt The deal.II Library: Integrators (dealii.org) 
 to fit 
my function but not sure how I can do it.

Best regards,
Oleg Kmechak

Inkedinit_river_with_boundary_cond_LI.jpg

--
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/e02ebf7a-c4f2-4202-8cd2-057900260914n%40googlegroups.com
 
.


--
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/79693521-76dc-ece1-e0eb-32fb0bf7c36c%40gmail.com.

Re: [deal.II] Facing problem while installation deal.ii

[Simon Sticko]

Hi,

Just to elaborate slightly on what Chen said: If your installation stops at 52%, it is NOT because 
of "dependencies not satisfied". This is not an error. Can you post the error you get 
from calling make in the terminal? E.g. copy the output you get from cmake and "make" to 
a txt-file and attach that.

Best,
Simon

On 10/06/2022 10:05, 陈敏 wrote:

Dear Abhishek Nath Thakur

These are many examples showing how to use deal.ii, these are corresponding documents on 
the official website . 
But some of these need external packages 

 not bundled in the deal.ii. you need to install them yourself. Finally, you don't need to 
install all of the external packages, only the necessary packages are needed according to 
your applications.

Best
Chen

Abhishek Nath Thakur mailto:abhisheknaththa...@gmail.com>> 于2022年6月10日周五 15:57写道：

Dear deal.ii users,

I tried to install deal.ii 9.2.0 several times but every time i'm getting 
"dependencies not satisfied" (*image attached below*) due to which my 
installation stops at 52% everytime. I'm not sure if its an error or something else. So 
any help is appreciated

Thank you in advance


-- 
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/5483bfa8-a269-4503-8624-557e8aa015edn%40googlegroups.com
 
.

--
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/CADr3OdL74ZJdfrwXhe1M%2BB0sZ4%2BwtNBBhV4bGUCqv1vspUvEqg%40mail.gmail.com
 
.


--
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/1bb4cf10-c237-938f-666a-7b87e175d13d%40gmail.com.

Re: [deal.II] Encountering error in installation of deal.ii

[Simon Sticko]

Hi,

From the error message, it looks like you don't have C++ compiler installed.
If you are using ubuntu, try

sudo apt-get install g++

and rerun cmake.

Best,
Simon

On 08/06/2022 16:32, Abhishek Nath Thakur wrote:

Respected deal.ii users,

I'm trying to install deal.ii 9.2.0 version but keep getting the error(*error 
image attached below*). I'm very new to this so i can really use someone's help 
to find me the way out of this.

Thank you very much in advance


--
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/e5dd7aeb-ade5-4c4a-9774-6ce10535b3b3n%40googlegroups.com
 
.


--
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/6e38f5a1-4032-47c5-e31c-0036353adc35%40gmail.com.

Re: [deal.II] How to use create_mass_matrix tools with FESystem

[Simon Sticko]

Hi,

In the stacktrace it looks like the error does not come from create_mass_matrix 
but from VectorTools::project. It complains that the passed Function has 1 
component when it should have dim components.

If you are projecting your own function, you need to make sure you set the 
number of components in the constructor, e.g. as for the body-force in step-18:

https://github.com/dealii/dealii/blob/39e29ba7e1cb8c7f711df6b89d19a76394b16f6a/examples/step-18/step-18.cc#L552

Best,
Simon



On 10/04/2022 01:51, Matthew Rich wrote:

Hi all,

I am trying to use the helper tools to assemble a system for a vector valued 
problem.

I am looking to have 2 or 3 displacement DoFs.

So I use

FESystem fe;

and initialize with

fe(FE_Q 
(1), dim)

So my code compiles, but when I try and run I get the following


An error occurred in line <182> of file 

 in function
     void dealii::VectorTools::internal::project_matrix_free(const dealii::Mapping&, const dealii::DoFHandler&, 
const dealii::AffineConstraints&, const dealii::Quadrature&, const dealii::Function::value_type>&, dealii::LinearAlgebra::distributed::Vector&, bool, const 
dealii::Quadrature<(dim - 1)>&, bool) [with int components = 2; int fe_degree = 1; int dim = 2; Number = double; int spacedim = 2; typename 
dealii::LinearAlgebra::distributed::Vector::value_type = double]
The violated condition was:
     dof.get_fe(0).n_components() == function.n_components
Additional information:
     Dimension 2 not equal to 1.

Stacktrace:
---
#0  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::internal::project_matrix_free<2, 1, 2, double, 2>(dealii::Mapping<2, 
2> const&, dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints const&, dealii::Quadrature<2> const&, 
dealii::Function<2, dealii::LinearAlgebra::distributed::Vector::value_type> const&, 
dealii::LinearAlgebra::distributed::Vector&, bool, dealii::Quadrature<(2)-(1)> const&, bool)
#1  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::internal::project_matrix_free_degree<2, 2, double, 2>(dealii::Mapping<2, 
2> const&, dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints const&, dealii::Quadrature<2> const&, 
dealii::Function<2, dealii::LinearAlgebra::distributed::Vector::value_type> const&, 
dealii::LinearAlgebra::distributed::Vector&, bool, dealii::Quadrature<(2)-(1)> const&, bool)
#2  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::internal::project_matrix_free_component<2, double, 2>(dealii::Mapping<2, 
2> const&, dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints const&, dealii::Quadrature<2> const&, 
dealii::Function<2, dealii::LinearAlgebra::distributed::Vector::value_type> const&, 
dealii::LinearAlgebra::distributed::Vector&, bool, dealii::Quadrature<(2)-(1)> const&, bool)
#3  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::internal::project_matrix_free_copy_vector<2, dealii::Vector, 
2>(dealii::Mapping<2, 2> const&, dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints::value_type> 
const&, dealii::Quadrature<2> const&, dealii::Function<2, dealii::Vector::value_type> const&, dealii::Vector&, bool, 
dealii::Quadrature<(2)-(1)> const&, bool)
#4  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::internal::project, 2>(dealii::Mapping<2, 2> 
const&, dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints::value_type> const&, dealii::Quadrature<2> 
const&, dealii::Function<2, dealii::Vector::value_type> const&, dealii::Vector&, bool, dealii::Quadrature<(2)-(1)> 
const&, bool)
#5  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::project<2, dealii::Vector, 2>(dealii::Mapping<2, 2> const&, 
dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints::value_type> const&, dealii::Quadrature<2> const&, 
dealii::Function<2, dealii::Vector::value_type> const&, dealii::Vector&, bool, dealii::Quadrature<(2)-(1)> const&, bool)
#6  /usr/lib/x86_64-linux-gnu/libdeal.ii.g.so.9.2.0: void dealii::VectorTools::project<2, dealii::Vector, 2>(dealii::DoFHandler<2, 2> 
const&, dealii::AffineConstraints::value_type> const&, dealii::Quadrature<2> const&, dealii::Function<2, 
dealii::Vector::value_type> const&, dealii::Vector&, bool, dealii::Quadrature<(2)-(1)> const&, bool)
#7  ./final_proj: final_proj::WaveEquation<2>::run()
#8  ./final_proj: main


Using some crude debugging it is complaining about the create_mass_matrix call. 
Now I know there are two variants. I was trying to create the matrix without 
passing a mapping but that is generating this error. My issue is I am not sure 
how to create a mapping that works.

How does one make a mapping for a system? I looked at some of the tutorials but 
I could not find what I am looking for


Thanks in advance.

[deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder

[Simon Sticko]

Hi,

> and the code seems to be working fine! Is there anything that I am 
> missing in your opinion?

No. If you're happy, I'm happy.

Best,
Simon

On Wednesday, July 1, 2020 at 7:51:11 AM UTC+2, Samuel Rodriguez wrote:
>
> Hi Simon,
>
> I am a slow learner but you were correct. Here is the first image of the 
> function plotted along the radius. As you can see there is only noise but 
> value of the solution is constant. I can also change the value of the 
> solution by changing the value of the boundary condition as you can also 
> see below:
>
> [image: Screen Shot 2020-06-30 at 2.20.11 PM.png]
>
> [image: Screen Shot 2020-06-30 at 10.22.07 PM.png]
>
> I was able to change my boundary condition to something ridiculous like:
>
> class NeumannBoundaryCondition : public Function<3>
> {
> public:
>   double
>   value(const Point<3> ,  const unsigned int component = 0) const
>   {
> double test = 1;
> for(int i = 0; i < 3; i++)
>   test += 20 + std::pow(point(i), 20) + std::pow(point(i), 10) + 
> exp(point(i));
> std::cout << test << std::endl;
> return test;
>   }
> };
>
>
>
> The reason this was not initially working was because when I implemented 
> the boundary condition below:
>
>   for(unsigned int face_number = 0; face_number < 
> GeometryInfo<3>::faces_per_cell; ++face_number)
>   {
>if(cell->face(face_number)->at_boundary() && 
> (cell->face(face_number)->boundary_id() == 0))
> {
>  fe_face_values.reinit(cell, face_number);
>  for (unsigned int q_point = 0; q_point < n_face_q_points; 
> ++q_point)
>   {
>double neumannvalue = 
> boundary_condition.value(fe_face_values.quadrature_point(q_point));
>for (unsigned int i = 0; i < dofs_per_cell; ++i)
> {
>  const unsigned int component_i = 
> fe.system_to_component_index(i).first;
>  cell_rhs(i) += (neumannvalue *
>  fe_face_values.shape_value(i, q_point) *
>  fe_face_values.JxW(q_point));
> }
>   }
> 
>}
>  }
>
>
> I never added it the right hand side (I think that I am using that 
> terminology correctly). I added the previous lines of code right above these
>
>   cell->get_dof_indices(local_dof_indices);
>
>   for (unsigned int i = 0; i < dofs_per_cell; ++i)
> for (unsigned int j = 0; j < dofs_per_cell; ++j)
>   system_matrix.add(local_dof_indices[i],
> local_dof_indices[j],
> cell_matrix(i, j));
>   for (unsigned int i = 0; i < dofs_per_cell; ++i)
> system_rhs(local_dof_indices[i]) += cell_rhs(i);
>
>
> and the code seems to be working fine! Is there anything that I am missing 
> in your opinion?
>
> Best,
>
> Samuel Rodriguez
>
>
>
>
>
>
>
>
>
>
>
> On Tuesday, June 30, 2020 at 3:33:25 AM UTC-7, Simon Sticko wrote:
>>
>> Hi,
>> From the figure you can't really tell whether the boundary condition is 
>> fulfilled or not. If you want to check if the normal derivative is 
>> reasonable (almost equal to 0) at the boundary, the easiest is probably to 
>> plot the solution on a line in the radial direction. This can be done by 
>> using
>> Filters -> Data Analysis -> Plot Over Line in Paraview.
>>
>> I don't quite understand what you mean when you say that you have tried 
>> changing the value of component. Component is intended to be used for 
>> vector valued problems. If you are solving Poisson, it doesn't make sense 
>> to set it to something different than 0. 
>>
>> Note that, right now you are calling the function 
>> NeumannBoundaryCondition::value with a single argument (of type Point<3>). 
>> This means that you get the default argument for component: component=0.
>>
>> Best,
>> Simon
>>
>> On Tuesday, June 30, 2020 at 2:59:20 AM UTC+2, Samuel Rodriguez wrote:
>>>
>>> Hi Simon,
>>>
>>> It seems I might not have understood my own question. You answered the 
>>> question I was trying to ask. Would you be able to help me out with one 
>>> more thing? I followed your advice but the boundary conditions do not seem 
>>> bedo not seem to be getting the correct boundary conditions, so I don't 
>>> think they are being applied correctly. Here is the simple piece of code 
>>> that I wrote:
>>>
>>> class NeumannBoundaryCondition : public Function<

[deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder

[Simon Sticko]

Hi,
>From the figure you can't really tell whether the boundary condition is 
fulfilled or not. If you want to check if the normal derivative is 
reasonable (almost equal to 0) at the boundary, the easiest is probably to 
plot the solution on a line in the radial direction. This can be done by 
using
Filters -> Data Analysis -> Plot Over Line in Paraview.

I don't quite understand what you mean when you say that you have tried 
changing the value of component. Component is intended to be used for 
vector valued problems. If you are solving Poisson, it doesn't make sense 
to set it to something different than 0. 

Note that, right now you are calling the function 
NeumannBoundaryCondition::value with a single argument (of type Point<3>). 
This means that you get the default argument for component: component=0.

Best,
Simon

On Tuesday, June 30, 2020 at 2:59:20 AM UTC+2, Samuel Rodriguez wrote:
>
> Hi Simon,
>
> It seems I might not have understood my own question. You answered the 
> question I was trying to ask. Would you be able to help me out with one 
> more thing? I followed your advice but the boundary conditions do not seem 
> bedo not seem to be getting the correct boundary conditions, so I don't 
> think they are being applied correctly. Here is the simple piece of code 
> that I wrote:
>
> class NeumannBoundaryCondition : public Function<3>
> {
> public:
>   double
>   value(const Point<3> ,  const unsigned int component = 0) const
>   {
> return component;
>   }
> };
>
>
> This neumann boundary condition is only meant to apply a constant value of 
> zero on the sides. Here is how I implemented it:
>
> QGauss<3> quadrature_formula(fe.degree + 1);
> QGauss<3 - 1> face_quadrature_formula(fe.degree + 1);
> FEFaceValues<3> fe_face_values(fe,
>face_quadrature_formula,
>update_values | update_quadrature_points |
>update_normal_vectors |
>update_JxW_values);const NeumannBoundaryCondition 
> boundary_condition;
>
> const unsigned int n_face_q_points = face_quadrature_formula.size();
>
>
> for(unsigned int face_number = 0; face_number < GeometryInfo<3>::
> faces_per_cell; ++face_number)
> {
>   if(cell->face(face_number)->at_boundary() && (cell->face(face_number)->
> boundary_id() == 0))
>   {
>fe_face_values.reinit(cell, face_number);
>for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point) {
>  const double neumannvalue = boundary_condition.value(fe_face_values.
> quadrature_point(q_point));
>  for (unsigned int i = 0; i < dofs_per_cell; ++i) {
>const unsigned int component_i = fe.system_to_component_index(i).
> first;
>cell_rhs(i) += (neumannvalue *
>fe_face_values.shape_value(i, q_point) *
>fe_face_values.JxW(q_point));
>}
>  }
>  
>   }
> }
>
>
> This seems correct to me, but the code is showing a cyldiner that looks 
> like this:
>
>
> [image: Screen Shot 2020-06-29 at 5.50.19 PM.png]
> which doesn't appear to be the right solution. Even when I change the 
> value of component to something ridiculous such as 10,000,000 the solution 
> on the cylinder does not change. Would you be able to tell me what I am 
> doing wrong in the implementation?
>
> - Samuel
>
>
>
> On Friday, June 26, 2020 at 2:21:19 AM UTC-7, Simon Sticko wrote:
>>
>> Hi,
>> I'm not sure I understand your question correctly, but if you want to 
>> impose an inhomogeneous Neumann boundary condition you would typically 
>> write a small class representing the boundary condition you want:
>>
>> class NeumannBoundaryCondition : public Function<3>
>> {
>> public:
>>   double
>>   value(const Point<3> , const unsigned int component = 0) const
>>   {
>> // Implement the boundary condition you want here.
>> return 1;
>>   }
>> };
>>
>> Then just use this when you assemble:
>>
>> const NeumannBoundaryCondition boundary_condition;
>>
>> const Point<3>  = fe_face_values.quadrature_point(q_point);
>>
>> const double neumann_value = boundary_condition.value(point);
>>
>> cell_rhs(i) += (neumann_value *
>> fe_face_values.shape_value(i, q_point) *
>> fe_face_values.JxW(q_point));
>>
>>
>> Does that answer your question?
>>
>> Best,
>> Simon
>>
>> On Friday, June 26, 2020 at 6:02:48 AM UTC+2, Samuel Rodriguez wrote:
>>>
>>> Hello Everyone,
>>>
>>> I am a masters student barely getting into working with the very heavy 
>>> theoreti

[deal.II] Re: Non Homogenous Boundary Conditions on a Cylinder

[Simon Sticko]

Hi,
I'm not sure I understand your question correctly, but if you want to 
impose an inhomogeneous Neumann boundary condition you would typically 
write a small class representing the boundary condition you want:

class NeumannBoundaryCondition : public Function<3>
{
public:
  double
  value(const Point<3> , const unsigned int component = 0) const
  {
// Implement the boundary condition you want here.
return 1;
  }
};

Then just use this when you assemble:

const NeumannBoundaryCondition boundary_condition;

const Point<3>  = fe_face_values.quadrature_point(q_point);

const double neumann_value = boundary_condition.value(point);

cell_rhs(i) += (neumann_value *
fe_face_values.shape_value(i, q_point) *
fe_face_values.JxW(q_point));


Does that answer your question?

Best,
Simon

On Friday, June 26, 2020 at 6:02:48 AM UTC+2, Samuel Rodriguez wrote:
>
> Hello Everyone,
>
> I am a masters student barely getting into working with the very heavy 
> theoretical foundations of using deal.ii. I have been working on a research 
> project that uses a code that creates a cylinder with Dirichlet boundary 
> conditions on the hull/bottom faces and a Neumann boundary condition on the 
> top face Recently I have been trying to change the Dirichlet boundary 
> condition on the hull of the cylinder to a Neumann boundary condition. Our 
> code can be found here inside of QuasistaticBrownianThermalNoise.cpp:
>
> Numerical Coating Thermal Noise 
> 
>
> I have applied the boundary condition but it does not seem to be working. 
> So what I did instead was go to the examples in the tutorial section to 
> learn ow to apply non-homogenous boundary conditions on a cylinder. I have 
> solved the laplacian on the this cylinder using step 3. I have tried to use 
> step 7 to implement non-homogeneous boundary conditions but step 7 requires 
> using a known solution. I would like to be able to set the derivative of 
> the function on the hull to zero and the derivative on the top face to some 
> known function. A Dirichlet boundary condition would be set on the bottom 
> face. The assemble function is down below. My question is how exactly do I 
> apply the Neumann boundary condition to the faces of the cylinder. In the 
> code below I have comments that explain my question in more detail. My 
> question begins where I loop over the faces of the cylinder. However, I 
> included the entire assembly function for completeness in case it was 
> needed. I also attached the whole cpp file in case that was also needed. 
> Any help would be greatly appreciated. If I can provide any more 
> information I would be happy to do so.
>
> QGauss<3> quadrature_formula(fe.degree + 1);
> QGauss<3 - 1> face_quadrature_formula(fe.degree + 1);
>
> FEValues<3> fe_values(fe,
> quadrature_formula,
> update_values | update_gradients | update_JxW_values);
> FEFaceValues<3> fe_face_values(fe,
> face_quadrature_formula,
> update_values | update_quadrature_points |
> update_normal_vectors |
> update_JxW_values);
>
> const unsigned int dofs_per_cell = fe.dofs_per_cell;
> const unsigned int n_q_points = quadrature_formula.size();
> const unsigned int n_face_q_points = face_quadrature_formula.size();
>
> FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
> Vector cell_rhs(dofs_per_cell);
>
> std::vector local_dof_indices(dofs_per_cell);
>
> for (const auto  : dof_handler.active_cell_iterators())
> {
> fe_values.reinit(cell);
>
> cell_matrix = 0;
> cell_rhs = 0;
>
> for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
> {
> for (unsigned int i = 0; i < dofs_per_cell; ++i)
> for (unsigned int j = 0; j < dofs_per_cell; ++j)
> cell_matrix(i, j) +=
> (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
> fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
> fe_values.JxW(q_index)); // dx
> for (unsigned int i = 0; i < dofs_per_cell; ++i)
> cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
> 1 * // f(x_q)
> fe_values.JxW(q_index)); // dx
> }
> cell->get_dof_indices(local_dof_indices);
>
> for (unsigned int i = 0; i < dofs_per_cell; ++i)
> for (unsigned int j = 0; j < dofs_per_cell; ++j)
> system_matrix.add(local_dof_indices[i],
> local_dof_indices[j],
> cell_matrix(i, j));
> for (unsigned int i = 0; i < dofs_per_cell; ++i)
> system_rhs(local_dof_indices[i]) += cell_rhs(i);
>
> // I know that the hull of the cylinder has boundary id 0, the top face 
> has boundary id 1, and the bottom of the face has boundary id 2.
> for(unsigned int face_number = 0; face_number < 
> GeometryInfo<3>::faces_per_cell; 
> ++face_number)
> {
> if(cell->face(face_number)->at_boundary() && 
> (cell->face(face_number)->boundary_id() == 0))
> {
> fe_face_values.reinit(cell, face_number);
> // If we come in here we have found a face that belongs to the boundary 
> condtion of the hull
> // I know that I am supposed to do something like the code in green below, 

Re: [deal.II] How to transform points on faces from unit cells to real cells?

[Simon Sticko]

I just want to clarify that Dougs suggestion also works,
if we use the QProjector class to lift the points
on the face from R^(dim-1) to R^dim before we transform them:

// Create a dummy quadrature from support points on face.
const Quadrature dummy_face_quadrature(
fe.get_unit_face_support_points());

// Lift support points on face from R^(dim-1) to R^dim
const Quadrature dummy_quadrature =
QProjector::project_to_face(dummy_face_quadrature, face_no);

// Now we can use transform_unit_to_real_cell
const Point real_point = mapping.transform_unit_to_real_cell(
cell, dummy_quadrature.point(quadrature_points_no));

Depending on how your code looks, this might be an easier solution.

Best,
Simon


On Monday, June 22, 2020 at 11:17:35 PM UTC+2, Lex Lee wrote:
>
> Doug, thanks for sharing me your ideas and the discussion. -Lex
>
>
> On Monday, June 22, 2020 at 2:13:25 PM UTC-7, Doug Shi-Dong wrote:
>>
>> Agreed. Sorry for the confusion. Looked back at my code and it's probably 
>> the most straightforward way.
>>
>> For some other unrelated reason I didn't want to use FEValues and 
>> directly deal with the Mapping class.
>>
>> Doug
>>
>> On Monday, June 22, 2020 at 5:05:57 PM UTC-4, Simon Sticko wrote:
>>>
>>> Hi,
>>>
>>> Lex, given your previous question:
>>>
>>> https://groups.google.com/forum/#!topic/dealii/xghVE7Km78o
>>>
>>> If you are already using an FEFaceValues object with a dummy-quadrature 
>>> (created from FiniteElement::get_unit_face_support_points())
>>> then you can use one of the following functions on FEFaceValues to get 
>>> the location of the support points in real space:
>>>
>>> quadrature_point()
>>> get_quadrature_points()
>>>
>>> Best,
>>> Simon
>>>
>>> On Monday, June 22, 2020 at 10:53:10 PM UTC+2, Lex Lee wrote:
>>>>
>>>> Hello Doug,
>>>>
>>>>
>>>> Thanks for your help.
>>>>
>>>> However, I need to say, I just had done exactly what you described 
>>>> before I posted this question.
>>>>
>>>> "FiniteElement::get_unit_face_support_points()” returns Points;
>>>>
>>>> “transform_unit_to_real_cell()” requires Point as input.
>>>>
>>>> The dimensions of points in cells is different from that on faces. 
>>>> That’s one of the reasons why I posted this question.
>>>>
>>>>
>>>> If I misunderstand you, please kindly let me know.
>>>>
>>>>
>>>>
>>>> Regards,
>>>> Lex 
>>>>
>>>>
>>>> On Jun 21, 2020, at 7:05 PM, Doug Shi-Dong  wrote:
>>>>
>>>> Hello Lex,
>>>>
>>>> transform_unit_to_real_cell() takes a point as unit cell's point an 
>>>> input and returns the physical point location. Therefore, you can give the 
>>>> unit face support points from the 
>>>> FiniteElement::get_unit_face_support_points() (assuming your FiniteElement 
>>>> has support points). And pass it to transform_unit_to_real_cell().
>>>>
>>>> Doug
>>>>
>>>> On Sunday, June 21, 2020 at 8:10:30 PM UTC-4, Lex Lee wrote:
>>>>>
>>>>> Hello Doug,
>>>>>
>>>>> Thanks for your kind help.
>>>>>
>>>>> I am trying to understand your suggestion in the original email.
>>>>>
>>>>> Are you suggesting me using the member function 
>>>>> "transform_unit_to_real_cell()" to get all support points in a unit cell? 
>>>>> In the step, sort out the face points I need? Am I correct?
>>>>>
>>>>>
>>>>> Regards,
>>>>>
>>>>> Lex
>>>>>
>>>>> On Sunday, June 21, 2020 at 5:00:35 PM UTC-7, Doug Shi-Dong wrote:
>>>>>>
>>>>>> Hello Lex,
>>>>>>
>>>>>> You were close.
>>>>>>
>>>>>>
>>>>>> https://www.dealii.org/current/doxygen/deal.II/classMapping.html#ae5df63553eb8ed170c3b90524853dd48
>>>>>>
>>>>>> The project_real_point_to_unit_point_on_face() is useful is you take 
>>>>>> a "volume" point within the cell. Since you are mapping from unit to 
>>>>>> real, 
>>>>>> you would always know whether your point is on the face or not.
>>>>>>
>>>>>> Therefore, tr

Re: [deal.II] How to transform points on faces from unit cells to real cells?

[Simon Sticko]

Hi,

Lex, given your previous question:

https://groups.google.com/forum/#!topic/dealii/xghVE7Km78o

If you are already using an FEFaceValues object with a dummy-quadrature 
(created from FiniteElement::get_unit_face_support_points())
then you can use one of the following functions on FEFaceValues to get the 
location of the support points in real space:

quadrature_point()
get_quadrature_points()

Best,
Simon

On Monday, June 22, 2020 at 10:53:10 PM UTC+2, Lex Lee wrote:
>
> Hello Doug,
>
>
> Thanks for your help.
>
> However, I need to say, I just had done exactly what you described before 
> I posted this question.
>
> "FiniteElement::get_unit_face_support_points()” returns Points;
>
> “transform_unit_to_real_cell()” requires Point as input.
>
> The dimensions of points in cells is different from that on faces. That’s 
> one of the reasons why I posted this question.
>
>
> If I misunderstand you, please kindly let me know.
>
>
>
> Regards,
> Lex 
>
>
> On Jun 21, 2020, at 7:05 PM, Doug Shi-Dong  > wrote:
>
> Hello Lex,
>
> transform_unit_to_real_cell() takes a point as unit cell's point an input 
> and returns the physical point location. Therefore, you can give the unit 
> face support points from the FiniteElement::get_unit_face_support_points() 
> (assuming your FiniteElement has support points). And pass it to 
> transform_unit_to_real_cell().
>
> Doug
>
> On Sunday, June 21, 2020 at 8:10:30 PM UTC-4, Lex Lee wrote:
>>
>> Hello Doug,
>>
>> Thanks for your kind help.
>>
>> I am trying to understand your suggestion in the original email.
>>
>> Are you suggesting me using the member function 
>> "transform_unit_to_real_cell()" to get all support points in a unit cell? 
>> In the step, sort out the face points I need? Am I correct?
>>
>>
>> Regards,
>>
>> Lex
>>
>> On Sunday, June 21, 2020 at 5:00:35 PM UTC-7, Doug Shi-Dong wrote:
>>>
>>> Hello Lex,
>>>
>>> You were close.
>>>
>>>
>>> https://www.dealii.org/current/doxygen/deal.II/classMapping.html#ae5df63553eb8ed170c3b90524853dd48
>>>
>>> The project_real_point_to_unit_point_on_face() is useful is you take a 
>>> "volume" point within the cell. Since you are mapping from unit to real, 
>>> you would always know whether your point is on the face or not.
>>>
>>> Therefore, transform_unit_to_real_cell() should do the trick. Unless you 
>>> somehow actually need to project a "volume" point, but projecting within 
>>> unit cell is easy, which you would then follow with 
>>> transform_unit_to_real_cell().
>>>
>>> Best regards,
>>>
>>> Doug
>>>
>>> On Saturday, June 20, 2020 at 8:16:22 PM UTC-4, Lex Lee wrote:

 Hello Deal.II Users,


 I want to get the physical (geometry) coordinates for support points on 
 cell faces. Also I know that there are such member functions that can me 
 map points between reference and real cells in this link: 
 https://www.dealii.org/current/doxygen/deal.II/classMapping.html . 

 However, via the link mentioned above, I only find a member function 
 called project_real_point_to_unit_point_on_face 
 
  (real 
 to unit) , no functions like project_unit_point_to_real_point_on_face 
 (unit to real).  

 Could someone kindly tell me how I can solve this problem?  


 Regards,

 Lex

>>>
> -- 
> 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/uglzE4d4Flo/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to 
> dea...@googlegroups.com .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/fd8f078d-f0d4-4319-8977-a9a8e4672fb5o%40googlegroups.com
>  
> 
> .
>
>
>

-- 
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/aedf08b6-71cf-4452-9619-2155a497f35bo%40googlegroups.com.

[deal.II] Re: Error calculating shape gradient

[Simon Sticko]

Hi,
Yes. Here, you add a number of scalar elements of different degree:

> for (unsigned int deg=1; deg <= max_fe_degree; ++deg)
> {
> fe_collection.push_back(FE_Q(deg));
> quadrature_collection.push_back(QGauss(deg+1));
> face_quadrature_collection.push_back(QGauss(deg+1));
> }

which makes sense for a scalar problem, like the Poisson equation.

But if you solve for displacement in a material the solution has dim 
components:

u = (u_x, u_y, u_z)

so you want to use a FESystem element with dim components:

fe_collection.push_back(FESystem(FE_Q(deg), dim));

Best,
Simon


On Wednesday, June 10, 2020 at 12:03:57 PM UTC+2, A.Z Ihsan wrote:
>
> Hi Simon, 
>
> i followed the tutorial step-27 about hp-fem, using the similar 
> constructor with dim = 3
> for (unsigned int deg=1; deg <= max_fe_degree; ++deg)
>   {
> fe_collection.push_back(FE_Q(deg));
> quadrature_collection.push_back(QGauss(deg+1));
> face_quadrature_collection.push_back(QGauss(deg+1));
>   }
>
> I tried to look into the fe_collection.n_components(), it yields 1.
> is that what you meant?
>
> BR, 
> Ihsan
>
> On Wednesday, June 10, 2020 at 9:07:58 AM UTC+2, Simon Sticko wrote:
>>
>> Hi,
>> from the error message you can see that the element you are using only 
>> has 1 component. You get an error because you are trying to access 
>> component 1, which doesn't exist. Since your element should have dim 
>> components, there is likely something wrong with how you create your 
>> element. It should probably be similar as in step-8.
>>
>> Best,
>> Simon
>>
>>
>> On Wednesday, June 10, 2020 at 8:29:03 AM UTC+2, A.Z Ihsan wrote:
>>>
>>> Hi all, 
>>>
>>> i am getting this error while calculating the local cell matrix for my 
>>> hp fem application
>>>
>>> dealii::Tensor<1, spacedim> dealii::FEValuesBase>>>> spacedim>::shape_grad_component(unsigned int, unsigned int, unsigned int) 
>>>>> const [with int dim = 3; int spacedim = 3]
>>>>
>>>> The violated condition was: 
>>>>
>>>> component < fe->n_components()
>>>>
>>>> Additional information: 
>>>>
>>>> Index 1 is not in the half-open range [0,1).
>>>>
>>>> Stacktrace:
>>>>
>>>> ---
>>>>
>>>> #0  ../build/local: dealii::FEValuesBase<3, 
>>>>> 3>::shape_grad_component(unsigned int, unsigned int, unsigned int) const
>>>>
>>>> #1  ../build/local: dealii::SymmetricTensor<2, 3, double> 
>>>>> LinearElasticity::get_strain<3>(dealii::FEValues<3, 3> const&, unsigned 
>>>>> int, unsigned int)
>>>>
>>>> #2  ../build/local: LinearElasticity::HPSolver<3, 
>>>>> dealii::Vector 
>>>>> >::assemble_cell_matrix(dealii::TriaActiveIterator>>>> > 
>>>>> 3>, false> > const&, dealii::FullMatrix&, dealii::hp::FEValues<3, 
>>>>> 3>&)
>>>>
>>>> #3  ../build/local: 
>>>>> LinearElasticity::HPSerialSolver<3>::assemble_linear_system(LinearElasticity::HPSerialSolver<3>::LinearSystem&)
>>>>
>>>> #4  ../build/local: LinearElasticity::HPSerialSolver<3>::solve_problem()
>>>>
>>>>  
>>>>
>>>>
>>> The idea is similar to the tutorial step-27, but here i use the 
>>> symmetric tensor for returning strain. Can someone explain to me what the 
>>> error is and how to resolve it? 
>>>
>>> BR, 
>>> Ihsan
>>>
>>

-- 
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/dac7f402-fb90-4f1c-940d-5450ddd0cf5do%40googlegroups.com.

[deal.II] Re: Error calculating shape gradient

[Simon Sticko]

Hi,
from the error message you can see that the element you are using only has 
1 component. You get an error because you are trying to access component 1, 
which doesn't exist. Since your element should have dim components, there 
is likely something wrong with how you create your element. It should 
probably be similar as in step-8.

Best,
Simon


On Wednesday, June 10, 2020 at 8:29:03 AM UTC+2, A.Z Ihsan wrote:
>
> Hi all, 
>
> i am getting this error while calculating the local cell matrix for my hp 
> fem application
>
> dealii::Tensor<1, spacedim> dealii::FEValuesBase>> spacedim>::shape_grad_component(unsigned int, unsigned int, unsigned int) 
>>> const [with int dim = 3; int spacedim = 3]
>>
>> The violated condition was: 
>>
>> component < fe->n_components()
>>
>> Additional information: 
>>
>> Index 1 is not in the half-open range [0,1).
>>
>> Stacktrace:
>>
>> ---
>>
>> #0  ../build/local: dealii::FEValuesBase<3, 
>>> 3>::shape_grad_component(unsigned int, unsigned int, unsigned int) const
>>
>> #1  ../build/local: dealii::SymmetricTensor<2, 3, double> 
>>> LinearElasticity::get_strain<3>(dealii::FEValues<3, 3> const&, unsigned 
>>> int, unsigned int)
>>
>> #2  ../build/local: LinearElasticity::HPSolver<3, dealii::Vector 
>>> >::assemble_cell_matrix(dealii::TriaActiveIterator>> > 
>>> 3>, false> > const&, dealii::FullMatrix&, dealii::hp::FEValues<3, 
>>> 3>&)
>>
>> #3  ../build/local: 
>>> LinearElasticity::HPSerialSolver<3>::assemble_linear_system(LinearElasticity::HPSerialSolver<3>::LinearSystem&)
>>
>> #4  ../build/local: LinearElasticity::HPSerialSolver<3>::solve_problem()
>>
>>  
>>
>>
> The idea is similar to the tutorial step-27, but here i use the symmetric 
> tensor for returning strain. Can someone explain to me what the error is 
> and how to resolve it? 
>
> BR, 
> Ihsan
>

-- 
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/9e6af2f7-695c-4a4f-bcd6-17c9a4f008cao%40googlegroups.com.

[deal.II] Re: Resolving with different RHS

[Simon Sticko]

Hi,
a better option than computing the inverse is to factorize the matrix. This 
can be done using the SparseDirectUMFPACK solver:

https://www.dealii.org/current/doxygen/deal.II/classSparseDirectUMFPACK.html

You might want to take a look at step 22, which uses this solver:

https://www.dealii.org/current/doxygen/deal.II/step_22.html

Best,
Simon

On Tuesday, June 2, 2020 at 12:04:38 PM UTC+2, Andreas Kyritsakis wrote:
>
> Dear all,
>
> I have a problem where a Poisson-like equation has to be solved again and 
> again in many time steps. At each timestep the LHS remains the same while 
> the RHS changes (slightly). My current implementation is to use the 
> standard SolverCG, passing the old solution as initial, which already 
> reduces the number of CG steps required. Yet, I wonder whether my approach 
> is naive and there is a faster way to implement this, taking better 
> advantage of the fact that the LHS stays always the same. I initially 
> thought about calculating the inverse matrix once and just doing matrix 
> multiplication, but since the mass matrix is a huge sparse matrix, only 
> storing the inverse (which in general is not sparse) would require huge 
> memory. I also thought about the LinearOperator concepts, but if I 
> understood correctly, they just implement a nice wrapper to call a solver 
> each time. Am I missing something?
>
> Cheers,
> Andreas
>

-- 
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/4bc6a09b-834c-485e-a26c-cbbbc6a45338%40googlegroups.com.

Re: [deal.II] Re: FE_Q

[Simon Sticko]

Hi,
Just to elaborate on what Daniel already said, Oliver, note that it is the 
constructor of the Step3-class that calls the constructor of FE_Q, 
that is

Step3::Step3()
  : fe(1)
  , dof_handler(triangulation)
{}

If you miss this call in the Step3-constructor, you will get the error you 
stated earlier.

Best,
Simon

-- 
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/633d2d99-c5fb-4720-94f4-38d4127d8b06%40googlegroups.com.

Re: [deal.II] Re: Bug Report in GridTools::get_active_neighbors()

[Simon Sticko]

Hi,
Fixed in #9108, already merged in. I should probably have written that here for 
completeness. Sorry.

Best,
Simon

-- 
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/eb765d3a-d9bc-4669-928d-01c281f73911%40googlegroups.com.

[deal.II] Re: Bug Report in GridTools::get_active_neighbors()

[Simon Sticko]

Hi,
I think it's just the documentation that can be slightly improved here.

The function get_active_child_cells has a similar signature but it has 
this note in its documentation:

* @note Since in C++ the MeshType template argument can not be deduced from
* a function call, you will have to specify it after the function name, as
* for example in
* @code
*   GridTools::get_active_child_cells > (cell)
* @endcode
*/
template 
std::vector
get_active_child_cells(const typename MeshType::cell_iterator );

The same note should maybe be added to get_active_neighbors for clarity.

Best,
Simon

On Thursday, November 28, 2019 at 10:03:17 AM UTC+1, Andreas Kyritsakis 
wrote:
>
> Hello,
>
> This is to report what I think is a bug in 
> GridTools::get_active_neighbors().
>
> The following code should compile but it does not.
>
> #include 
> #include 
> #include 
> #include 
> using namespace dealii;
>
> template
> class Test{
> public: 
> Test(Triangulation *tria) : triangulation(triangulation), 
> dof_handler(*tria){}
> 
> void run(){
> typename DoFHandler::active_cell_iterator cell;
>   
> for (cell = dof_handler.begin_active(); cell != dof_handler.end(); 
> ++cell){
> std::vector::active_cell_iterator> 
> neighbors;
> GridTools::get_active_neighbors(cell, neighbors); 
> }
> }
> private:
> Triangulation* triangulation;
> DoFHandler dof_handler;
> };
>
> int main(){
> Triangulation<2> triangulation;
> GridGenerator::hyper_cube(triangulation, -1, 1);
> Test<2> test();
> test.run();
> }
>
> I think somehow the compiler cannot resolve the generic type 
> MeshType::active_cell_iterator into the specific one 
> DoFHandler::active_cell_iterator, although from what I understand in 
> the documentation it should. Probably some template specification missing 
> somewhere, but I can't find it.
>
> In my code I fixed it easily by copying the function directly to my own 
> files and specifying appropriately, but I think it would be nice to report 
> the bug.
>
> Cheers,
> Andreas 
>

-- 
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/0a78a60f-89d1-4c6a-9ba3-724d7112c355%40googlegroups.com.

Re: [deal.II] Re: fe_enriched and step-47

[Simon Sticko]

Hi.
Regarding this:

> On Thursday, April 6, 2017 at 1:51:41 PM UTC+2, Denis Davydov wrote:
> 
> The issue with integrating and quadrature would still remain.


I've been implementing an algorithm to generate bulk and surface quadrature 
given a level set function as boundary description. It's mainly an 
implementation of this algorithm: dx.doi.org/10.1137/140966290. 
I'm hoping to finish this within the coming weeks and contribute with this 
to dealii.

Best 
Simon

-- 
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.