from:"Marco Feder"

[deal.II] Re: Issue with intersection of two meshes with curved boundaries

2023-12-05 Thread Marco Feder

After debugging this particular instance, it turns out the issue is related 
to the CGAL kernels used by deal.II. In particular, with CGAL it's possible 
to rely on the so-called "exact computation paradigm" (a brief explanation 
is available here: https://www.cgal.org/exact.html). In the 
dealii::CGALWrappers namespace this kernel is **not** used for the 
quad-quad intersections, which is the one relevant for this example. 
Switching to exact kernels allows to keep the error around 1e-13 for most 
of the refinement cycles with this particular configuration. 

Maybe it could be nice to provide a policy that allows the used to simply 
tweak the CGAL kernels that are employed.


Best,
Marco

Il giorno domenica 26 novembre 2023 alle 09:58:51 UTC+1 najwaa...@gmail.com 
ha scritto:

> Dear team,
>
> I am adding the two domain meshes (cycle 0) in .vtk format for a better 
> understanding of what type of grid I am working with
>
>
> [image: mesh_intersection.png]
>
>
> On Saturday, November 25, 2023 at 10:50:25 AM UTC+3 Najwa Alshehri wrote:
>
>> Dear developers and users,
>>
>>  
>>
>> I have two meshes one is immersed in the other. I wanted to find the 
>> intersection between the two meshes, so I used the following function.
>>
>>  
>>
>> NonMatching::compute_intersection(omega_grid_tools_cache,
>>
>>   omega2_grid_tools_cache,
>>
>>   4,  // degree
>>
>>   1e-20); // to
>>
>>
>> include header 
>> 
>>
>> source code 
>> 
>>
>>
>> This function uses CGAL to find a quadrature formula to integrate exactly 
>> on the intersection of two meshes which neglecte intersections of areas 
>> with tolerance smaller than “tol“ that one chooses and gives a quadrature 
>> formula on the triangulation of the intersection area that integrates 
>> exactly polynomials of a specific degree ( which allows maximum degree of 
>> 4).
>>
>>
>> Say that I want to integrate, on the intersection area, a polynomial of 
>> order 4.
>>
>>  
>>
>> I noticed that If I am considering a circular immersed domain (unlike 
>> square or L-shaped domains), after a few cycles, the quadrature formula is 
>> not accurate enough. To be precise, I find that the sum of the weights of 
>> the quadrature formula defined on the triangulation of the exact 
>> intersection of the two meshes does not sum up (up to a tolerance) to the 
>> measure of the domain. When this occurs, the solution that is solved by 
>> evaluating the integral considering the quadrature formula on the exact 
>> intersection is no longer correct and the error starts to diverge in the 
>> later cycles after this point.
>>
>>  
>>
>> Moreover, the difference gets large suddenly, in one cycle, the 
>> difference was relatively smaller (1e-13), and in the next, it is much 
>> larger (1e-8) as can be seen in the attached plot (Plot shows the 
>> difference between the sum of the weights on the whole domain defined on 
>> the triangulation of the exact intersection of the two meshes and the 
>> measure of the domain under uniform refinement of the mesh). 
>>
>>  
>>
>> Any suggestion on what could be the issue and what should I do to fix it?
>>
>>  
>>
>> Thanks
>>
>> Najwa
>>
>

-- 
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/30085ad8-de31-4513-9370-6f03cd6f37dfn%40googlegroups.com.

[deal.II] Re: Trying to interpolate deal.ii solution to a Cartesian grid

2023-03-08 Thread Marco Feder

Hi,

I've recently used the following function for a similar task: 
https://www.dealii.org/developer/doxygen/deal.II/classFunctions_1_1FEFieldFunction.html#ac25d965867c71d1139262cf383f9f593
 
(see in particular the first snippet). This works under the assumption that 
the cartesian grid is entirely included in the original one. 

Hope that helps,
Marco
Il giorno mercoledì 8 marzo 2023 alle 08:17:17 UTC+1 elyn...@gmail.com ha 
scritto:

> Hello,
>
> I have a solution on an unstructured grid obtained by deal.ii in VTU 
> format, and now I want to interpolate it to a cartesian grid for further 
> processing in other softwares. It seems like Paraview can do the job, but 
> I'm wondering if this can be done in deal.ii directly? Thank you!
>
> Best,
> Elyn
>
>
>

-- 
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/c4afa8b7-3165-440a-9cb1-e5b207beec39n%40googlegroups.com.

[deal.II] Re: Spack error with cgal and arbox

2023-02-22 Thread Marco Feder

Hi Wayne,

CGAL requires C++17, hence you should add a proper compiler flag 
(https://github.com/dealii/dealii/wiki/deal.II-in-Spack#compiler-flags) if 
you want to have it.

Adding *cppflags="-std=c++17"* should fix that. 

For what concerns your last message, here's a related 
discussion: https://groups.google.com/g/dealii/c/nodblNqEjd4

Best,
Marco
Il giorno mercoledì 22 febbraio 2023 alle 13:36:12 UTC+1 yy.wayne ha 
scritto:

> Since spack is calling error, I install directly with the depending 
> libraries installed by spack.
> Deal.II cmake and 'make install' successfully. But shows error when 'make 
> test':
>
> ...
> An error occurred in line <170> of file 
> 
>  
> in function
> void dealii::SparsityTools::{anonymous}::partition_metis(const 
> dealii::SparsityPattern&, const std::vector&, unsigned int, 
> std::vector&)
> The violated condition was:
> ierr == 1
> Additional information: 
> An error with error number -3 occurred while calling a METIS function
> ...
>
> The METIS_DIR was set to metis-5.1.0. I changed METIS_DIR to parmetis, but 
> the error remains.
>
> Besides, I fail to build step1 after this fault installation. cmake passed 
> but make gives following error:
> [image: Snipaste_2023-02-22_20-36-58.png]
>
> 在2023年2月22日星期三 UTC+8 15:16:40 写道：
>
>> The spacl-build-out.txt with middle content deleted. The original file 
>> was too big.
>> 在2023年2月22日星期三 UTC+8 15:12:25 写道：
>>
>>> I set variant "~arborx". There's only ninja warning but installation 
>>> exit with ERROR. However in spack-build-out.txt every libraries seem 
>>> connected successfully. Why ninja return an ERROR?
>>
>>

-- 
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/dda5265f-eabf-4d3d-8223-b7b47800604dn%40googlegroups.com.

[deal.II] Re: Finding cut cells

2022-12-15 Thread Marco Feder

Hi, 
I think Step 85 does this by using a (global) level set function. 

Otherwise, you can use RTrees of bounding boxes 
with 
https://www.dealii.org/developer/doxygen/deal.II/classGridTools_1_1Cache.html#aca2782d6e93b5a0033c046b57904c67f
 
combined with boost::geometry::index to identify intersections, in case 
where you don't want to use a level-set and you just work with the mesh 
itself. The following snippet does this:

*namespace bgi = boost::geometry::index;*
*const auto &tree =*
*space_cache -> get_locally_owned_cell_bounding_boxes_rtree();*
* // Bounding boxes of the embedded grid*
*const auto &embedded_tree = embedded_cache -> 
get_cell_bounding_boxes_rtree();*

*for (const auto &[embedded_box, embedded_cell] : embedded_tree)*
* {*
*for (const auto &[space_box, space_cell] :*
* tree | bgi::adaptors::queried(bgi::intersects(embedded_box)))*
* {*
* // do what you need here*
* }*
* }*

Best,
Marco

Il giorno giovedì 15 dicembre 2022 alle 05:55:51 UTC+1 Praveen C ha scritto:

> Dear all
>
> I am looking for available methods to find cells cut by some given curve.
>
> I found such a case in step-60 and made this example
>
> https://bitbucket.org/cpraveen/deal_ii/src/master/embedded_curve/
>
> As shown in the figures, some cells are not identified though they are cut 
> by the curve. This is a limitation which is already explained in step-60.
>
> Are there other methods available in deal.II for this purpose ? Or do you 
> know of any external library than can be used for this, and which can be 
> called from a deal.II code ?
>
> I am assuming that the curve is given as piecewise polynomials. I only 
> need 2-D case.
>
> Thanks
> praveen

-- 
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/a909081d-b3d5-48bc-86a8-dacf2ed049d9n%40googlegroups.com.

Re: [deal.II] Tutorial 40 , error: use of deleted function

2022-11-07 Thread Marco Feder

Ah, I hadn’t noticed you were compiling a tutorial program. My apologies.

Il giorno lunedì 7 novembre 2022 alle 17:53:10 UTC+1 d.arnd...@gmail.com ha 
scritto:

> You would get this compile-time error if you haven't configured deal.II 
> with `DEAL_II_WITH_MPI=ON`.
>
> Best,
> Daniel
>
> On Mon, Nov 7, 2022 at 9:51 AM ztdep...@gmail.com  
> wrote:
>
>> when i compile the tutorial 40, i met the following error. It semms that 
>> the function has been deleted. 
>>
>>  error: use of deleted function 
>> ‘dealii::parallel::distributed::Triangulation> spacedim>::Triangulation(const MPI_Comm&, typename 
>> dealii::Triangulation::MeshSmoothing, 
>> dealii::parallel::distributed::Triangulation::Settings) 
>> [with int dim = 2; int spacedim = 2; MPI_Comm = int; typename 
>> dealii::Triangulation::MeshSmoothing = 
>> dealii::Triangulation<2, 2>::MeshSmoothing]’
>>
>> -- 
>> 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+un...@googlegroups.com.
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/6682d965-46f7-4ca7-b3e6-a5a365e138ffn%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/1fefb1cf-d954-4c43-acd2-752cfc6e633dn%40googlegroups.com.

[deal.II] Re: Tutorial 40 , error: use of deleted function

2022-11-07 Thread Marco Feder

Hi,

the copy constructor for Triangulation is indeed deleted. If you do need a 
copy, you may want to use copy_triangulation() 
,
 
that is also implemented for p::d::T.


Best,
Marco
Il giorno lunedì 7 novembre 2022 alle 15:51:17 UTC+1 ztdep...@gmail.com ha 
scritto:

> when i compile the tutorial 40, i met the following error. It semms that 
> the function has been deleted. 
>
>  error: use of deleted function 
> ‘dealii::parallel::distributed::Triangulation spacedim>::Triangulation(const MPI_Comm&, typename 
> dealii::Triangulation::MeshSmoothing, 
> dealii::parallel::distributed::Triangulation::Settings) 
> [with int dim = 2; int spacedim = 2; MPI_Comm = int; typename 
> dealii::Triangulation::MeshSmoothing = 
> dealii::Triangulation<2, 2>::MeshSmoothing]’
>

-- 
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/3c20e3a3-1c5e-410a-a7ed-329836c2f19cn%40googlegroups.com.

Re: [deal.II] How to implemente the Neumann bc for laplacian

2022-08-26 Thread Marco Feder

With your u(x,y) you get as (manufactured) rhs:

f(x,y) \approx 1.8 e^(-1.9 x)

I don’t think that with these data your condition

\int_\Omega f + \int_{\partial \Omega} \partial_n u = 0

is satisfied.


Best,
Marco

> On 26 Aug 2022, at 14:42, ztdep...@gmail.com  wrote:
> 
> I solve Laplacian problem with "-0.5*std::exp(-2.0*0.963740544195800*p(0))" 
> as exact solution, I apply Neumaan bc on all 4 boundaries. 
> the solver feedback the error. no solution is obtained. 
> 
> On Friday, August 26, 2022 at 8:01:11 PM UTC+8 marco@gmail.com wrote:
> What are the data (rhs, Neumann condition,...) of your problem ?
> What is the value of the residual ?
> 
> Best,
> Marco
> 
> 
>> On 26 Aug 2022, at 13:48, ztdep...@gmail.com <http://gmail.com/> 
>> > 
>> wrote:
>> 
> 
>> it seems I cann't get a converged solution . which Solver  should i use 
>> to treat this situation?The error message is as follows .
>> 
>> 
>>This error message can indicate that you have simply not allowed a
>> sufficiently large number of iterations for your iterative solver to
>> converge. This often happens when you increase the size of your
>> problem. In such cases, the last residual will likely still be very
>> small, and you can make the error go away by increasing the allowed
>> number of iterations when setting up the SolverControl object that
>> determines the maximal number of iterations you allow.
>> 
>> The other situation where this error may occur is when your matrix is
>> not invertible (e.g., your matrix has a null-space), or if you try to
>> apply the wrong solver to a matrix (e.g., using CG for a matrix that
>> is not symmetric or not positive definite). In these cases, the
>> residual in the last iteration is likely going to be large.
>> 
>> 
>> On Wednesday, August 24, 2022 at 12:56:59 AM UTC+8 marco@gmail.com 
>> <http://gmail.com/> wrote:
>> Actually, it turns out there's a tutorial program doing precisely what you 
>> are looking for: step-11 
>> (https://www.dealii.org/developer/doxygen/deal.II/step_11.html 
>> <https://www.dealii.org/developer/doxygen/deal.II/step_11.html>)
>> 
>> Best,
>> Marco
>> 
>> Il giorno martedì 23 agosto 2022 alle 18:54:24 UTC+2 Marco Feder ha scritto:
>> 
>> Hi,
>> 
>> step-7 (https://www.dealii.org/current/doxygen/deal.II/step_7.html 
>> <https://www.dealii.org/current/doxygen/deal.II/step_7.html>) shows this 
>> procedure. Be careful that if you apply Neumann BCs on all your boundary, 
>> then the solution is determined up to a constant and you need to add an 
>> additional constraint to your solution. Also, you need to satisfy a 
>> compatibility condition between the r.h.s. and the neumann data.
>> 
>> Best,
>> Marco
>> Il giorno martedì 23 agosto 2022 alle 15:54:19 UTC+2 ztdep...@gmail.com <> 
>> ha scritto:
>> In step 3 , I want to implemente the Neumann bc for all the boundary 
>> condition. which tutorial show this procedure.
>> 
> 
>> -- 
>> 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 a topic in the 
>> Google Groups "deal.II User Group" group.
>> To unsubscribe from this topic, visit 
>> https://groups.google.com/d/topic/dealii/0KK8UDeFgRg/unsubscribe 
>> <https://groups.google.com/d/topic/dealii/0KK8UDeFgRg/unsubscribe>.
>> To unsubscribe from this group and all its topics, send an email to 
>> dealii+un...@googlegroups.com 
>> .
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/8262a85f-01b0-4b9e-b1f8-ee5a0b3361b3n%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/8262a85f-01b0-4b9e-b1f8-ee5a0b3361b3n%40googlegroups.com?utm_medium=email&utm_source=footer>.
> 
> 
> 
> -- 
> 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 a topic in the Google 
> Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit 
> https://groups.google.com/d/topic/dealii/0K

Re: [deal.II] How to implemente the Neumann bc for laplacian

2022-08-26 Thread Marco Feder

What are the data (rhs, Neumann condition,...) of your problem ?
What is the value of the residual ?

Best,
Marco

> On 26 Aug 2022, at 13:48, ztdep...@gmail.com  wrote:
> 
> it seems I cann't get a converged solution . which Solver  should i use 
> to treat this situation?The error message is as follows .
> 
> 
>This error message can indicate that you have simply not allowed a
> sufficiently large number of iterations for your iterative solver to
> converge. This often happens when you increase the size of your
> problem. In such cases, the last residual will likely still be very
> small, and you can make the error go away by increasing the allowed
> number of iterations when setting up the SolverControl object that
> determines the maximal number of iterations you allow.
> 
> The other situation where this error may occur is when your matrix is
> not invertible (e.g., your matrix has a null-space), or if you try to
> apply the wrong solver to a matrix (e.g., using CG for a matrix that
> is not symmetric or not positive definite). In these cases, the
> residual in the last iteration is likely going to be large.
> 
> 
> On Wednesday, August 24, 2022 at 12:56:59 AM UTC+8 marco@gmail.com wrote:
> Actually, it turns out there's a tutorial program doing precisely what you 
> are looking for: step-11 
> (https://www.dealii.org/developer/doxygen/deal.II/step_11.html 
> <https://www.dealii.org/developer/doxygen/deal.II/step_11.html>)
> 
> Best,
> Marco
> 
> Il giorno martedì 23 agosto 2022 alle 18:54:24 UTC+2 Marco Feder ha scritto:
> 
> Hi,
> 
> step-7 (https://www.dealii.org/current/doxygen/deal.II/step_7.html 
> <https://www.dealii.org/current/doxygen/deal.II/step_7.html>) shows this 
> procedure. Be careful that if you apply Neumann BCs on all your boundary, 
> then the solution is determined up to a constant and you need to add an 
> additional constraint to your solution. Also, you need to satisfy a 
> compatibility condition between the r.h.s. and the neumann data.
> 
> Best,
> Marco
> Il giorno martedì 23 agosto 2022 alle 15:54:19 UTC+2 ztdep...@gmail.com <> ha 
> scritto:
> In step 3 , I want to implemente the Neumann bc for all the boundary 
> condition. which tutorial show this procedure.
> 
> -- 
> 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 a topic in the Google 
> Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit 
> https://groups.google.com/d/topic/dealii/0KK8UDeFgRg/unsubscribe 
> <https://groups.google.com/d/topic/dealii/0KK8UDeFgRg/unsubscribe>.
> To unsubscribe from this group and all its topics, 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/8262a85f-01b0-4b9e-b1f8-ee5a0b3361b3n%40googlegroups.com
>  
> <https://groups.google.com/d/msgid/dealii/8262a85f-01b0-4b9e-b1f8-ee5a0b3361b3n%40googlegroups.com?utm_medium=email&utm_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/E7EF3FF0-DC88-4B9D-B49E-7CB8FDF8C32F%40gmail.com.

[deal.II] Re: How to implemente the Neumann bc for laplacian

2022-08-23 Thread Marco Feder

Actually, it turns out there's a tutorial program doing precisely what you 
are looking for: step-11 
(https://www.dealii.org/developer/doxygen/deal.II/step_11.html)

Best,
Marco

Il giorno martedì 23 agosto 2022 alle 18:54:24 UTC+2 Marco Feder ha scritto:

>
> Hi,
>
> step-7 (https://www.dealii.org/current/doxygen/deal.II/step_7.html) shows 
> this procedure. Be careful that if you apply Neumann BCs on all your 
> boundary, then the solution is determined up to a constant and you need to 
> add an additional constraint to your solution. Also, you need to satisfy a 
> compatibility condition between the r.h.s. and the neumann data.
>
> Best,
> Marco
> Il giorno martedì 23 agosto 2022 alle 15:54:19 UTC+2 ztdep...@gmail.com 
> ha scritto:
>
>> In step 3 , I want to implemente the Neumann bc for all the boundary 
>> condition. which tutorial show this procedure.
>
>

-- 
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/e87f2511-a0e9-4aa8-8f78-8adc143ac563n%40googlegroups.com.

[deal.II] Re: How to implemente the Neumann bc for laplacian

2022-08-23 Thread Marco Feder


Hi,

step-7 (https://www.dealii.org/current/doxygen/deal.II/step_7.html) shows 
this procedure. Be careful that if you apply Neumann BCs on all your 
boundary, then the solution is determined up to a constant and you need to 
add an additional constraint to your solution. Also, you need to satisfy a 
compatibility condition between the r.h.s. and the neumann data.

Best,
Marco
Il giorno martedì 23 agosto 2022 alle 15:54:19 UTC+2 ztdep...@gmail.com ha 
scritto:

> In step 3 , I want to implemente the Neumann bc for all the boundary 
> condition. which tutorial show this procedure.

-- 
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/5b956f60-a0a3-4991-966d-2544373dfb8dn%40googlegroups.com.

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

2022-08-18 Thread Marco Feder

Thanks Simon,

As you might suspect, I was using your step-85 to solve an interface problem 
also in the “outside” region with cutFEM, and  “hardcoding" the jump terms like 
in DG as you suggested finally solved the issue. Sorry for the late reply. 

Best, 
Marco

> On 9 Aug 2022, at 17:59, Simon Sticko  wrote:
> 
> 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/ 
>> <http://www.dealii.org/> <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> 
>> <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> 
>> <mailto:dealii+unsubscr...@googlegroups.com 
>> <mailto: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
>>  
>> <https://groups.google.com/d/msgid/dealii/2d71b802-132c-4f8a-a20d-03542939036an%40googlegroups.com><https://groups.google.com/d/msgid/dealii/2d71b802-132c-4f8a-a20d-03542939036an%40googlegroups.com?utm_medium=email&utm_source=footer
>>  
>> <https://groups.google.com/d/msgid/dealii/2d71b802-132c-4f8a-a20d-03542939036an%40googlegroups.com?utm_medium=email&utm_source=footer>>.
> 
> -- 
> 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 a topic in the 
> Google Groups "deal.II User Group" group.
> To unsubscribe from this topic, visit 
> https://groups.google.com/d/topic/dealii/4nmXUhsSNfc/unsubscribe 
> <https://groups.google.com/d/topic/dealii/4nmXUhsSNfc/unsubscribe>.
> To unsubscribe from this group and all its topics, 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/23e81ebd-cf4e-60ad-ea66-ad825581b533%40gmail.com
>  
> <https://groups.google.com/d/msgid/dealii/23e81ebd-cf4e-60ad-ea66-ad825581b533%40gmail.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/B0BCA4D3-48AD-4D75-9AF1-1E1134B9A5B9%40gmail.com.

[deal.II] FEInterfaceValues and hp::FECollection

2022-08-08 Thread Marco Feder

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.

[deal.II] Re: dealii-9.4.0 release candidate 1

2022-07-01 Thread Marco Feder

This may be of interest to 
you: https://github.com/dealii/dealii/issues/13574#issuecomment-1166472620

Best,
Marco

Il giorno giovedì 30 giugno 2022 alle 13:55:34 UTC+2 Alberto Salvadori ha 
scritto:

> Matthias,
>
> I wonder if this release is  fully compatible with Apple M1 (and M2 as 
> well) or not. 
> Thank you,
>
> Alberto
>
> Il giorno giovedì 16 giugno 2022 alle 23:04:32 UTC+2 Matthias Maier ha 
> scritto:
>
>> Dear all, 
>>
>> We have tagged a first release candidate for the upcoming deal.II 9.4.0 
>> release. 
>> You can find the sources and a generated offline documentation here: 
>>
>> https://github.com/dealii/dealii/releases/tag/v9.4.0-rc1 
>>
>> Alternatively, you can switch to the dealii-9.4 branch in your local git 
>> repository 
>>
>> git remote update 
>> git checkout v9.4.0-rc1 
>>
>> It would be great if you could test it on your machine with your typical 
>> configuration! 
>> If no further regressions show up, we will release Friday June 24. The 
>> release 
>> progress is tracked in the following GitHub issue: 
>>
>> https://github.com/dealii/dealii/issues/13574 
>>
>> Thanks! 
>>
>> Matthias 
>> on behalf of the deal.II developer team 
>>
>

-- 
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/a1859612-f976-4696-97df-21c6d39a1fb9n%40googlegroups.com.

Re: [deal.II] convert linear_operator to SparseMatrix or get the inverse of linear_operator

2022-04-09 Thread Marco Feder

Hi Chen,

You can work directly with LinearOperator(s). What you need is the *inverse 
operator* of op_M , and you can get it by giving a Solver and a 
corresponding preconditioner to inverse_operator() 
[https://www.dealii.org/current/doxygen/deal.II/group__LAOperators.html#ga87e38fbde431397c069a88692bd24ae7],
 
as explained in step-20. 
(https://www.dealii.org/current/doxygen/deal.II/step_20.html#TheLinearOperatorframeworkindealII)


Best,
Marco
Il giorno sabato 9 aprile 2022 alle 05:33:22 UTC+2 hkch...@gmail.com ha 
scritto:

> Hello team,
> First, I got the linear_operator op_M(the matrix M). I want to get the 
> inverse of M named as M_inv. Because of the difficult of inverting the M 
> and that M is not symmetry, I used such codes to get M_inv:
>
>> SparseMatrix M = op_M;
>
> solver SolverGMRES> gmres(solver_control); 
>
> const auto M_inv = linear_operator(M, gmres, somepreconditioner);
>
>  There are some questions on it：
>
>1. I don't know how to convert linear_operator to 
>SparseMatrix, the code "SparseMatrix M = op_M;" get 
>error.
>2. Whether could I get the M_inv directly by op_M using iterator 
>method  and don't need to convert linear_operator to SparseMatrix?
>3. Where I can find the support preconditioner for the iterator? it 
>seems that the SolverGMRES don't supports the preconditoner of 
>SparseILU and SparseDirectUMFPACK. 
>
> best regards
> chen
>

-- 
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/ffc215e1-742b-4074-8b80-7e3ebfb365b5n%40googlegroups.com.

[deal.II] Re: Integral over different lines

2022-03-12 Thread Marco Feder

Hi Giuseppe,

if your goal is to compute the integral of an arbitrary function f(x,y) 
over a random segment inside your rectangle (you can assume the segment is 
not parallel to x-axis or y-axis), you could just consider its 
parametrisation \phi(t), t \in [0,1] so that you have a 1D integral that 
you can compute. Weights and points in the 1D reference domain can then be 
asked to dealii.

If your line is parallel to one of the axes, you may want to use 
compute_affine_transformation() (
https://www.dealii.org/current/doxygen/deal.II/classQSimplex.html#ab670894e45080a2c39d0698c0f582c9d)
 
to get points and weights in the "real" line, since one coordinate is fixed.

More advanced users probably have nicer solutions :-)


Best,
Marco

Il giorno venerdì 11 marzo 2022 alle 16:26:05 UTC+1 gius...@gmail.com ha 
scritto:

> Good morning everyone,
> I would have a question about the computation of integrals. Suppose I have 
> a rectangular domain, is it there any way to compute the integral of a 
> certain function over different heights (not only on the bottom or top 
> boundary)?
>
> Thanks in advance,
> Giuseppe
>

-- 
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/9713bc9d-65c9-4008-8b96-169f92a87ea9n%40googlegroups.com.

[deal.II] Re: Integral over different lines

2022-03-12 Thread Marco Feder

Hi Giuseppe,

if your goal is to compute the integral an arbitrary function f(x,y) over a 
random segment inside your rectangle (you can assume the segment is not 
parallel to x-axis or y-axis), you could just consider its parametrisation 
\phi(t), t \in [0,1] so that you have a 1D integral that you can compute by 
quadrature. Weights and points in the 1D reference domain can then be asked 
to dealii.

If your line is parallel to one of the axes, you may want to use 
compute_affine_transformation() 
(https://www.dealii.org/current/doxygen/deal.II/classQSimplex.html#ab670894e45080a2c39d0698c0f582c9d)
 
to get points and weights in the "real" line, since one coordinate is fixed.

More advanced users probably have nicer solutions :-)


Best,
Marco

Il giorno venerdì 11 marzo 2022 alle 16:26:05 UTC+1 gius...@gmail.com ha 
scritto:

> Good morning everyone,
> I would have a question about the computation of integrals. Suppose I have 
> a rectangular domain, is it there any way to compute the integral of a 
> certain function over different heights (not only on the bottom or top 
> boundary)?
>
> Thanks in advance,
> Giuseppe
>

-- 
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/0ab40b18-4b0e-4ea3-8b7f-6fa589431f32n%40googlegroups.com.

Re: [deal.II] Extractors for local_rhs in step-22

2021-11-19 Thread Marco Feder

Thanks Wolfgang,

indeed, as we have incompressibility, the last component of the rhs will 
always be 0. I wrote a possible patch, I'll open the pr tomorrow so that I 
can take another look at it!

Best,
Marco

Il giorno venerdì 19 novembre 2021 alle 23:15:50 UTC+1 Wolfgang Bangerth ha 
scritto:

>
> Marco,
>
> > I was going through step-22 and I saw that since you're using primitive 
> > elements, then you compute the local_rhs contribution using 
> > /fe.system_to_component_index(i).first/
> > 
> > However, as written in the program, we could as well multiply the dim+1 
> tensor 
> > having the values of the i-th shape function with the whole rhs vector 
> (f,0).
> > 
> > So, I have two questions:
> > 1) Are you using /fe.system_to_component_index(i).first /just to make 
> the 
> > program run faster? Why not using extractors?
>
> I think it is an oversight. I would accept a patch to correct that :-)
>
>
> > 2) Taking the scalar product between a rank-1 Tensor and a Vector with 
> the 
> > same size is not "conceptually" allowed. Did you have in mind something 
> like 
> > the following?
> > 
> > //inside assemble_system()
> > std::vector> phi_u(dofs_per_cell);
> > for (unsigned int q = 0; q < n_q_points; ++q)
> > {
> > for (unsigned int k = 0; k < dofs_per_cell; ++k)
> > {
> > //store symgrad_phi_u, div_phi_u, etc. using extractors
> > phi_u[k] = fe_values[velocities].value(k, q);
> > }
> > // compute local_matrix and local_preconditioner
> > 
> > double sum = 0.0; //store scalar product between first dim components of 
> i-th 
> > shape function (corresponding to velocity) and first dim components of 
> the rhs
> > for (unsigned int s = 0; s < dim; ++s)
> > {
> > sum += phi_u[i][s] * rhs_values[q][s];
> > }
> > local_rhs(i) +=
> > (sum + rhs_values[q][dim] * phi_p[i]) * fe_values.JxW(q);
> > }
> > }
>
> That would work. I think a better solution would be to say outright that 
> only 
> the first equation can have force terms, and then make the RightHandSide 
> class 
> derived from TensorFunction<1,dim> and let it return its values as 
> Tensor<1,dim> objects for which you can take the dot product
> fe_values[u].shape_value(q) * rhs_values[q]
>
> Want to give this a try and write a patch? We'd be happy to walk you 
> through 
> what is necessary for such a patch!
>
> 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.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/d83ed03f-5efd-4d94-a97b-b250daa51954n%40googlegroups.com.

[deal.II] Extractors for local_rhs in step-22

2021-11-18 Thread Marco Feder

Dear all,

I was going through step-22 and I saw that since you're using primitive 
elements, then you compute the local_rhs contribution using 
*fe.system_to_component_index(i).first*

However, as written in the program, we could as well multiply the dim+1 
tensor having the values of the i-th shape function with the whole rhs 
vector (f,0). 

So, I have two questions: 
1) Are you using *fe.system_to_component_index(i).first *just to make the 
program run faster? Why not using extractors?

2) Taking the scalar product between a rank-1 Tensor and a Vector with the 
same size is not "conceptually" allowed. Did you have in mind something 
like the following?

//inside assemble_system()
std::vector> phi_u(dofs_per_cell);
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int k = 0; k < dofs_per_cell; ++k)
{
//store symgrad_phi_u, div_phi_u, etc. using extractors
phi_u[k] = fe_values[velocities].value(k, q);
}
// compute local_matrix and local_preconditioner

double sum = 0.0; //store scalar product between first dim components of 
i-th shape function (corresponding to velocity) and first dim components of 
the rhs
for (unsigned int s = 0; s < dim; ++s)
{
sum += phi_u[i][s] * rhs_values[q][s];
}
local_rhs(i) +=
(sum + rhs_values[q][dim] * phi_p[i]) * fe_values.JxW(q);
}
}


or is there a better way to achieve 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/10319c6f-57c3-44ab-af4e-ca1a6cab10d6n%40googlegroups.com.

Re: [deal.II] Inaccurate convergence rate using ParsedConvergenceTable class

2021-08-09 Thread Marco Feder

You're right, I totally missed the default parameter h in the constructor.

Thanks,
Marco
Il giorno domenica 8 agosto 2021 alle 22:26:01 UTC+2 Timo Heister ha 
scritto:

> FunctionParser uses a finite difference approximation for the gradient. I 
> think this is explained in the class documentation.
> This explains the results you see...
>
> On Sun, Aug 8, 2021, 15:15 Marco Feder  wrote:
>
>> Ciao Luca!
>>
>> > Can you check if you get the same results with this order?
>> Yes, the rates are the same. 
>>
>> > Are you calling `error_from_exact` with mapping as a first argument?
>> Yes
>>
>> Actually, I've fixed the issue after reading this:
>> > If the Solution class implements the Gradient, then 
>> ParsedConvergenceTable should use that. 
>>
>> Indeed, I realised I wrote:
>> error_table.error_from_exact(mapping, dof_handler, solution, 
>> *exact_solution*);
>>
>> where exact_solution is a FunctionParser<*dim*> which will be 
>> initialised with the expression of the analytical solution, so it doesn't 
>> give any info about the Gradient. Instead, if I replace exact_solution in 
>> the last argument with the "classical" Solution() (so that there's 
>> also an implementation of Gradient) I have my expected rates. Apparently it 
>> wasn't using the Gradient function. 
>>
>> However, looking at the source code and in the documentation of 
>> integrate_difference, I don't see how the H^1 norm was computed without 
>> having an implementation of the Gradient. I would expect to see an 
>> exception asking for its implementation. Was it using some FD-like 
>> approximation or am I missing something?
>>
>>
>> Thanks,
>> Marco
>>
>>
>> Il giorno domenica 8 agosto 2021 alle 18:52:03 UTC+2 luca@gmail.com 
>> ha scritto:
>>
>>> The other option is that you are not calling the `error_from_exact` 
>>> function that takes a mapping argument, but the one without it, and your 
>>> boundary cells are introducing some error due to a worse mapping…. Are you 
>>> calling `error_from_exact` with mapping as a first argument?
>>>
>>> Luca
>>>
>>> Il giorno 8 ago 2021, alle ore 6:49 PM, Luca Heltai  
>>> ha scritto:
>>>
>>> 
>>>
>>> Ciao Marco,
>>>
>>> I was expanding step-12 
>>> <https://www.dealii.org/current/doxygen/deal.II/step_12.html> using a 
>>> manufactured solution to check the order p in H^1 (and p+1 in L^2) norm on 
>>> a uniformly refined mesh for the DG upwind method. I've used the 
>>> ParsedConvergenceTable with a parameter file as described in the docs. As 
>>> Rate_key I'm using the DoFs, while as Rate_mode I have reduction_rate_log2. 
>>>
>>> With p=1 and p=2 everything is fine. But if I set the finite element 
>>> degree to 3, then the H^1 convergence rate decreases, as you can see in the 
>>> attached image.
>>> 
>>>
>>>
>>> Is it possible that you are using different quadrature rules in the two 
>>> cases? Your image shows a deterioration of the error on the order of 1e-8 
>>> for H1, and 1e-12 for L2, which is very close to machine precision.
>>>
>>> Internally the parsed convergence table does the exact same thing you 
>>> wrote explicitly. (If you check the source code, you’ll see that it calls 
>>> integrate_difference and compute_global_error for each error type you 
>>> specify in the parameter file).
>>>
>>> This, however, doesn't happen if I use a classical ConvergenceTable. 
>>> Namely, I first compute the local error in each cell, and then the global 
>>> error in the classical way:
>>>
>>> VectorTools::integrate_difference(mapping,dof_handler,solution, Solution<
>>> *dim*>(),H1_error_per_cell, QGauss<*dim*
>>> >(fe->tensor_degree()+1),VectorTools::H1_norm);
>>>
>>> *const* *double* H1_error = 
>>> VectorTools::compute_global_error(triangulation, H1_error_per_cell,  
>>> VectorTools::H1_norm); //assuming I provided also the gradient method for 
>>> the Solution class
>>>
>>>
>>> Does anyone have any idea why this is happening? My guess was that while 
>>> computing the H^1 *semi-*norm the ParsedConvergenceTable class does 
>>> some approximation to compute the gradient from the exact solution 
>>> expression and hence that could be the source of the issue. Conversely, in 
>>> the &

Re: [deal.II] Inaccurate convergence rate using ParsedConvergenceTable class

2021-08-08 Thread Marco Feder

Ciao Luca!

> Can you check if you get the same results with this order?
Yes, the rates are the same. 

> Are you calling `error_from_exact` with mapping as a first argument?
Yes

Actually, I've fixed the issue after reading this:
> If the Solution class implements the Gradient, then 
ParsedConvergenceTable should use that. 

Indeed, I realised I wrote:
error_table.error_from_exact(mapping, dof_handler, solution, 
*exact_solution*);

where exact_solution is a FunctionParser<*dim*> which will be 
initialised with the expression of the analytical solution, so it doesn't 
give any info about the Gradient. Instead, if I replace exact_solution in 
the last argument with the "classical" Solution() (so that there's 
also an implementation of Gradient) I have my expected rates. Apparently it 
wasn't using the Gradient function. 

However, looking at the source code and in the documentation of 
integrate_difference, I don't see how the H^1 norm was computed without 
having an implementation of the Gradient. I would expect to see an 
exception asking for its implementation. Was it using some FD-like 
approximation or am I missing something?


Thanks,
Marco


Il giorno domenica 8 agosto 2021 alle 18:52:03 UTC+2 luca@gmail.com ha 
scritto:

> The other option is that you are not calling the `error_from_exact` 
> function that takes a mapping argument, but the one without it, and your 
> boundary cells are introducing some error due to a worse mapping…. Are you 
> calling `error_from_exact` with mapping as a first argument?
>
> Luca
>
> Il giorno 8 ago 2021, alle ore 6:49 PM, Luca Heltai  
> ha scritto:
>
> 
>
> Ciao Marco,
>
> I was expanding step-12 
>  using a 
> manufactured solution to check the order p in H^1 (and p+1 in L^2) norm on 
> a uniformly refined mesh for the DG upwind method. I've used the 
> ParsedConvergenceTable with a parameter file as described in the docs. As 
> Rate_key I'm using the DoFs, while as Rate_mode I have reduction_rate_log2. 
>
> With p=1 and p=2 everything is fine. But if I set the finite element 
> degree to 3, then the H^1 convergence rate decreases, as you can see in the 
> attached image.
> 
>
>
> Is it possible that you are using different quadrature rules in the two 
> cases? Your image shows a deterioration of the error on the order of 1e-8 
> for H1, and 1e-12 for L2, which is very close to machine precision.
>
> Internally the parsed convergence table does the exact same thing you 
> wrote explicitly. (If you check the source code, you’ll see that it calls 
> integrate_difference and compute_global_error for each error type you 
> specify in the parameter file).
>
> This, however, doesn't happen if I use a classical ConvergenceTable. 
> Namely, I first compute the local error in each cell, and then the global 
> error in the classical way:
>
> VectorTools::integrate_difference(mapping,dof_handler,solution, Solution<
> *dim*>(),H1_error_per_cell, QGauss<*dim*
> >(fe->tensor_degree()+1),VectorTools::H1_norm);
>
> *const* *double* H1_error = 
> VectorTools::compute_global_error(triangulation, H1_error_per_cell,  
> VectorTools::H1_norm); //assuming I provided also the gradient method for 
> the Solution class
>
>
> Does anyone have any idea why this is happening? My guess was that while 
> computing the H^1 *semi-*norm the ParsedConvergenceTable class does some 
> approximation to compute the gradient from the exact solution 
> expression and hence that could be the source of the issue. Conversely, in 
> the "ConvergenceTable" way I do define explicitly the gradient of the exact 
> solution in the Solution class.
>
>
> If the Solution class implements the Gradient, then 
> ParsedConvergenceTable should use that. 
>
> You are calling “error_from_exact” of that class, right?
>
> This is where it is called:
>
>
> https://www.dealii.org/current/doxygen/deal.II/parsed__convergence__table_8h_source.html
>
> Line 621. As you see, the quadrature rule used is a Gauss formula of order 
> (dh.get_fe().degree+1)*2.
>
> Can you check if you get the same results with this order?
>
> L.
>
>

-- 
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/dc72f870-1549-471c-b93b-854bfe3b67d9n%40googlegroups.com.

[deal.II] Re: Issue with intersection of two meshes with curved boundaries

[deal.II] Re: Trying to interpolate deal.ii solution to a Cartesian grid

[deal.II] Re: Spack error with cgal and arbox

[deal.II] Re: Finding cut cells

Re: [deal.II] Tutorial 40 , error: use of deleted function

[deal.II] Re: Tutorial 40 , error: use of deleted function

Re: [deal.II] How to implemente the Neumann bc for laplacian

Re: [deal.II] How to implemente the Neumann bc for laplacian

[deal.II] Re: How to implemente the Neumann bc for laplacian

[deal.II] Re: How to implemente the Neumann bc for laplacian

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

[deal.II] FEInterfaceValues and hp::FECollection

[deal.II] Re: dealii-9.4.0 release candidate 1

Re: [deal.II] convert linear_operator to SparseMatrix or get the inverse of linear_operator

[deal.II] Re: Integral over different lines

[deal.II] Re: Integral over different lines

Re: [deal.II] Extractors for local_rhs in step-22

[deal.II] Extractors for local_rhs in step-22

Re: [deal.II] Inaccurate convergence rate using ParsedConvergenceTable class

Re: [deal.II] Inaccurate convergence rate using ParsedConvergenceTable class

20 matches

Site Navigation

Mail list logo

Footer information