Re: [deal.II] Re: Tips on writing "versatile" assembly function

[Corbin Foucart]

Hello everyone, 

I"m encountering a similar question for large 3D incompressible fluid 
solver with many assembly options (which I implement via inheritance). I'd 
like to investigate how much class/templated inheritance plays a role in 
performance---Bruno, how do you go about profiling something like this as 
you suggest (without having to implement an entirely new class free of 
inheritance for comparison)?

Thank you,
Corbin

On Tuesday, January 5, 2021 at 6:33:11 PM UTC-5 Doug Shi-Dong wrote:

> That's interesting. Seems like it was more about inlining than branch 
> prediction. Surprising how much difference it made.
>
> On Tuesday, January 5, 2021 at 6:18:38 PM UTC-5 Timo Heister wrote:
>
>> What I forgot to say: 
>> We used to have something like 
>>
>> if (use_anisotropic_viscosity==true) 
>> cell(i,j) += viscosity_tensor *  
>> else 
>> cell(i,j) += viscosity_constant *  
>>
>> and improved the performance by making two separate assemblers (note 
>> that there is no function call/vtable lookup here, just an "if"). I 
>> don't remember how big the difference was, so I went back and found 
>> the PR: 
>> https://github.com/geodynamics/aspect/pull/2139 
>> Rene claims 25% fewer instructions. 
>>
>>
>>
>> On Tue, Jan 5, 2021 at 4:48 PM blais...@gmail.com  
>> wrote: 
>> > 
>> > Hi Timo, 
>> > 
>> > I understand. It makes a lot of sense. 
>> > Thanks! 
>> > Bruno 
>> > 
>> > On Tuesday, January 5, 2021 at 4:34:19 p.m. UTC-5 Timo Heister wrote: 
>> >> 
>> >> Hi Bruno, 
>> >> 
>> >> We mitigate the performance problem by making the decision per cell in 
>> ASPECT: 
>> >> We have a set of "assemblers" that are executed one after each other 
>> >> per cell. This means the vtable access cost is small compared to the 
>> >> actual work. 
>> >> See 
>> >> 
>> https://github.com/geodynamics/aspect/blob/b9add5f53172aac18bdbb19d14ca266e88c491dd/include/aspect/simulator/assemblers/interface.h#L493-L515
>>  
>> >> 
>> >> On Tue, Jan 5, 2021 at 4:28 PM blais...@gmail.com  
>> wrote: 
>> >> > 
>> >> > Bruno, 
>> >> > 
>> >> > Thanks, you are right. As always, measure first and then optimize 
>> after. No point in optimising stuff that costs nothing... 
>> >> > 
>> >> > 
>> >> > On Tuesday, January 5, 2021 at 3:15:06 p.m. UTC-5 
>> bruno.t...@gmail.com wrote: 
>> >> >> 
>> >> >> Bruno, 
>> >> >> 
>> >> >> If you are worry about the cost of looking up though the vtable, I 
>> think that you are stuck using template. So either use 2 or 3 and CRTP. But 
>> first of all, I think that you should profile your code and verify that 
>> this is a problem. There is no point in spending time refactoring your 
>> code, if your are going to gain less than 1%... 
>> >> >> 
>> >> >> Best, 
>> >> >> 
>> >> >> Bruno 
>> >> >> 
>> >> >> On Monday, January 4, 2021 at 3:31:45 PM UTC-5 blais...@gmail.com 
>> wrote: 
>> >> >>> 
>> >> >>> Dear all, 
>> >> >>> I wish you all an happy new year! 
>> >> >>> 
>> >> >>> One problem we always end up facing with FEM problems is that, as 
>> program grow, more and more features are added to the equations. This leads 
>> to multiple variation of the same equations (for example, Navier-Stokes 
>> with Newtonian and non-Newtonian viscosity, etc.). In turn, this leads to 
>> assembly functions for the non-linear systems of equations which branch out 
>> to multiple possibilities. 
>> >> >>> 
>> >> >>> I was wondering what is the best approach to deal with multiple 
>> "switches" in an assembly function. The naïve approach would be to use if 
>> conditions, but I have a feeling that if they appear down in the assembly 
>> (for example the loop on dofs), this would lead to significantly higher 
>> assembly cost because the loops would spend time just checking if. 
>> >> >>> 
>> >> >>> From experience, I have seen the following approaches: 
>> >> >>> 1- If or switches in the assembly routine. Simplest/most versatile 
>> way, but adds significant overhead 
>> >> >>> 2- Template the assembly function with the parameters. I think 
>> this would lead to zero additional cost, but as the number of parameters 
>> grows, this can become more and more complex since the various 
>> possibilities have to be known at compile time. 
>> >> >>> 3- Use generic interface objects for common elements (for example 
>> a viscosity class to calculate viscosity, etc.) and use inheritance to 
>> specialise the model. I think this could be inefficient because of the need 
>> to extensively look up through the vtable. 
>> >> >>> 
>> >> >>> What is the general consensus in the deal.II community on this 
>> problem? What's the best way to deal with multiple possibilities in the 
>> same assembly function? I would be very interested in hearing the 
>> perspective of the ASPECT author since it has such a large range of sub 
>> models. 
>> >> >>> 
>> >> >>> 
>> >> > -- 
>> >> > 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=

Re: [deal.II] Re: Mean Value Constraints

[Corbin Foucart]

Thank you very much for the pointers, Daniel!

- Including the templates header did the trick. I no longer need to copy 
the data over
- I was able to implement the penalty method and the DoF constraint method; 
the subspace projection method seems unambiguously best in terms of 
accuracy, iterations, and wall clock time. I benchmarked all three 
approaches against an MS for the pure Neumann problem.
- Re: spectrum--no, I only need the max and min eigenvalues to get at the 
conditioning of the system. However, this isn't crucial as I can just 
measure how the number of conjugate gradient iterations scale as a proxy.

Cheers,
Corbin

On Friday, May 19, 2023 at 1:47:42 PM UTC-4 d.arnd...@gmail.com wrote:

> Corbin,
>
> - Using ARPACK might indeed be an option. Do you really need to know the 
> whole spectrum, though?
> - Yes, you would just call AffineConstraints::add_line() to add constrain 
> a given DoF to zero. Note that it's a no-op if the DoF is already 
> constrained. This doesn't change the size of the linear system.
> - You will need to include deal.II/lac/affine_constraints.templates.h, 
> deal.II/lac/affine_constraints.h only contains the declarations but not the 
> definitions.
>
> Best,
> Daniel
>
> On Fri, May 19, 2023 at 12:04 PM Corbin Foucart  
> wrote:
>
>> I'm working on a similar pure Neumann problem with a rank-1 deficiency 
>> (pressure known only up to a constant). I've implemented the subspace 
>> projection method Konrad mentioned, with good results. However, I'm 
>> interested in comparing it to the single DoF constraint method, as well as 
>> an alternative penalization method alluded to in the paper by Bochev. I had 
>> a few questions:
>>
>>1. For a Sparse or ChunkSparse matrix object, is there a simple way 
>>of obtaining the eigenvalues? I know I can rebuild deal.ii with ARPACK, 
>> but 
>>this seems like overkill 
>>2. Exactly how would I use an AffineConstraints object to "pin" the 
>>degree of freedom? Is it a matter of a single call to 
>>AffineConstraints::add_line()? It's unclear to me from the documentation 
>>whether this would add an additional "row" to the linear system, or 
>> replace 
>>an existing one (I'd be looking to do the second)
>>3. I've implemented the subspace projection by wrapping the 
>>ChunkSparse matrix class  and overriding the vmult() function; however, 
>> at 
>>the moment, I'm forced to copy an original ChunkSparse matrix object into 
>>the wrapped object because the 
>>AffineConstraints::distribute_local_to_global() function doesn't 
>> recognize 
>>the wrapped matrix --I receive a linker error
>>
>> libHDG_experimental.so: error: undefined reference to 'void 
>> dealii::AffineConstraints::distribute_local_to_global>  
>> dealii::Vector >(dealii::FullMatrix const&, dealii::Vector const&, 
>> std::vector const&, 
>> HDG::LinearSystem::ChunkSparseWrapper&, dealii::Vector&, bool, 
>> std::integral_constant) const'
>>
>> I find that surprising, given that my wrapper class inherits 
>> publicly from ChunkSparseMatrix.
>>
>> Any guidance would be appreciated!
>> Corbin
>>
>> On Monday, December 28, 2020 at 1:23:19 PM UTC-5 Wolfgang Bangerth wrote:
>>
>>> On 12/26/20 3:06 AM, Konrad Simon wrote: 
>>> > What one can also do is just constrain one DoF to a specific value 
>>> (this would 
>>> > also remove rigid motion in elasticity). But think about your solution 
>>> > variable: If it is in the Sobolev space H^1 then point evaluations may 
>>> not be 
>>> > defined for dimension larger than 2. Similarly if, for example, the 
>>> pressure 
>>> > in a mechanical or fluid problem is often just in L^2. Point 
>>> evaluations do 
>>> > not make sense there at all. 
>>>
>>> Right, this is the correct approach: Constrain a single degree of 
>>> freedom to 
>>> zero (or any other value you choose) and solve the problem. Then you can 
>>> subtract the mean value of the solution *after* solving the linear 
>>> system. 
>>> (See VectorTools::subtract_mean_value and 
>>> VectorTools::compute_mean_value.) 
>>>
>>> If you're uncomfortable with the ill-posedness of taking a point value, 
>>> you 
>>> can also take the mean value along a small segment of the boundary 
>>> (step-1) or 
>>> a small part of the domain. But in practice, this is not necessary a

Re: [deal.II] Re: Mean Value Constraints

[Corbin Foucart]

I'm working on a similar pure Neumann problem with a rank-1 deficiency 
(pressure known only up to a constant). I've implemented the subspace 
projection method Konrad mentioned, with good results. However, I'm 
interested in comparing it to the single DoF constraint method, as well as 
an alternative penalization method alluded to in the paper by Bochev. I had 
a few questions:

   1. For a Sparse or ChunkSparse matrix object, is there a simple way of 
   obtaining the eigenvalues? I know I can rebuild deal.ii with ARPACK, but 
   this seems like overkill 
   2. Exactly how would I use an AffineConstraints object to "pin" the 
   degree of freedom? Is it a matter of a single call to 
   AffineConstraints::add_line()? It's unclear to me from the documentation 
   whether this would add an additional "row" to the linear system, or replace 
   an existing one (I'd be looking to do the second)
   3. I've implemented the subspace projection by wrapping the ChunkSparse 
   matrix class  and overriding the vmult() function; however, at the moment, 
   I'm forced to copy an original ChunkSparse matrix object into the wrapped 
   object because the AffineConstraints::distribute_local_to_global() function 
   doesn't recognize the wrapped matrix --I receive a linker error

libHDG_experimental.so: error: undefined reference to 'void 
dealii::AffineConstraints::distribute_local_to_global(dealii::FullMatrix const&, dealii::Vector const&, 
std::vector const&, 
HDG::LinearSystem::ChunkSparseWrapper&, dealii::Vector&, bool, 
std::integral_constant) const'

I find that surprising, given that my wrapper class inherits publicly from 
ChunkSparseMatrix.

Any guidance would be appreciated!
Corbin

On Monday, December 28, 2020 at 1:23:19 PM UTC-5 Wolfgang Bangerth wrote:

> On 12/26/20 3:06 AM, Konrad Simon wrote:
> > What one can also do is just constrain one DoF to a specific value (this 
> would 
> > also remove rigid motion in elasticity). But think about your solution 
> > variable: If it is in the Sobolev space H^1 then point evaluations may 
> not be 
> > defined for dimension larger than 2. Similarly if, for example, the 
> pressure 
> > in a mechanical or fluid problem is often just in L^2. Point evaluations 
> do 
> > not make sense there at all.
>
> Right, this is the correct approach: Constrain a single degree of freedom 
> to 
> zero (or any other value you choose) and solve the problem. Then you can 
> subtract the mean value of the solution *after* solving the linear system. 
> (See VectorTools::subtract_mean_value and VectorTools::compute_mean_value.)
>
> If you're uncomfortable with the ill-posedness of taking a point value, 
> you 
> can also take the mean value along a small segment of the boundary 
> (step-1) or 
> a small part of the domain. But in practice, this is not necessary and 
> Konrad's solution is what everyone seems to be doing.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/ae39fbcc-c76b-4c03-bfbb-93d19ab2e3d7n%40googlegroups.com.

Re: [deal.II] Understanding MeshWorker::mesh_loop order with adaptive refinement

[Corbin Foucart]

Hi Timo, I replied to this from a different email address, and it doesn't 
seem to have made it to the mailing list; most likely an error on my end. 
Here's what I wrote a few weeks ago:

> Do you mean a pull request to the documentation? Or a pull request for 
mesh_loop? I ask because I'm unsure if what I'm seeing with 
AssembleFlags::assemble_own_interior_faces_both is a bug or not.

> I've since re-written my code to assemble all the face data using 
AssembleFlags::assemble_own_interior_faces_once, then computing and 
applying the residuals in a separate mesh_loop call, and I'm once again 
recovering the correct residual in the global refinement and adaptive 
refinement cases; no unexpected behavior there.

Thank you,
Corbin 

On Friday, April 7, 2023 at 2:49:48 PM UTC-4 Timo Heister wrote:

> > to my eye, in both cases, the work on the face is done, followed by the 
> work on the cell. I think, however, the default behavior is to work on the 
> cells first, followed by work on the faces.
>
> That is correct. Would you be able to submit a pull request to fix this?
>
> On Thu, Apr 6, 2023 at 8:30 AM Corbin Foucart  
> wrote:
> >
> > I also think there may be a small typo in the documentation; "If the 
> flag AssembleFlags: : assemble_own_cells is passed, then the default 
> behavior is to first loop over faces and do the work there, and then 
> compute the actual work on the
> > ZjQcmQRYFpfptBannerStart
> > This Message Is From an External Sender
> > Use caution when opening links or attachments if you do not recognize 
> the sender.
> >
> > ZjQcmQRYFpfptBannerEnd
> > I also think there may be a small typo in the documentation;
> >
> > "If the flag AssembleFlags::assemble_own_cells is passed, then the 
> default behavior is to first loop over faces and do the work there, and 
> then compute the actual work on the cell. It is possible to perform the 
> integration on the cells after working on faces, by adding the flag 
> AssembleFlags::cells_after_faces."
> >
> > to my eye, in both cases, the work on the face is done, followed by the 
> work on the cell. I think, however, the default behavior is to work on the 
> cells first, followed by work on the faces.
> >
> > On Mon, Apr 3, 2023 at 6:35 PM Corbin Foucart  
> wrote:
> >>
> >> Hello everyone,
> >>
> >> I'm solving a 1D explicit DG-FEM problem and I've encountered behavior 
> that I don't understand using MeshWorker::mesh_loop.
> >>
> >> I'm using the MeshWorker::assemble_own_interior_faces_both flag since I 
> want to do work on each face corresponding to the same cell interior
> >> If the grid is created using typical GridGenerator calls, the order is 
> exactly as I'd expect: first, work is done on the cell, then the faces 
> (boundary or interior), followed by the next cell (further, all mass 
> matrices are the same)
> >> However, if I manually adapt the mesh by refining some cells, the order 
> seems to change; the face work is done without respect to the cell interior 
> last worked on
> >>
> >> I've attached a stripped-down program that illustrates this behavior on 
> a toy mesh in 1D, as well as the output.
> >>
> >> Ultimately, my goal is to assemble an inverse mass matrix on each cell, 
> and apply it to a residual vector containing interior and face 
> contributions (which can be done element-wise since the elements are FE_DGQ)
> >> I was attempting to store the inverse via CopyData and then apply it in 
> the copy worker.
> >> However, I'm finding that due to the order of execution, I can't rely 
> on the face work being done immediately after the cell work, and the 
> inverse mass matrix stored to CopyData is often from another cell than the 
> faces being worked on.
> >> I could do the assembly and application of the inverse mass matrices 
> separately, or store the inverse mass matrices in a map to the cell 
> iterators, but I'm curious why this ordering occurs.
> >>
> >> Am I misunderstanding how mesh_loop is supposed to work? Any guidance 
> would be appreciated!
> >>
> >> Corbin
> >>
> >> --
> >> 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 v

Re: [deal.II] Converting/viewing Vector objects as Tensor<1, dim> objects

[Corbin Foucart]

Thanks for the fast reply, Wolfgang!

I have done some preliminary benchmarking of the code and indeed this 
operation is not anywhere close to a bottleneck--I was asking primarily out 
of curiosity, as I'm refactoring the code and was considering what return 
type of an abstract representation of the flux makes sense in terms of 
class design. In light of your answer it seems that a view-based design 
doesn't make a lot of sense, and I might as well return a Tensor object 
directly, as that's the easiest for user code to work with.

Best,
Corbin

On Monday, May 1, 2023 at 7:20:05 PM UTC-4 Wolfgang Bangerth wrote:

> On 5/1/23 16:15, Corbin Foucart wrote:
> > 
> > Is there an developer-intended way to view a Vector object with 
> > dim entries as a Tensor<1, dim> object for ease of multiplication with 
> > other Tensor<1, dim> instances such as normal vectors?
> > I can do something like:
> > 
> > // dim = 3
> > Vector vec({4,5,6});
> > ArrayView view(&a[0], 3);
> > Tensor<1, 3, double> tensor(view);
> > 
> > but from the documentation this seems to make a copy, and I wondered if 
> > there was a deal.ii way to "view" the Vector as a Tensor object. 
> > This operation is the result of a numerical flux computation and is 
> > called many times during assembly.
>
> Correct, this makes a copy. There really isn't a different way (in 
> deal.II): Tensor objects own their data, they are not "views", and so 
> they cannot avoid the copy. Of course, if you don't want to do the copy, 
> then you can always write out the tensor operation you are interested in 
> from the original vector elements by hand.
>
> That said: It would greatly surprise me if the copying of elements is 
> limiting the speed of your program. Have you benchmarked this? Are you 
> sure that you aren't trying to optimize a part of the program that 
> doesn't actually need to be optimized?
>
> Best
> W.
>
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/f0ef3b0e-34c0-45bf-aaff-0f0b0692266fn%40googlegroups.com.

[deal.II] Converting/viewing Vector objects as Tensor<1, dim> objects

[Corbin Foucart]

Hello everyone,

Is there an developer-intended way to view a Vector object with dim 
entries as a Tensor<1, dim> object for ease of multiplication with other 
Tensor<1, dim> instances such as normal vectors?
I can do something like:

// dim = 3
Vector vec({4,5,6});
ArrayView view(&a[0], 3); 
Tensor<1, 3, double> tensor(view);

but from the documentation this seems to make a copy, and I wondered if 
there was a deal.ii way to "view" the Vector as a Tensor object. 
This operation is the result of a numerical flux computation and is called 
many times during assembly. 

Thank you,
Corbin

-- 
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/722fb438-322b-4b40-ac21-d76c9cf749d4n%40googlegroups.com.

Re: [deal.II] Understanding MeshWorker::mesh_loop order with adaptive refinement

[Corbin Foucart]

I also think there may be a small typo in the documentation
<https://www.dealii.org/current/doxygen/deal.II/group__MeshWorker.html#ga76ec61fbd188fb320fe8ca166a79b322>
;

"If the flag AssembleFlags::assemble_own_cells is passed, then the default
behavior is to first loop over faces and do the work there, and then
compute the actual work on the cell. It is possible to perform the
integration on the cells after working on faces, by adding the flag
AssembleFlags::cells_after_faces."

to my eye, in both cases, the work on the face is done, followed by the
work on the cell. I think, however, the default behavior is to work on the
cells first, followed by work on the faces.

On Mon, Apr 3, 2023 at 6:35 PM Corbin Foucart 
wrote:

> Hello everyone,
>
> I'm solving a 1D explicit DG-FEM problem and I've encountered behavior
> that I don't understand using MeshWorker::mesh_loop.
>
>- I'm using the MeshWorker::assemble_own_interior_faces_both flag
>since I want to do work on each face corresponding to the same cell 
> interior
>- If the grid is created using typical GridGenerator calls, the order
>is exactly as I'd expect: first, work is done on the cell, then the faces
>(boundary or interior), followed by the next cell (further, all mass
>matrices are the same)
>- However, if I manually adapt the mesh by refining some cells, the
>order seems to change; the face work is done without respect to the cell
>interior last worked on
>
> I've attached a stripped-down program that illustrates this behavior on a
> toy mesh in 1D, as well as the output.
>
>- Ultimately, my goal is to assemble an inverse mass matrix on each
>cell, and apply it to a residual vector containing interior and face
>contributions (which can be done element-wise since the elements are
>FE_DGQ)
>- I was attempting to store the inverse via CopyData and then apply it
>in the copy worker.
>- However, I'm finding that due to the order of execution, I can't
>rely on the face work being done immediately after the cell work, and the
>inverse mass matrix stored to CopyData is often from another cell than the
>faces being worked on.
>- I could do the assembly and application of the inverse mass matrices
>separately, or store the inverse mass matrices in a map to the cell
>iterators, but I'm curious why this ordering occurs.
>
> Am I misunderstanding how mesh_loop is supposed to work? Any guidance
> would be appreciated!
>
> Corbin
>
> --
> 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/ac4d4550-b014-4440-b2c2-227dba71fffen%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/ac4d4550-b014-4440-b2c2-227dba71fffen%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/CAE%2BMPodXQ5tZWiOkhq1j-p4uuvQiWOg3n57R2z7cWCKMz1pKFQ%40mail.gmail.com.

[deal.II] Understanding MeshWorker::mesh_loop order with adaptive refinement

[Corbin Foucart]

Hello everyone, 

I'm solving a 1D explicit DG-FEM problem and I've encountered behavior that 
I don't understand using MeshWorker::mesh_loop.

   - I'm using the MeshWorker::assemble_own_interior_faces_both flag since 
   I want to do work on each face corresponding to the same cell interior
   - If the grid is created using typical GridGenerator calls, the order is 
   exactly as I'd expect: first, work is done on the cell, then the faces 
   (boundary or interior), followed by the next cell (further, all mass 
   matrices are the same)
   - However, if I manually adapt the mesh by refining some cells, the 
   order seems to change; the face work is done without respect to the cell 
   interior last worked on

I've attached a stripped-down program that illustrates this behavior on a 
toy mesh in 1D, as well as the output.

   - Ultimately, my goal is to assemble an inverse mass matrix on each 
   cell, and apply it to a residual vector containing interior and face 
   contributions (which can be done element-wise since the elements are 
   FE_DGQ) 
   - I was attempting to store the inverse via CopyData and then apply it 
   in the copy worker. 
   - However, I'm finding that due to the order of execution, I can't rely 
   on the face work being done immediately after the cell work, and the 
   inverse mass matrix stored to CopyData is often from another cell than the 
   faces being worked on. 
   - I could do the assembly and application of the inverse mass matrices 
   separately, or store the inverse mass matrices in a map to the cell 
   iterators, but I'm curious why this ordering occurs.

Am I misunderstanding how mesh_loop is supposed to work? Any guidance would 
be appreciated!

Corbin

-- 
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/ac4d4550-b014-4440-b2c2-227dba71fffen%40googlegroups.com.
// mesh_loop_order.cc

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

using namespace dealii;

struct CopyData
{
  FullMatrix  inverse_mass_matrix;
  Vector  residual_contribution;
  std::vectorlocal_dof_indices;

  template
  void reinit(const Iterator &cell, unsigned int dofs_per_cell)
  {
inverse_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
residual_contribution.reinit(dofs_per_cell);
local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
  }
};


template
class Wrapper
{
  public:
Wrapper(DoFHandler&  dof_handler,
FE_DGQ&  fe)
  : dofh(&dof_handler)
  , fe(&fe)
  , cell_quadrature(fe.tensor_degree() + 1)
  , face_quadrature(fe.tensor_degree() + 1)
{ }

void do_mesh_loop();

  protected:
const QGauss cell_quadrature;
const QGauss   face_quadrature;
SmartPointer> dofh;
const SmartPointer> fe;
};


template
void Wrapper::do_mesh_loop()
{
  using Iterator = typename DoFHandler::active_cell_iterator;
  using ScratchData = MeshWorker::ScratchData;

  const auto cell_worker = [&](const Iterator&cell,
   ScratchData   &sd,
   CopyData  ©_data) {

std::cout << "cell worker entered on cell " << cell << std::endl;
const FEValues &fev = sd.reinit(cell);
const unsigned int n_cell_dofs = fev.get_fe().n_dofs_per_cell();

std::cout << "\tcopy data reinitialized" << std::endl;
copy_data.reinit(cell, n_cell_dofs);
// cell work here (assemble inverse mass matrix)
  };

  const auto face_worker = [&](const Iterator  &cell,
   const unsigned int  &f,
   const unsigned int  &sf,
   const Iterator  &ncell,
   const unsigned int  &nf,
   const unsigned int  &nsf,
   ScratchData &sd,
   CopyData©_data) {

std::cout << "face worker called on cell " << cell << std::endl;
const FEInterfaceValues &fe_iv = sd.reinit(cell, f, sf, ncell, nf, nsf);
  };

  const auto boundary_worker = [&](const Iterator  &cell,
   const unsigned int  &face_no,
   ScratchData &sd,
   CopyData©_data) {

std::cout << "boundary worker entered on cell " << cell << std::endl;
const FEFaceValues &fe_fv = sd.reinit(cell, face_no);
  };

  AffineConstraints constraints;
  constraints.close();

  const auto copier = [&](const CopyDa

[deal.II] Factorization of ChunkSparseMatrix

[Corbin Foucart]

Hello everyone,

I'm looking to LU factorize a ChunkSparseMatrix instance (the same 
type as is used in step-51). However, it seems that ChunkSparseMatrix does 
not inherit from SparseMatrix, so I receive the predictable:

undefined reference to 'void 
dealii::SparseDirectUMFPACK::factorize 
>(dealii::ChunkSparseMatrix const&)'

I have a few ideas around this, but I'm not sure which is simplest to do 
(or whether I'm missing an easy templating/casting trick that would avoid 
the problem):

   1. Implement a way to copy the ChunkSparseMatrix to a SparseMatrix as an 
   intermediate before factorization
   2. Implement SparseDirectUMFPACK::factorize(const Matrix &matrix) 
   directly (this is perhaps tougher, since the function is templated by 
   Matrix, and I can't seem to find the explicit template instantiations in 
   the source code)
   3. Do away with use of ChunkSparseMatrix entirely, and simply assemble 
   into a SparseMatrix throughout (since I'm factorizing it anyway--perhaps 
   this is the truest to the intentions of the library)


Any guidance would be appreciated!


-- 
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/6636e5b6-a09a-4bc2-8ec1-7bae018adfdcn%40googlegroups.com.

[deal.II] Adding two finite element fields of unequal order

[Corbin Foucart]

Hello everyone,

   - I'm looking to add a finite element solution of one polynomial order 
   to another (say, adding a p=2 field to a p=3 field, without loss of 
   generality)
   - I can do this by hand the usual way
  - declare a nodal quadrature rule (at p=3)
  - iterate over the active cells in DoFHandler, using 
  FEValues::get_function_values to add the results appropriately, either 
  manually or with MeshWorker
   - However, I suspect that deal.ii must have a more idiomatic way of 
   doing so, since this must be a common operation

Is there a library function to do operations like this? Any guidance would 
be appreciated.

Thank you,
Corbin

-- 
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/48d14cd4-a252-4332-bac3-7fd4f9dec1c5n%40googlegroups.com.

[deal.II] Iterating over mesh cells in a custom order

[Corbin Foucart]

Hello everyone,

I have a problem in which I'm propagating information downwards in depth by 
solving the same local finite element problem on each element in an 
adaptive grid. The only condition is that the cells above the current cell 
must have already been worked on. 

I'm looking for a way to loop over cells using something akin to 
MeshWorker::mesh_loop or WorkStream::run, but with a custom order, namely, 
according to the z-coordinate of each cell center. Is there a simple way to 
do this?

If not, I can always loop over the cells manually by recording a 
depth-sorted list, but then I'd lose the multi-thread capabilities of 
mesh_loop() and the like. Any advice would be appreciated.

Thank you,
Corbin

-- 
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/a599de3d-86b8-4681-8204-bb0881b4a35bn%40googlegroups.com.

[deal.II] Using MatrixFree functionality without distributed arrays

[Corbin Foucart]

Hello everyone,

I've been reading through some of the matrix-free tutorials---step-67 in 
particular. I'd like to write something similar, an RKDG solver for 
hyperbolic problems for which I can use explicit RK time integration.

   - Does it make sense to use the MatrixFree functionality with regular 
   Vector operands? The tutorials often make use of 
   LinearAlgebra::distributed::Vector objects to store the vectors.
   - I don't want to use the latter because I'm driving the code from 
   Python, and any use of MPI or distributed memory could be a headache in 
   communicating data back and forth between the two languages. Additionally, 
   I'm solving many small problems as opposed to a single large one.
   - Ultimately, I'm looking for a fast way to apply the elemental inverse 
   mass matrices and face integrations, so MatrixFree functionality seemed the 
   logical choice.

Best, 
Corbin

-- 
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/09075ee8-93e7-412a-9b86-8a220785ea2fn%40googlegroups.com.

Re: [deal.II] Solving separate FEM problems on subdomains

[Corbin Foucart]

I never replied to this, but the approach of using assembling only the 
interface terms worked very well, and was painless to code. 

Thanks for the great idea, this saved me a ton of time :-)

Corbin

On Wednesday, March 23, 2022 at 11:17:15 PM UTC-4 Corbin Foucart wrote:

> Bruno and Wolfgang,
>
> Thank you for the helpful suggestions--I had not thought of the 
> block-diagonal linear system approach in particular, but I think it is a 
> clever way to do what I had previously planned. I don't anticipate it being 
> too expensive, even unordered, but if so, I will implement one of the other 
> two approaches.
>
> Corbin
>
> On Tuesday, March 22, 2022 at 10:38:32 PM UTC-4 Wolfgang Bangerth wrote:
>
>> On 3/21/22 16:47, Corbin Foucart wrote:
>> > 
>> > I have an extruded mesh, and I would like to solve a finite element 
>> problem 
>> > separately on each 'column' as part of a larger finite element solve in 
>> which 
>> > all DoF are coupled through a discontinuous Galerkin scheme for an 
>> elliptic 
>> > equation--so multiply defined values on the boundaries of each column 
>> are not 
>> > a problem.
>> > 
>> > * I already have the columns tagged by material ID
>> > * I imagine I should instantiate a list of new Triangulations, 
>> DoFHandlers,
>> > and solver classes for each column and then map the solutions back to 
>> the
>> > original domain
>> > * *Is there a good way to extract a separate triangulation for each 
>> column*?
>> > I was looking for something similar to extract_boundary_mesh but with a
>> > material_id
>> > 
>> > Alternatively, it seems that I could assign the DoFs in the globally 
>> coupled 
>> > DoFHandler with DoFTools::get_subdomain_association(). *Would this be a 
>> better 
>> > way to set up each finite element problem by column?* It's unclear to 
>> me if I 
>> > could use that information to set up an independent linear system for 
>> each 
>> > subdomain.
>>
>> Corbyn,
>> you could go the way of creating triangulations, DoFHandlers, etc, for 
>> each of 
>> the columns. But you might also want to consider one of the following two 
>> approaches:
>>
>> * Since you're already doing a DG approach, consider building your matrix 
>> with 
>> only the interface terms that correspond to the "bottom" and "top" of 
>> each of 
>> your cells. This means that you don't build any term that corresponds to 
>> coupling columns to each other. Solving such a linear system then 
>> corresponds 
>> to solving all of the 1d problems at the same time. (I will note that 
>> such a 
>> matrix would also make for a good preconditioner; it corresponds to one 
>> of the 
>> directions one would solve for in the ADI method.)
>>
>> Conceptually, the matrix you obtain using the approach above is block 
>> diagonal 
>> where each block corresponds to one column. In practice, unless you do 
>> something special, the matrix will not actually look block diagonal 
>> because 
>> the order of DoFs does not respect the column-by-column order. But it 
>> wouldn't 
>> be very difficult to implement a function that renumbers DoFs in such a 
>> way 
>> that they are numbered column by column, and then you would indeed have a 
>> block diagonal matrix.
>>
>> * If you don't like this approach, you can create a mapping from the DoFs 
>> of 
>> each column to a column-local numbering from zero to n_dofs_per_column. 
>> Instead of writing into global vectors and matrices in the method above 
>> where 
>> you omit the lateral face terms, you can translate writing into smaller 
>> matrices and vectors that correspond to just one column. You can then 
>> solve 
>> the resulting linear system. The mapping of DoF indices then also allows 
>> you 
>> to get back at what the global index of these DoFs was.
>>
>> Best
>> W.
>>
>> -- 
>> 
>> Wolfgang Bangerth email: bang...@colostate.edu
>> www: http://www.math.colostate.edu/~bangerth/
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/5752745c-8e3d-4724-b835-2c362ef89d59n%40googlegroups.com.

[deal.II] Accessing and modifying modal properties of a FE_DGQ solution

[Corbin Foucart]

Hello everyone, 

I have a discontinuous velocity field composed of FE_DGQ elements; I would 
like to access a modal representation of the solution coefficients on each 
element and modify them (e.g., to perform slope limiting). 

I was planning on 

   - manually constructing the Vandermonde matrix on each element to 
   convert the solution representation to a modal one
   - modifying the solution coefficients
   - transforming them back to a nodal representation
   - updating the entries in the global solution vector
   
Is there a better way of doing this? I want to make sure I'm not 
reinventing the wheel in the event that this functionality is already 
implemented. 

I was looking at the documentation for FE_DGP but it seems that these 
elements may not span the same space as their corresponding nodal 
counterparts of the same polynomial order, and the two finite element 
objects are not necessarily meant to be used in conjunction (for example, 
by defining two FEValues objects with FE_DGP and FE_DGQ finite elements and 
a nodal quadrature rule).

Thank you,
Corbin

-- 
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/fe346eae-8f32-4953-9bff-ac1806016edcn%40googlegroups.com.

Re: [deal.II] Re: Python bindings

[Corbin Foucart]

Bruno, 

Ok, I've added my fix--please see pr13631 
<https://github.com/dealii/dealii/pull/13631>.

Best,
Corbin

On Thursday, April 21, 2022 at 9:04:43 AM UTC-4 bruno.t...@gmail.com wrote:

> Corbin,
>
> Le mer. 20 avr. 2022 à 19:10, Corbin Foucart  a 
> écrit :
>
>>
>> I'm hesitant to propose a change, as doing so may break the build process 
>> for other people using different Python binding configurations. I'm also by 
>> no means a CMake expert :-) However, if it's a safe assumption that anyone 
>> using the bindings is using Python3, sure.
>>
> Python2 reached end-of-life two years ago, I think it's fair to require 
> Python3.
>  
>
>> As for what type of bindings, I was planing on writing some for the 
>> mesh-refinement indicators and executing mesh refinement. Most of the 
>> bindings I'm planning to write would be specific to solvers I've developed 
>> on top of deal ii, but some of the functionality might be useful to the 
>> community. The idea is to drive the adaptive refinement of a solver with a 
>> Python code that has access to machine learning libraries such as PyTorch.
>>
> That sounds very interesting. Note that you can already flag a cell for 
> refinement and then execute the refinement (see here 
> <https://github.com/dealii/dealii/blob/master/contrib/python-bindings/notebooks/tutorial-1.ipynb>
> )
>
> Best,
>
> Bruno
>  
>

-- 
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/0c973efa-4041-43c1-af63-45ac8beb987bn%40googlegroups.com.

Re: [deal.II] Re: Interacting with Python data // external codes

[Corbin Foucart]

I was able to get a minimal deal ii wrapper example working with Cython as 
suggested above; I'm including the code here in case others find it useful.

   - This wraps some of the basic Triangulation functionality, inspired by 
   step-1
   - The CMake build creates a standalone shared object file which can be 
   imported by Python
   - All the functionality wrapped can be driven entirely from Python in a 
   Conda environment

The difficult part was the CMake configuration; it seems somewhat 
non-portable (due to the particulars of the conda environment), and in 
order to talk to Cython, I had to include the CMake modules for cython 
discovery 
<https://github.com/thewtex/cython-cmake-example/tree/master/cmake> 
provided by Kitware. However, once configured, I agree with Alex that this 
seems like a fairly painless way to drive deal.ii from Python. I'll post 
again if there are issues communicating data back-and-forth between the two 
languages.

Corbin
On Wednesday, April 20, 2022 at 6:59:25 PM UTC-4 Corbin Foucart wrote:

> Hi Alex,
>
> I've done the communication between the two codes with IO transfer, and 
> I'd like to explore the Cython idea. How are you able to build the Cython 
> wrapper? 
>
>- In order to wrap a class making use of, say, #include 
> , I can generate a shared object file using CMake
>- however, the setup.py file for the Cython compilation knows nothing 
>about the complicated Deal II build process and can't resolve the 
>dependencies
>- Specifying them directly without CMake would be like compiling a 
>deal.ii tutorial program by hand with gcc
>
> Were you able to build the Cython wrapper in CMake directly? That seems 
> like it would be much better
>
> Thank you,
> Corbin
>
> On Thursday, December 3, 2020 at 3:22:12 AM UTC-5 Alex Cobb wrote:
>
>> Hi Corbin,
>>
>> On Wed, 2 Dec 2020 at 14:55, Corbin Foucart  wrote:
>>
>>> The code will have to be driven in Python, but hopefully I can find a 
>>> good way to access the deal.ii data structures without writing to file. I'm 
>>> using a time-dependent HDG scheme on an adaptive grid. From a python 
>>> session, is it possible to access:
>>>
>>>1. The state of the triangulation ( this seems doable from the 
>>>python bindings, as in tutorial 1 in the notebooks Bruno linked)
>>>2. For every cell in the mesh
>>>   1. the nodal DG coefficients u_h
>>>   2. the solution jumps [[u_h]] on all faces of the cell (between 
>>>   the solution on this cell and its neighbors)
>>>   3. the right hand side (source term and time-dependent forcing), 
>>>   ideally at the nodal points of the cell
>>>3. The coarsen or refine flags for the triangulation (seems possible 
>>>as in 1)?
>>>
>>> I think the best way to do this would be with the Cython approach 
>>> suggested by Alex, however, it's unclear to me the best way to pass the 
>>> cell-based data to Python. Is it possible to define a struct on from 
>>> deal.ii containing the data in (2) and pass an array of these structs (one 
>>> per cell) to the python session? 
>>>
>>
>> You always have the option of retrieving primitives as their Python 
>> equivalents (e,g., double -> Python floats) via some kind of getter 
>> function; this means that you will have the overhead of a Python function 
>> call and creation of an object on the heap for each such retrieval (from 
>> Python).  If you have large amounts of data to pass back and forth via 
>> memory, Numpy arrays are a good option (you could also use the Python 
>> Standard Library struct module, but I can't see any advantages to this). 
>> Cython can give you the pointer to the block of memory where a Numpy 
>> array's data reside, and you can use that to read and write data to 
>> contiguous arrays.  If you want to store data in structs (or the data you 
>> want are already packed in structs) you could store the data in a Numpy 
>> record array (array of structs).  Your preferred strategy will probably 
>> depend on how exactly you will use the data on the Python side.
>>
>> Once you start passing pointers to memory around, an important 
>> consideration is always who owns the memory and is responsible for 
>> deallocating it. The strategy I have mostly used (and been happy with) is 
>> to allocate all the memory I need from the Python side, either by creating 
>> Numpy arrays or by calling PyMem_Malloc in the constructor and PyMem_Free 
>> in the destructor of a Cython extension class (so that the memory is tied 
>&

Re: [deal.II] Re: Python bindings

[Corbin Foucart]

Hi Bruno,

I'm hesitant to propose a change, as doing so may break the build process 
for other people using different Python binding configurations. I'm also by 
no means a CMake expert :-) However, if it's a safe assumption that anyone 
using the bindings is using Python3, sure.

As for what type of bindings, I was planing on writing some for the 
mesh-refinement indicators and executing mesh refinement. Most of the 
bindings I'm planning to write would be specific to solvers I've developed 
on top of deal ii, but some of the functionality might be useful to the 
community. The idea is to drive the adaptive refinement of a solver with a 
Python code that has access to machine learning libraries such as PyTorch.

Corbin

On Wednesday, April 20, 2022 at 8:15:57 AM UTC-4 bruno.t...@gmail.com wrote:

> Corbin,
>
> That's great. Let us know if you encounter any other problems. What kind 
> of bindings do you want to write? BTW, do you want to open a PR with the 
> FindPython3 change? 
>
> Best,
>
> Bruno
>
> Le mer. 20 avr. 2022 à 06:56, Corbin Foucart  a 
> écrit :
>
>> It looks like in the python-bindings CMakeLists.txt file here: 
>> https://github.com/dealii/dealii/blob/master/contrib/python-bindings/CMakeLists.txt
>>
>> There's a call to 
>> > INCLUDE(FindPythonInterp)
>> on line 24. 
>>
>> It seems that for cmake versions > 3.12, FindPythonInterp has been 
>> deprecated, with FindPython3 replacing it... modifying that line of the 
>> file manually, as well as manually setting the Boost version and python 
>> versions in the same CMakeLists.txt file to 1.71 and 3.8, respectively, I 
>> am able to get past the cmake errors. If the build process succeeds, I'll 
>> try writing some bindings.
>>
>> Corbin
>>
>> On Wednesday, April 20, 2022 at 6:06:21 AM UTC-4 Corbin Foucart wrote:
>>
>>> Hi Bruno,
>>>
>>> Thanks for the advice. 
>>>
>>> It's my understanding that compiling boost.python code outputs a shared 
>>> object file (in this case, I'm guessing something like PyDealII/Debug.so) 
>>> which can be loaded as a module by a python interpreter (in this case by 
>>> adding the shared object location to the PYTHONPATH). Running 'find' in my 
>>> build directory, I was unable to find the shared object file.
>>>
>>> Digging deeper, I found that the actual flag is 
>>> -DDEAL_II_COMPONENT_PYTHON_BINDING*S*=ON, so I hadn't built the library 
>>> with the correct configuration. Attempting to build the python bindings, 
>>> the build of deal.ii fails, with the following error
>>>
>>> ```
>>> -- Setting up python bindings
>>> -- Found PythonInterp: /usr/bin/python (found version "2.7.18") 
>>> -- Found PythonLibs: /usr/lib/x86_64-linux-gnu/libpython2.7.so (found 
>>> version "2.7.18") 
>>> CMake Error at 
>>> /usr/share/cmake-3.16/Modules/FindPackageHandleStandardArgs.cmake:146 
>>> (message):
>>>   Could NOT find Boost (missing: python27) (found suitable version 
>>> "1.71.0",
>>>   minimum required is "1.67")
>>> Call Stack (most recent call first):
>>>   /usr/share/cmake-3.16/Modules/FindPackageHandleStandardArgs.cmake:393 
>>> (_FPHSA_FAILURE_MESSAGE)
>>>   /usr/share/cmake-3.16/Modules/FindBoost.cmake:2179 
>>> (find_package_handle_standard_args)
>>>   contrib/python-bindings/CMakeLists.txt:46 (_FIND_PACKAGE)
>>> ```
>>>
>>> I'm confused by this--I've installed boost and the build process finds 
>>> it, but I don't know why it's looking for Python 2.7, I have more recent 
>>> system versions (namely, 3.8) in the same directory. I imagine this might 
>>> be causing the issue with Boost, and I wouldn't be developing bindings for 
>>> python 2 anyway, as it's not supported. I'm using an Ubuntu 18.04 machine 
>>> and attempting to install deal.ii-9.3.3. Am I building this the wrong way? 
>>>
>>> Thank you,
>>> Corbin
>>>
>>> On Tuesday, April 19, 2022 at 4:37:11 PM UTC-4 bruno.t...@gmail.com 
>>> wrote:
>>>
>>>> Corbin,
>>>>
>>>> On Tuesday, April 19, 2022 at 3:04:39 PM UTC-4 corbin@gmail.com 
>>>> wrote:
>>>>
>>>>> I've built deal.ii version 9.3.0 with the 
>>>>> -DDEAL_II_COMPONENT_PYTHON_BINDING=ON configuration, yet when I run the 
>>>>> notebook, the Jupyter kernel is unable to locate the P

Re: [deal.II] Re: Interacting with Python data // external codes

[Corbin Foucart]

Hi Alex,

I've done the communication between the two codes with IO transfer, and I'd 
like to explore the Cython idea. How are you able to build the Cython 
wrapper? 

   - In order to wrap a class making use of, say, #include 
, I can generate a shared object file using CMake
   - however, the setup.py file for the Cython compilation knows nothing 
   about the complicated Deal II build process and can't resolve the 
   dependencies
   - Specifying them directly without CMake would be like compiling a 
   deal.ii tutorial program by hand with gcc

Were you able to build the Cython wrapper in CMake directly? That seems 
like it would be much better

Thank you,
Corbin

On Thursday, December 3, 2020 at 3:22:12 AM UTC-5 Alex Cobb wrote:

> Hi Corbin,
>
> On Wed, 2 Dec 2020 at 14:55, Corbin Foucart  wrote:
>
>> The code will have to be driven in Python, but hopefully I can find a 
>> good way to access the deal.ii data structures without writing to file. I'm 
>> using a time-dependent HDG scheme on an adaptive grid. From a python 
>> session, is it possible to access:
>>
>>1. The state of the triangulation ( this seems doable from the python 
>>bindings, as in tutorial 1 in the notebooks Bruno linked)
>>2. For every cell in the mesh
>>   1. the nodal DG coefficients u_h
>>   2. the solution jumps [[u_h]] on all faces of the cell (between 
>>   the solution on this cell and its neighbors)
>>   3. the right hand side (source term and time-dependent forcing), 
>>   ideally at the nodal points of the cell
>>3. The coarsen or refine flags for the triangulation (seems possible 
>>as in 1)?
>>
>> I think the best way to do this would be with the Cython approach 
>> suggested by Alex, however, it's unclear to me the best way to pass the 
>> cell-based data to Python. Is it possible to define a struct on from 
>> deal.ii containing the data in (2) and pass an array of these structs (one 
>> per cell) to the python session? 
>>
>
> You always have the option of retrieving primitives as their Python 
> equivalents (e,g., double -> Python floats) via some kind of getter 
> function; this means that you will have the overhead of a Python function 
> call and creation of an object on the heap for each such retrieval (from 
> Python).  If you have large amounts of data to pass back and forth via 
> memory, Numpy arrays are a good option (you could also use the Python 
> Standard Library struct module, but I can't see any advantages to this). 
> Cython can give you the pointer to the block of memory where a Numpy 
> array's data reside, and you can use that to read and write data to 
> contiguous arrays.  If you want to store data in structs (or the data you 
> want are already packed in structs) you could store the data in a Numpy 
> record array (array of structs).  Your preferred strategy will probably 
> depend on how exactly you will use the data on the Python side.
>
> Once you start passing pointers to memory around, an important 
> consideration is always who owns the memory and is responsible for 
> deallocating it. The strategy I have mostly used (and been happy with) is 
> to allocate all the memory I need from the Python side, either by creating 
> Numpy arrays or by calling PyMem_Malloc in the constructor and PyMem_Free 
> in the destructor of a Cython extension class (so that the memory is tied 
> to the life cycle of the object).  Another option would be to allocate the 
> memory in a C++ class and wrap that with Cython.
>
> Good luck!
>
> Alex
>
> On Tuesday, November 24, 2020 at 11:02:21 PM UTC-5 Alex Cobb wrote:
>>
>>> On Wed, 25 Nov 2020 at 00:05, Bruno Turcksin  
>>> wrote:
>>>
>>>> deal.II has some limited support for python mainly for mesh 
>>>> manipulation. We have some python notebooks here 
>>>> <https://github.com/dealii/dealii/blob/master/contrib/python-bindings/notebooks/index.ipynb>.
>>>>  
>>>> I think what you want to do is similar to the step-62 notebook. Right now, 
>>>> the only way to interact with numpy is to print the data to a file and 
>>>> then 
>>>> load it (see here 
>>>> <https://dealii.org/current/doxygen/deal.II/classLinearAlgebra_1_1Vector.html#a2fadcc595d3e7e9ba44de8c85fd7595c>
>>>>  
>>>> and here 
>>>> <https://dealii.org/current/doxygen/deal.II/classSparseMatrix.html#a3141075e3ad6362fce005d2f1c8da699>).
>>>>  
>>>> If you want to manipulate the mesh directly in python, you need 
>>>> boost.python and you need to configure deal.II with 
>

[deal.II] Re: Python bindings

[Corbin Foucart]

It looks like in the python-bindings CMakeLists.txt file here: 
https://github.com/dealii/dealii/blob/master/contrib/python-bindings/CMakeLists.txt

There's a call to 
> INCLUDE(FindPythonInterp)
on line 24. 

It seems that for cmake versions > 3.12, FindPythonInterp has been 
deprecated, with FindPython3 replacing it... modifying that line of the 
file manually, as well as manually setting the Boost version and python 
versions in the same CMakeLists.txt file to 1.71 and 3.8, respectively, I 
am able to get past the cmake errors. If the build process succeeds, I'll 
try writing some bindings.

Corbin

On Wednesday, April 20, 2022 at 6:06:21 AM UTC-4 Corbin Foucart wrote:

> Hi Bruno,
>
> Thanks for the advice. 
>
> It's my understanding that compiling boost.python code outputs a shared 
> object file (in this case, I'm guessing something like PyDealII/Debug.so) 
> which can be loaded as a module by a python interpreter (in this case by 
> adding the shared object location to the PYTHONPATH). Running 'find' in my 
> build directory, I was unable to find the shared object file.
>
> Digging deeper, I found that the actual flag is 
> -DDEAL_II_COMPONENT_PYTHON_BINDING*S*=ON, so I hadn't built the library 
> with the correct configuration. Attempting to build the python bindings, 
> the build of deal.ii fails, with the following error
>
> ```
> -- Setting up python bindings
> -- Found PythonInterp: /usr/bin/python (found version "2.7.18") 
> -- Found PythonLibs: /usr/lib/x86_64-linux-gnu/libpython2.7.so (found 
> version "2.7.18") 
> CMake Error at 
> /usr/share/cmake-3.16/Modules/FindPackageHandleStandardArgs.cmake:146 
> (message):
>   Could NOT find Boost (missing: python27) (found suitable version 
> "1.71.0",
>   minimum required is "1.67")
> Call Stack (most recent call first):
>   /usr/share/cmake-3.16/Modules/FindPackageHandleStandardArgs.cmake:393 
> (_FPHSA_FAILURE_MESSAGE)
>   /usr/share/cmake-3.16/Modules/FindBoost.cmake:2179 
> (find_package_handle_standard_args)
>   contrib/python-bindings/CMakeLists.txt:46 (_FIND_PACKAGE)
> ```
>
> I'm confused by this--I've installed boost and the build process finds it, 
> but I don't know why it's looking for Python 2.7, I have more recent system 
> versions (namely, 3.8) in the same directory. I imagine this might be 
> causing the issue with Boost, and I wouldn't be developing bindings for 
> python 2 anyway, as it's not supported. I'm using an Ubuntu 18.04 machine 
> and attempting to install deal.ii-9.3.3. Am I building this the wrong way? 
>
> Thank you,
> Corbin
>
> On Tuesday, April 19, 2022 at 4:37:11 PM UTC-4 bruno.t...@gmail.com wrote:
>
>> Corbin,
>>
>> On Tuesday, April 19, 2022 at 3:04:39 PM UTC-4 corbin@gmail.com 
>> wrote:
>>
>>> I've built deal.ii version 9.3.0 with the 
>>> -DDEAL_II_COMPONENT_PYTHON_BINDING=ON configuration, yet when I run the 
>>> notebook, the Jupyter kernel is unable to locate the PyDealII module, that 
>>> is, the notebook cell
>>>
>>> > import PyDealII.Release as dealii
>>>
>>> fails with an import error.
>>>
>> I havent' use the python binding in a while but you probably need to set 
>> the PYTHONPATH to where pydealii is located. You may also need to add the 
>> path to libdealii.so to LD_LIBRARY_PATH
>>
>> Where is this library located, and how should I go about setting up a 
>>> (conda) environment to be able to use PyDealII? Are there dependencies?
>>>
>> The library should be located in the same directory where you installed 
>> deal.II. Look somewhere in /path/to/install/lib . I don't use conda so I 
>> can't help you with that. What kind of dependencies are you talking about? 
>> PyDealII only depends on deal.II and Boost.Python, nothing else
>>
>> Best,
>>
>> Bruno 
>>
>

-- 
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/deac2b5e-0449-40f6-aeaf-ff83e536e614n%40googlegroups.com.

[deal.II] Re: Python bindings

[Corbin Foucart]

Hi Bruno,

Thanks for the advice. 

It's my understanding that compiling boost.python code outputs a shared 
object file (in this case, I'm guessing something like PyDealII/Debug.so) 
which can be loaded as a module by a python interpreter (in this case by 
adding the shared object location to the PYTHONPATH). Running 'find' in my 
build directory, I was unable to find the shared object file.

Digging deeper, I found that the actual flag is 
-DDEAL_II_COMPONENT_PYTHON_BINDING*S*=ON, so I hadn't built the library 
with the correct configuration. Attempting to build the python bindings, 
the build of deal.ii fails, with the following error

```
-- Setting up python bindings
-- Found PythonInterp: /usr/bin/python (found version "2.7.18") 
-- Found PythonLibs: /usr/lib/x86_64-linux-gnu/libpython2.7.so (found 
version "2.7.18") 
CMake Error at 
/usr/share/cmake-3.16/Modules/FindPackageHandleStandardArgs.cmake:146 
(message):
  Could NOT find Boost (missing: python27) (found suitable version "1.71.0",
  minimum required is "1.67")
Call Stack (most recent call first):
  /usr/share/cmake-3.16/Modules/FindPackageHandleStandardArgs.cmake:393 
(_FPHSA_FAILURE_MESSAGE)
  /usr/share/cmake-3.16/Modules/FindBoost.cmake:2179 
(find_package_handle_standard_args)
  contrib/python-bindings/CMakeLists.txt:46 (_FIND_PACKAGE)
```

I'm confused by this--I've installed boost and the build process finds it, 
but I don't know why it's looking for Python 2.7, I have more recent system 
versions (namely, 3.8) in the same directory. I imagine this might be 
causing the issue with Boost, and I wouldn't be developing bindings for 
python 2 anyway, as it's not supported. I'm using an Ubuntu 18.04 machine 
and attempting to install deal.ii-9.3.3. Am I building this the wrong way? 

Thank you,
Corbin

On Tuesday, April 19, 2022 at 4:37:11 PM UTC-4 bruno.t...@gmail.com wrote:

> Corbin,
>
> On Tuesday, April 19, 2022 at 3:04:39 PM UTC-4 corbin@gmail.com wrote:
>
>> I've built deal.ii version 9.3.0 with the 
>> -DDEAL_II_COMPONENT_PYTHON_BINDING=ON configuration, yet when I run the 
>> notebook, the Jupyter kernel is unable to locate the PyDealII module, that 
>> is, the notebook cell
>>
>> > import PyDealII.Release as dealii
>>
>> fails with an import error.
>>
> I havent' use the python binding in a while but you probably need to set 
> the PYTHONPATH to where pydealii is located. You may also need to add the 
> path to libdealii.so to LD_LIBRARY_PATH
>
> Where is this library located, and how should I go about setting up a 
>> (conda) environment to be able to use PyDealII? Are there dependencies?
>>
> The library should be located in the same directory where you installed 
> deal.II. Look somewhere in /path/to/install/lib . I don't use conda so I 
> can't help you with that. What kind of dependencies are you talking about? 
> PyDealII only depends on deal.II and Boost.Python, nothing else
>
> Best,
>
> Bruno 
>

-- 
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/c39ec853-1838-4d02-84f6-15589fc1ecd8n%40googlegroups.com.

[deal.II] Python bindings

[Corbin Foucart]

Hello everyone,

I'd like to write some custom Python bindings for deal.ii, and as a first 
step, I've been trying to reproduce the ipython notebook for step-49. 

I've built deal.ii version 9.3.0 with the 
-DDEAL_II_COMPONENT_PYTHON_BINDING=ON configuration, yet when I run the 
notebook, the Jupyter kernel is unable to locate the PyDealII module, that 
is, the notebook cell

> import PyDealII.Release as dealii

fails with an import error. Where is this library located, and how should I 
go about setting up a (conda) environment to be able to use PyDealII? Are 
there dependencies?

Any guidance would be appreciated!
Corbin

-- 
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/f2148581-7815-404c-b1f9-783026235308n%40googlegroups.com.

Re: [deal.II] Solving separate FEM problems on subdomains

[Corbin Foucart]

Bruno and Wolfgang,

Thank you for the helpful suggestions--I had not thought of the 
block-diagonal linear system approach in particular, but I think it is a 
clever way to do what I had previously planned. I don't anticipate it being 
too expensive, even unordered, but if so, I will implement one of the other 
two approaches.

Corbin

On Tuesday, March 22, 2022 at 10:38:32 PM UTC-4 Wolfgang Bangerth wrote:

> On 3/21/22 16:47, Corbin Foucart wrote:
> > 
> > I have an extruded mesh, and I would like to solve a finite element 
> problem 
> > separately on each 'column' as part of a larger finite element solve in 
> which 
> > all DoF are coupled through a discontinuous Galerkin scheme for an 
> elliptic 
> > equation--so multiply defined values on the boundaries of each column 
> are not 
> > a problem.
> > 
> > * I already have the columns tagged by material ID
> > * I imagine I should instantiate a list of new Triangulations, 
> DoFHandlers,
> > and solver classes for each column and then map the solutions back to the
> > original domain
> > * *Is there a good way to extract a separate triangulation for each 
> column*?
> > I was looking for something similar to extract_boundary_mesh but with a
> > material_id
> > 
> > Alternatively, it seems that I could assign the DoFs in the globally 
> coupled 
> > DoFHandler with DoFTools::get_subdomain_association(). *Would this be a 
> better 
> > way to set up each finite element problem by column?* It's unclear to me 
> if I 
> > could use that information to set up an independent linear system for 
> each 
> > subdomain.
>
> Corbyn,
> you could go the way of creating triangulations, DoFHandlers, etc, for 
> each of 
> the columns. But you might also want to consider one of the following two 
> approaches:
>
> * Since you're already doing a DG approach, consider building your matrix 
> with 
> only the interface terms that correspond to the "bottom" and "top" of each 
> of 
> your cells. This means that you don't build any term that corresponds to 
> coupling columns to each other. Solving such a linear system then 
> corresponds 
> to solving all of the 1d problems at the same time. (I will note that such 
> a 
> matrix would also make for a good preconditioner; it corresponds to one of 
> the 
> directions one would solve for in the ADI method.)
>
> Conceptually, the matrix you obtain using the approach above is block 
> diagonal 
> where each block corresponds to one column. In practice, unless you do 
> something special, the matrix will not actually look block diagonal 
> because 
> the order of DoFs does not respect the column-by-column order. But it 
> wouldn't 
> be very difficult to implement a function that renumbers DoFs in such a 
> way 
> that they are numbered column by column, and then you would indeed have a 
> block diagonal matrix.
>
> * If you don't like this approach, you can create a mapping from the DoFs 
> of 
> each column to a column-local numbering from zero to n_dofs_per_column. 
> Instead of writing into global vectors and matrices in the method above 
> where 
> you omit the lateral face terms, you can translate writing into smaller 
> matrices and vectors that correspond to just one column. You can then 
> solve 
> the resulting linear system. The mapping of DoF indices then also allows 
> you 
> to get back at what the global index of these DoFs was.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/3719937b-8b27-4f3c-a6be-705dfa3157e7n%40googlegroups.com.

[deal.II] Solving separate FEM problems on subdomains

[Corbin Foucart]

Hello all,

I have an extruded mesh, and I would like to solve a finite element problem 
separately on each 'column' as part of a larger finite element solve in 
which all DoF are coupled through a discontinuous Galerkin scheme for an 
elliptic equation--so multiply defined values on the boundaries of each 
column are not a problem.

   - I already have the columns tagged by material ID
   - I imagine I should instantiate a list of new Triangulations, 
   DoFHandlers, and solver classes for each column and then map the solutions 
   back to the original domain
   - *Is there a good way to extract a separate triangulation for each 
   column*? I was looking for something similar to extract_boundary_mesh 
   but with a material_id
   
Alternatively, it seems that I could assign the DoFs in the globally 
coupled DoFHandler with DoFTools::get_subdomain_association(). *Would this 
be a better way to set up each finite element problem by column?* It's 
unclear to me if I could use that information to set up an independent 
linear system for each subdomain.

Any guidance would be appreciated!
Corbin

-- 
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/6f82d2ba-4e43-4169-8511-a504031e4ce6n%40googlegroups.com.

Re: [deal.II] system_to_component_index() usage in face_worker for DG method

[Corbin Foucart]

Sylvain,

Oh, I think I know what those values appears and it is my bad ...
> It is probably not 22116 but 22 and 116 and it is probably not 89130 but 
> rather 89 and 130. I used something like *std::cout << dof_index[0] << 
> std::endl* to get those values and sometimes they just get attached to 
> each other like if the *std::endl* was skipped (not sure why, maybe 
> because I have several threads ?). I did not realise it before reading your 
> message but it actually makes sense because I was not getting an error for 
> those values and my Assert was not triggered. I should have seen it, sorry.
>

 
I ran into this exact issue a few weeks ago. For printing statements while 
debugging in assembly loops, check the signatures of functions like 
MeshWorker::run and MeshWorker::mesh_loop. You should be able to pass 
parameters that force non-multithreaded execution, and your print 
statements will make sense. I found this helpful for debugging.

Best,
Corbin

-- 
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/7ec902ce-e1fa-4baa-a669-4b2b288479dfn%40googlegroups.com.

Re: [deal.II] Assembling explicit advection term independently

[Corbin Foucart]

> That would be fantastic! There are some good resources linked to from 
> here, if 
> you're unsure about the process: 
> https://dealii.org/participate.html


I'm ready to start working on this. I've read the participation FAQ, forked 
the repo, etc. This might be a basic question, but how should I go about 
re-compiling the library code once I've made changes to the source? Is 
there a makefile that makes and runs all the unit tests? I presume the new 
code that I use to test the functions should be added in 
dealii/tests/feinterface?

Thank you for the guidance in advance!
Corbin

 

-- 
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/75248a20-1bc9-4cdf-b540-359787b1dd21n%40googlegroups.com.

Re: [deal.II] Looping over cells and faces: MeshWorker vs WorkStream

[Corbin Foucart]

I'm writing a similar code to Praveen, and had some additional questions, 
having looked at the documentation for mesh_loop as well as step-47, 
step-51, and the code gallery entry for the second-order LDG solver.

For a scheme that has a mixed-formulation element-local system (Q,U) as 
does the LDG scheme in the code gallery, it is desirable to instead 
eliminate the gradient Q on each cell and assemble a linear system only in 
U. This amounts to a Schur complement taken on each cell--that is, 
[[ A   B ]]  [Q]  = [ 0 ]
[[ -B^T  C ]]  [U] [ F ]

and a cell contribution to the LHS of the global linear system KU = L, 
where K = C + B^T A^{-1} B. My question is the following:

   - It seems that since cell interior data (A inverse) is needed to 
   premultiply the B operator (which contains face terms), it makes more sense 
   to loop over each cell and handle the faces of the cell manually, like in 
   step-51, which performs a similar element-local static condensation.
   - In the spirit of step-51, my idea is to use mesh_loop, but only use 
   the cell-assembly part and pass dummy functions to the face_worker and 
   boundary_worker arguments. Is that a sensible way to do it, or should I be 
   using a different MeshWorker function?
   - Or alternatively, is there a way to use the face_worker function, but 
   preserve the data on the cell to the left and right? It seems expensive to 
   recompute the inverse mass matrix on every interface, for each adjacent 
   cell.

Best,
Corbin

On Tuesday, June 1, 2021 at 12:06:42 AM UTC-4 Wolfgang Bangerth wrote:

> On 5/31/21 9:36 PM, Praveen C wrote:
> > 
> > It this how it should be used ?
> > 
>
> Something like this. You might want to look at step-47 which assembles 
> face 
> terms for interior faces.
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/914360a9-7c11-4779-bdab-e3f4d137c58bn%40googlegroups.com.

Re: [deal.II] Re: Extracting a component of FESystem for use elsewhere

[Corbin Foucart]

Hi Paras,

Thanks for the reply!

```
> dealii::FEValuesExtractorScalar pressureExtractor(2); //assuming you have 
> three scalar fields and pressure is the last one
> dealii::Vector solVecPerCell(nDOFsPerCell);
> feValues[pressureExtractor].get_function_values(totalSolution, 
> solVecPerCell); // here only the values corresponding to pressure are 
> filled.
> ```
>
 
Suppose I have 2 DG elements of order 3, i.e., 16 DoFs per element. U, v, p 
-> 2*(16)*3 = 96 DoF in totalSolution.

It's my understanding that this will extract the pressure values from 
totalSolution into the *cell-local* Vector solVecPerCell of length 
16 (assuming FEValues was instantiated with a quadrature rule that 
contained the nodal locations, rather than the usual quadrature points). 
However, this isn't what I want to do. I want to take all 32 of the p nodal 
values and add them to another field with a DoFHandler instantiated using a 
single FE_DGQ element, rather than a system--therefore with a solution 
Vector containing 32 DoF.

Maybe the way to go is a BlockVector object as you alternatively suggest, 
but my problem doesn't inherently have blocks. I'm simply looking to add, 
for example, the pressure field to another scalar DG field unrelated to the 
system; I'd rather keep the FESystem organized by (u,v,p) on each cell.

Best,
Corbin  

-- 
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/49866be9-6017-4dde-8f10-8b421863e6a4n%40googlegroups.com.

[deal.II] Re: Extracting a component of FESystem for use elsewhere

[Corbin Foucart]

I've found one way to do this, but I imagine it's not optimal:

   - call DoFRenumbering::component_wise(fesys_dofh)
   - Do the finite element solve
   - then copy from one vector to another using the number of dofs in that 
   component and an offset
   - This isn't great, however, because I'd rather store the fesys dofh 
   cell-by-cell for memory-friendliness, as almost all other operations are 
   elemental

If there's a better way, I'd be glad to know :)

Best,
Corbin 

On Wednesday, June 16, 2021 at 11:21:56 PM UTC-4 Corbin Foucart wrote:

> Hi everyone,
>
> This might be a simple question, but I haven't found a straightforward 
> answer in the "handling vector-valued problems" documentation.
>
> If I have a Vector which stores the solution to an FESystem, say 3 
> FE_DGQ fields for u, v, p, and I want to extract a new Vector which 
> corresponds to only the pressure, what's the best way to do this?
>
> I'm not looking to write the data out for visualization, how to use the 
> DataComponentInterpretation classes is clear to me. What I want to do is 
> extract a separate Vector which I can then add to another DG scalar 
> field of the same type as p in the (u,v,p) system.
>
> I don't need p at the quadrature points, only the nodal values to add to 
> another nodal field.
>
> Any advice would be appreciated!
> Corbin
>

-- 
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/2bb1442a-830e-46a2-8639-8eec8b101e15n%40googlegroups.com.

[deal.II] Extracting a component of FESystem for use elsewhere

[Corbin Foucart]

Hi everyone,

This might be a simple question, but I haven't found a straightforward 
answer in the "handling vector-valued problems" documentation.

If I have a Vector which stores the solution to an FESystem, say 3 
FE_DGQ fields for u, v, p, and I want to extract a new Vector which 
corresponds to only the pressure, what's the best way to do this?

I'm not looking to write the data out for visualization, how to use the 
DataComponentInterpretation classes is clear to me. What I want to do is 
extract a separate Vector which I can then add to another DG scalar 
field of the same type as p in the (u,v,p) system.

I don't need p at the quadrature points, only the nodal values to add to 
another nodal field.

Any advice would be appreciated!
Corbin

-- 
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/6d96604a-cd77-4532-aa16-a6fb638a6ff9n%40googlegroups.com.

[deal.II] Deal.II site appears to be down

[Corbin Foucart]

Hello everyone,

It seems dealii.org is experiencing some problems or is down (I'm on the 
east coast of the US)--checking 
on https://downforeveryoneorjustme.com/dealii.org, it seems the server is 
not responding and the issue isn't local to me.

I'm posting this for the site admins just in case I'm the first to notice.

Best,
Corbin 

-- 
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/1c320f3d-c5ba-4dd9-8cc1-511aee07c0ddn%40googlegroups.com.

Re: [deal.II] Assembling explicit advection term independently

[Corbin Foucart]

> I forgot that earlier. Could you come up with a small program that shows 
> the 
> segfault? We should have caught the wrong vector size via an assertion, 
> and 
> I'd like to add that assertion. 
>

Certainly--please find attached. 

To clarify, deal.II *did* throw some type of exception, I just wasn't able 
to interpret it and opted to use GDB instead; I'm a Python programmer and 
not a C++ expert by any stretch of the imagination :)

The program terminates with:

An error occurred in line <3196> of file 
 in function
void dealii::internal::do_function_values(const typename 
VectorType::value_type*, const dealii::Table<2, double>&, const 
dealii::FiniteElement&, const std::vector&, 
dealii::ArrayView, bool, unsigned int) [with int dim = 2; int 
spacedim = 2; VectorType = dealii::Vector; typename 
VectorType::value_type = double]
The violated condition was: 
static_cast::type)>::type>(values.size()) == 
static_cast::type)>::type>(n_quadrature_points)
Additional information: 
Dimension 0 not equal to 3.

Stacktrace:
---
#0  
/home/foucartc/CorbinFoucart_repos/dealii_applications/lib/libdeal_II.g.so.9.2.0:
 
void dealii::internal::do_function_values<2, 2, dealii::Vector 
>(dealii::Vector::value_type const*, dealii::Table<2, double> 
const&, dealii::FiniteElement<2, 2> const&, std::vector > const&, 
dealii::ArrayView, dealii::MemorySpace::Host>, bool, 
unsigned int)
#1  
/home/foucartc/CorbinFoucart_repos/dealii_applications/lib/libdeal_II.g.so.9.2.0:
 
void dealii::FEValuesBase<2, 2>::get_function_values 
>(dealii::Vector const&, 
std::vector::value_type>, 
std::allocator::value_type> > >&) 
const
#2  ./step: DGVelocityField<2>::access_interface_values()
#3  ./step: main

Aborted (core dumped)


As for implementing FEInterfaceValues::get_function_jumps et al: I have a 
deadline by which I'd really like these set of codes to be working, so I 
can't work on this in the next few weeks. However, I will have to implement 
that functionality anyway, so I wouldn't mind working on this sometime in 
July.

Best,
Corbin

-- 
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/1bfe5faa-2f0e-4b8a-a0b6-015c73aaa469n%40googlegroups.com.
// unification to a depth / surface mesh class
//
#include 
#include 


// deal ii imports
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
struct DGData
{
  FEValues fe_values;
  FEFaceValues fe_face_values;
  FEInterfaceValues fe_interface_values;

  // stores mapping between faces and the cell interior basis 
  // functions which have support on each face
  std::vector> fe_local_support_on_face;

  DGData(
const FiniteElement &fe,
const Quadrature&cell_quadrature,
const Quadrature  &face_quadrature,
const UpdateFlags  cell_update_flags = update_values 
| update_gradients 
| update_quadrature_points 
| update_JxW_values,
const UpdateFlags interface_update_flags = update_values 
| update_quadrature_points
| update_JxW_values 
| update_gradients
| update_normal_vectors)
: fe_values(fe, cell_quadrature, cell_update_flags)
, fe_face_values(fe, face_quadrature, interface_update_flags)
, fe_interface_values(fe, face_quadrature, interface_update_flags)
, fe_local_support_on_face(GeometryInfo::faces_per_cell)
{ }

};

template 
class DGVelocityField
{
  public:
DGVelocityField(
  const unsigned int degree,
  Triangulation &grid) 
: triangulation(grid)
, fe_sys(FE_DGQ(degree), dim)
, dofh(triangulation)
, quadrature(degree + 1)
, face_quadrature(degree + 1)
, dgd(fe_sys, quadrature, face_quadrature)
{ 
  // system DOF
  dofh.distribute_dofs(fe_sys);
  velocity_field.reinit(dofh.n_dofs());
}

Triangulation &triangulation;
FESystem fe_sys;
DoFHandler dofh;
Vector velocity_field;
const QGauss quadrature;
const QGauss face_quadrature;
DGData dgd;

void access_interface_values();


};

template 
void DGVelocityField::access_interface_values()
{
  for (const auto &cell : dofh.active_cell_iterators())
  {
dgd.fe_values.reinit(cell);
for (const auto face_index : GeometryInfo::face_indices())
{
  if (!cell->face(face_index)->at_boundary())
  {
dgd.fe_interface_values.reinit(cell,
  face_index,
   

Re: [deal.II] Assembling explicit advection term independently

[Corbin Foucart]

Upon reading the documentation closer, I think I've found the way to access 
u^k (+/-) across each interface


>- I can't seem to be able to access u^k(+) and u^k(-) using a 
>FEInterfaceValues fiv, because it's vector-valued
>- this same procedure works fine for scalars. On any face, I can use 
>something like fiv.get_face_values(1).get_function_values(scalar_field, 
>qpt_vector) to get scalar(-) at the face quadrature points (and use 
>get_face_values(0) for the scalar(+) value)
>
> Am I missing something? There must be an easy way to access jumps and 
> averages of vector-valued DG data on the right-hand side. Step-47 almost 
> does this, but all the vector-valued jumps appear in the bilinear form 
> rather than the right-hand side.
>

Indeed,
std::vector> velocity_qpt_vals(n_face_qpts, 
Vector(dim));
fiv.get_fe_face_values(1).get_function_values(velocity_field, 
velocity_qpt_vals); 

accesses the values; the issue was I hadn't explicitly instantiated the 
Vectors with length dim, which was causing a segfault.  However, 
I'm curious -- is there a better way to compute the jumps and averages of 
DG data other than requesting them via fe_face_values and computing them 
"by hand"?

Thank you,
Corbin

-- 
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/789d59fa-2bc9-4dca-a577-504713b29e87n%40googlegroups.com.

Re: [deal.II] Assembling explicit advection term independently

[Corbin Foucart]

Thanks for the detailed and helpful answer, Prof. Bangerth!

Yes, I thought about it in the past few days and arrived at the same 
conclusion--that it does not make sense to represent the explicit advective 
term as a single finite element field. I settled on a separate RHSAssembler 
abstract base class with virtual methods to be overloaded for the volume 
and face integrals on each cell. The main solvers simply request these 
integrals from this class in the course of WorkStream::run.

The subtlety with respect to these advection body/face terms arises in the 
following: as I'm sure you're well-aware, HDG methods involve an FESystem 
(Q, U) on each element which is parametrized to a global system in the 
skeleton space flux U_hat only (I used an FE_FaceQ to represent it). So all 
of the right-hand side forcing is *not* communicated to the global linear 
system directly, but indirectly through the U equation of the (Q, U) local 
system (which is inverted to eliminate everything in terms of U_hat). So I 
can't assemble the arbitrary RHS contributions directly into the linear 
system like in the tutorials. I have to compute them and add them to the 
RHS of the element-local system, like in step-51.

To that end, I'm having some trouble. I'd like to be able to compute 
numerical fluxes of the known discontinuous field u^k (i.e., sum jumps and 
averages of u^k) at the quadrature points on each face. However,

   - I can't seem to be able to access u^k(+) and u^k(-) using a 
   FEInterfaceValues fiv, because it's vector-valued
   - this same procedure works fine for scalars. On any face, I can use 
   something like fiv.get_face_values(1).get_function_values(scalar_field, 
   qpt_vector) to get scalar(-) at the face quadrature points (and use 
   get_face_values(0) for the scalar(+) value)

Am I missing something? There must be an easy way to access jumps and 
averages of vector-valued DG data on the right-hand side. Step-47 almost 
does this, but all the vector-valued jumps appear in the bilinear form 
rather than the right-hand side.

If not, I could maybe extract the jumps and averages component by 
component, using a storage scheme like step-35, but I feel like there must 
be a better way.

Many thanks,
Corbin

On Wednesday, June 9, 2021 at 9:55:41 PM UTC-4 Wolfgang Bangerth wrote:

> On 6/8/21 3:33 AM, Corbin Foucart wrote: 
> > 
> > I've written generic, implicit HDG/DG solvers that solve a second-order 
> > diffusion type equation $\frac{du}{dt} - \nabla^2 u = f$ (where f is a 
> generic 
> > DG finite element field). As I'm expanding the code to a set of 
> projection 
> > methods, I'd like this $f$ to represent an explicitly computed advection 
> term 
> > on the right-hand side $\nabla\cdot (u^k \otimes u^k)$ at time $k+1$. 
> > 
> > I'd like to assemble the weak form of the advection term, something like 
> > $-(∇v, u^k\otimes u^k)_K +(v, F^*(u)\cdot n)_{\partial K}$ (where F^*(u) 
> is 
> > the convective numerical flux) on the right hand side for every element 
> K and 
> > store that as a DG field. 
>
> This is conceptually not the right approach. A finite element field is a 
> function of the form 
> u_h(x) = \sum_j U_j \psi_j(x) 
> where U_j is a vector of coefficients. But you don't have a field of this 
> form: You have a vector 
> F_j 
> = sum_K -(∇psi_j, u^k\otimes u^k)_K +(\psi_j, F^*(u)\cdot n)_{\partial K} 
> but this is a vector of (weighted) integrals of the previous solution, not 
> a 
> vector of expansion coefficients. 
>
> In other words, your vector F_j is a vector alright, but it does *not* 
> corresponding to a finite element field, whether continuous or 
> discontinuous. 
>
>
> > To avoid breaking the abstraction of the implicit 
> > solver classes, is it possible to compute this term outside the main 
> solver 
> > classes and pass it in as a DG Field represented as a Vector? 
>
> What you plan to do is no different, really, than any of the other right 
> hand 
> side vectors we assemble. For example, in step-4, we assemble a right hand 
> side vector ("system_rhs") in assemble_system(). The only difference is 
> that 
> in your case, this right hand side vector depends on the previous 
> solution, 
> but conceptually you're doing the same thing. 
>
> Of course, it may be the case that you're not doing any other assembly 
> operations in each time step, for example because you keep using the same 
> matrix every time. But you can always write a function of the form 
> assemble_rhs() that only assembles a rhs vector. As just one example, the 
> new 
> step-77 program has a function compute_residual() that does essentially 
> this. 
>
> Best 
> Wolfgang 
>
> -- 
> --

[deal.II] Assembling explicit advection term independently

[Corbin Foucart]

Hello all, 

I've written generic, implicit HDG/DG solvers that solve a second-order 
diffusion type equation $\frac{du}{dt} - \nabla^2 u = f$ (where f is a 
generic DG finite element field). As I'm expanding the code to a set of 
projection methods, I'd like this $f$ to represent an explicitly computed 
advection term on the right-hand side $\nabla\cdot (u^k \otimes u^k)$ at 
time $k+1$. 

I'd like to assemble the weak form of the advection term, something like 
$-(∇v, u^k\otimes u^k)_K +(v, F^*(u)\cdot n)_{\partial K}$ (where F^*(u) is 
the convective numerical flux) on the right hand side for every element K 
and store that as a DG field. To avoid breaking the abstraction of the 
implicit solver classes, is it possible to compute this term outside the 
main solver classes and pass it in as a DG Field represented as a 
Vector?

I suspect this can be done, but I'm unsure how to populate a Vector 
directly with terms like these, since there's no global assembly to be 
done. Perhaps it can be filled directly with MeshWorker::mesh_loop?

Any advice would be appreciated!
Corbin


-- 
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/c2501b58-79fc-4840-ad96-4ccabff73c11n%40googlegroups.com.

Re: [deal.II] Extracting a component from Block FESystem and DofHandler

[Corbin Foucart]

Hi all,

I'm doing something very similar to Spencer. The answers provided by Prof. 
Bangerth and Bruno make sense, up to one question:

The first two yes. If you do it right, i.e., ensure that the velocity 
> DoFHandler numbers degrees of freedom in the same order as the system 
> DoFHandler does, then you won't need a separate solution vector -- you can 
> just use a block of your block vector.
>

 How would one go about ensuring that the velocity DoFHandler numbers the 
DOF in the same order as the system DoFHandler does? Is there a way to 
build the former from the latter?

Best,
Corbin

-- 
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/20e3e511-4368-4e2f-84e5-e9e96b813b6an%40googlegroups.com.

Re: [deal.II] Accessing cell data from its ID

[Corbin Foucart]

To those who might be following in my footsteps and encounter similar 
questions:

   - For understanding how accessors/iterators work, the following page was 
   indispensible: 
   https://dealii.org/developer/doxygen/deal.II/group__Iterators.html -- I 
   don't know how I missed it the first time around
   - For transferring bulk fields to surface fields and vice-versa, the 
   GridGenerator::extract_boundary_mesh does *not* keep the dof orientation 
   between the surface cells and the bulk cell faces from which they were 
   extracted
   - However, the map that it returns can be coupled with a DoFHandler loop 
   to construct a DoF map between the extracted surface field and bulk fields, 
   which makes data transfer between the two painless once it's done
   - If all else fails, you can do a one-time comparison of support points 
  to build the DOF map, see the documentation for 
  DoFTools::map_to_support_points
  
I found that this approach was easier than the periodic boundary condition 
approach, although both worked.

Thanks, Simon and David!
Corbin

On Saturday, May 1, 2021 at 9:36:02 PM UTC-4 Corbin Foucart wrote:

> Thank you, Simon and David. Both answers were helpful. I will try 
> implementing it both ways for practice.
>
> Best,
> Corbin
>
> On Thursday, April 29, 2021 at 8:40:48 AM UTC-4 d.arnd...@gmail.com wrote:
>
>> Corbin,
>>
>> if you only ever need the DoF mapping periodic 
>> DoFTools::make_periodicity_constrainst (
>> https://www.dealii.org/current/doxygen/deal.II/namespaceDoFTools.html#a03324e6e060c6a55191b2c60ad871eab)
>>  
>> sound like the right tool for you to use.
>> Just create a separate AffineConstraints object to be filled and hand it 
>> to that function. Then, you can loop through the filled object to either 
>> use it as a map directly or you can create a separate data structure. Note 
>> that the constraints might not be in the order (bottom -> top) you need 
>> them to be. So you might need to adjust that yourself.
>>
>> Best,
>> Daniel
>>
>>
>> Am Do., 29. Apr. 2021 um 04:48 Uhr schrieb Corbin Foucart <
>> corbin@gmail.com>:
>>
>>> Hi all,
>>>
>>> I'm looking to map DoF from the bottom of an extruded volumetric mesh to 
>>> the corresponding DOF on the surface of the mesh. Right now, I'm saving the 
>>> cell "number"/ID of each cell in the bottom layer of the volume mesh, and 
>>> want to match it with a surface cell by comparing their centers
>>>
>>>- I could loop over all the volume cells until I find the one 
>>>matching the surface cell center, but as I have the cell numbers, I 
>>> would 
>>>rather access the cell center from the cell id (which looks to me like 
>>> an 
>>>unsigned int)
>>>- something like get_cell(cell_id)->center  --- is this possible?
>>>- Alternatively, I was thinking that a periodic boundary condition 
>>>linking the bottom of the volume mesh with the top might be able to 
>>> achieve 
>>>the same thing, but I wasn't sure how to go about setting it up
>>>
>>> Any pointers would be appreciated!
>>> Corbin
>>>
>>> -- 
>>> 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/122e1b62-e672-446b-8c25-ce6012c1fe3fn%40googlegroups.com
>>>  
>>> <https://groups.google.com/d/msgid/dealii/122e1b62-e672-446b-8c25-ce6012c1fe3fn%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/c56ca4f3-4436-435b-be5b-9c16dc8c5287n%40googlegroups.com.

Re: [deal.II] Accessing cell data from its ID

[Corbin Foucart]

Thank you, Simon and David. Both answers were helpful. I will try 
implementing it both ways for practice.

Best,
Corbin

On Thursday, April 29, 2021 at 8:40:48 AM UTC-4 d.arnd...@gmail.com wrote:

> Corbin,
>
> if you only ever need the DoF mapping periodic 
> DoFTools::make_periodicity_constrainst (
> https://www.dealii.org/current/doxygen/deal.II/namespaceDoFTools.html#a03324e6e060c6a55191b2c60ad871eab)
>  
> sound like the right tool for you to use.
> Just create a separate AffineConstraints object to be filled and hand it 
> to that function. Then, you can loop through the filled object to either 
> use it as a map directly or you can create a separate data structure. Note 
> that the constraints might not be in the order (bottom -> top) you need 
> them to be. So you might need to adjust that yourself.
>
> Best,
> Daniel
>
>
> Am Do., 29. Apr. 2021 um 04:48 Uhr schrieb Corbin Foucart <
> corbin@gmail.com>:
>
>> Hi all,
>>
>> I'm looking to map DoF from the bottom of an extruded volumetric mesh to 
>> the corresponding DOF on the surface of the mesh. Right now, I'm saving the 
>> cell "number"/ID of each cell in the bottom layer of the volume mesh, and 
>> want to match it with a surface cell by comparing their centers
>>
>>- I could loop over all the volume cells until I find the one 
>>matching the surface cell center, but as I have the cell numbers, I would 
>>rather access the cell center from the cell id (which looks to me like an 
>>unsigned int)
>>- something like get_cell(cell_id)->center  --- is this possible?
>>- Alternatively, I was thinking that a periodic boundary condition 
>>linking the bottom of the volume mesh with the top might be able to 
>> achieve 
>>the same thing, but I wasn't sure how to go about setting it up
>>
>> Any pointers would be appreciated!
>> Corbin
>>
>> -- 
>> 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/122e1b62-e672-446b-8c25-ce6012c1fe3fn%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/122e1b62-e672-446b-8c25-ce6012c1fe3fn%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/602b1076-c4ec-4451-bd86-9c391b7644d6n%40googlegroups.com.

[deal.II] Re: Segfault upon printing Point in GDB

[Corbin Foucart]

And the file is here.

On Saturday, May 1, 2021 at 12:07:13 AM UTC-4 Corbin Foucart wrote:

> Hello all,
>
> When I attempt to use the GDB 'print' command to print information about a 
> Point object, I'm having some trouble. For example, suppose I've defined 
>
> const Point p_ex(1, 1);
>
>- If I run GDB and attempt 'print p_ex[0]' the output is 'Could not 
>find operator[]'
>- If I attempt 'print p_ex', GDB segfaults
>
> I imagine printing these sorts of objects is a common procedure, so I'm 
> hoping someone might be able to help me shed light on this!
>
> More details:
> I've attached a minimal code example which reproduces the issue. I compile 
> the program without issue on Ubuntu 18.04.5 LTS and the exact GDB commands 
> I use are:
> > gdb ./step
> >> break step.cc:26
> >> r
> Breakpoint 1, main () at src/step.cc:27
> 27  return 0;
> >> print bottom_left[0]
> Could not find operator[].
> >> print bottom_left
> $1 = {> = {
> Segmentation fault (core dumped)
>
> Best,
> Corbin
>

-- 
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/04b63e17-e26a-4ace-844c-c5666d318c6en%40googlegroups.com.
// unification to a depth / surface mesh class
//
#include 
#include 
#include 

// deal ii imports
#include 

// local imports
//#include "ocean_mesh.h"

using namespace dealii;

int main()
{
  const int dim = 2;

  // create an [NN]^2 grid
  unsigned int n_x_subdivisions = 2;
  unsigned int n_y_subdivisions = 2;
  std::vector repetitions{ n_x_subdivisions, n_y_subdivisions };
  
  const Point bottom_left(-1.0, -1.0);
  const Point   top_right( 0.0,  0.0);

  return 0;
}

[deal.II] Segfault upon printing Point in GDB

[Corbin Foucart]

Hello all,

When I attempt to use the GDB 'print' command to print information about a 
Point object, I'm having some trouble. For example, suppose I've defined 

const Point p_ex(1, 1);

   - If I run GDB and attempt 'print p_ex[0]' the output is 'Could not find 
   operator[]'
   - If I attempt 'print p_ex', GDB segfaults

I imagine printing these sorts of objects is a common procedure, so I'm 
hoping someone might be able to help me shed light on this!

More details:
I've attached a minimal code example which reproduces the issue. I compile 
the program without issue on Ubuntu 18.04.5 LTS and the exact GDB commands 
I use are:
> gdb ./step
>> break step.cc:26
>> r
Breakpoint 1, main () at src/step.cc:27
27  return 0;
>> print bottom_left[0]
Could not find operator[].
>> print bottom_left
$1 = {> = {
Segmentation fault (core dumped)

Best,
Corbin

-- 
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/31456099-c733-4523-b9d0-ee0ac7947eebn%40googlegroups.com.

[deal.II] Accessing cell data from its ID

[Corbin Foucart]

Hi all,

I'm looking to map DoF from the bottom of an extruded volumetric mesh to 
the corresponding DOF on the surface of the mesh. Right now, I'm saving the 
cell "number"/ID of each cell in the bottom layer of the volume mesh, and 
want to match it with a surface cell by comparing their centers

   - I could loop over all the volume cells until I find the one matching 
   the surface cell center, but as I have the cell numbers, I would rather 
   access the cell center from the cell id (which looks to me like an unsigned 
   int)
   - something like get_cell(cell_id)->center  --- is this possible?
   - Alternatively, I was thinking that a periodic boundary condition 
   linking the bottom of the volume mesh with the top might be able to achieve 
   the same thing, but I wasn't sure how to go about setting it up

Any pointers would be appreciated!
Corbin

-- 
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/122e1b62-e672-446b-8c25-ce6012c1fe3fn%40googlegroups.com.

Re: [deal.II] Bug in FEInterfaceValues?

[Corbin Foucart]

I took a look at the unit tests for FEInterfaceValues in
tests/feinterface/fe_interface_values_0*.cc, and found that the following
modification to the code in the face loop, changing how reinit and the
boundary check is done:

for ( cell loop)
  for ( face loop)
  ...

  // detect top face
  if (nrml[0][1] > 0 && !(cell->at_boundary(face_index)))
  {
std::cout << std::endl;
sd.fe_interface_values.reinit(cell,
  face_index,
  numbers::invalid_unsigned_int,
  cell->neighbor(face_index),

cell->neighbor_of_neighbor(face_index),
  numbers::invalid_unsigned_int);
   ...

sd.fe_interface_values.get_fe_face_values(1).get_function_values(field,
uhat_qpt);
   }

Works fine at getting the values from the other side of the interface.
Perhaps I used the reinit(...) function incorrectly.

Best,
Corbin


On Tue, Apr 13, 2021 at 9:56 PM Corbin Foucart 
wrote:

> Hello all,
>
> I have found an unexpected behavior in using the FEInterfaceValues object.
> Namely, it seems that a line like,
>
> fe_interface_values.get_fe_face_values(1).get_function_values(...);
>
> returns the following error
>
> > The violated condition was:
> > cell_index == 0 || !at_boundary()
> > Additional information:
> > You are on a boundary, so you can only ask for the first
> FEFaceValues object.
>
> ...if the cell is on a boundary, but the face is an interior face.
>
> I've attached code for a small test case illustrating the issue. It
> consists of a 2x2 2D mesh, and attempts to use FEInterfaceValues to get
> field values from the cell "above".
>
> Is this indeed a bug, or am I misunderstanding the intended usage of
> FEInterfaceValues? I'm using deal.II, version 9.2.0.
>
> Best,
> Corbin
>
> --
> 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/vGc16ulvglE/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/22813e93-5959-4244-aa45-1652d3e8c7b8n%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/22813e93-5959-4244-aa45-1652d3e8c7b8n%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/CAE%2BMPofcZ5nizy0XyfDtiXpn5S%2BXiyhNoC3CW-EgO1YC4xaDqA%40mail.gmail.com.

Re: [deal.II] Solving a local DG problem on each cell

[Corbin Foucart]

Prof. Bangerth,

Thank you for the advice. I've spent some time with step-61 (and step-51 as 
well) and I think I understand the distinction between the two classes. I 
am using FEFaceValues to assemble the edge term from the interior faces, 
and attempting to use FEInterfaceValues to access the values from the cell 
face above.

One thing I am still unclear on is if FEInterfaceValues is the best (or 
only) way to access the face values of a field on a neighboring cell, or if 
this can be done manually. Right now, I'm using an approach like 

std::vector 
field_at_qpts(sd.fe_interface_values.n_quadrature_points);
sd.fe_interface_values.get_fe_face_values(1).get_function_values(field, 
field_at_qpts);

where 'sd' is a typical scratchdata object with an FEInterfaceValues 
instance as a class attribute. Is this typically how neighboring data is 
accessed for DG-like methods?

Thank you,
Corbin

On Friday, April 9, 2021 at 7:17:54 PM UTC-4 Wolfgang Bangerth wrote:

>
> Corbin,
>
> > I'm attempting to solve a local discontinuous Galerkin problem on each 
> cell in 
> > my mesh. I've attached a screenshot with the cell-local equation.
> > 
> > Screen Shot 2021-04-09 at 13.56.15.png
> > 
> > I am solving these cell problems in a very specific order to explicitly 
> > propagate information downwards in depth, so assembling an entire system 
> like 
> > in step-12 is not what I am looking to do.
> > 
> > I'm able to assemble the interior terms without a problem, but I'm a bit 
> > unclear on how to assemble the edge terms.
> > 
> > Should I be using FEFaceValues or FEInterfaceValues? On the top of the 
> cell, I 
> > need the face data from the cell above (for u_hat), whereas on the 
> bottom 
> > cell, I need to integrate over the bottom face. In either case, it's 
> unclear 
> > to me how to map the integrals involving face basis functions to the 
> volume 
> > basis functions in the cell.
>
> The distinction is that FEFaceValues evaluates the values/derivatives of 
> solutions and shape functions on faces *as seen from one cell*, whereas 
> FEInterfaceValues does so for averages/differences between adjacent cells. 
> In 
> your case, it seems to me that you are only interested in evaluating \hat 
> u as 
> seen from one side of the face, so FEFaceValues is the right tool.
>
> You might be interested in spending an hour or two on step-61, which also 
> solves local problems (in a different context, but maybe there is a thing 
> of 
> two for you to learn from).
>
> Best
> W.
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/a25dffdb-74da-438a-b576-521858337b60n%40googlegroups.com.

[deal.II] Bug in FEInterfaceValues?

[Corbin Foucart]

Hello all,

I have found an unexpected behavior in using the FEInterfaceValues object. 
Namely, it seems that a line like, 

fe_interface_values.get_fe_face_values(1).get_function_values(...);

returns the following error 

> The violated condition was: 
> cell_index == 0 || !at_boundary()
> Additional information: 
> You are on a boundary, so you can only ask for the first FEFaceValues 
object.

...if the cell is on a boundary, but the face is an interior face. 

I've attached code for a small test case illustrating the issue. It 
consists of a 2x2 2D mesh, and attempts to use FEInterfaceValues to get 
field values from the cell "above".

Is this indeed a bug, or am I misunderstanding the intended usage of 
FEInterfaceValues? I'm using deal.II, version 9.2.0. 

Best,
Corbin

-- 
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/22813e93-5959-4244-aa45-1652d3e8c7b8n%40googlegroups.com.
// empty integration test
//
#include 
#include 

// deal ii imports
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

using namespace dealii;

template
struct ScratchData
{
  FEValues fe_values;
  FEInterfaceValues fe_interface_values;
  std::vector> fe_local_support_on_face;

  ScratchData(const Mapping &mapping,
  const FiniteElement &fe,
  const Quadrature&cell_quadrature,
  const Quadrature  &face_quadrature,
  const UpdateFlags  cell_update_flags = update_values 
  | update_gradients 
  | update_quadrature_points 
  | update_JxW_values,
  const UpdateFlags interface_update_flags = update_values 
  | update_quadrature_points
  | update_JxW_values 
  | update_normal_vectors)
  : fe_values(mapping, fe, cell_quadrature, cell_update_flags)
  , fe_interface_values(mapping, fe, face_quadrature, interface_update_flags)
  , fe_local_support_on_face(GeometryInfo::faces_per_cell)
  {
  for (unsigned int face_no : GeometryInfo::face_indices())
for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
  {
if (fe.has_support_on_face(i, face_no))
  fe_local_support_on_face[face_no].push_back(i);
  }
  }
};

int main()
{
  
  const int dim = 2;

  // create an [NN]^2 grid
  std::vector repetitions{2,2};
  
  const Point bottom_left(0.0, 0.0);
  const Point   top_right(1.0, 1.0);

  Triangulation volume_mesh;
  GridGenerator::subdivided_hyper_rectangle(volume_mesh, repetitions, bottom_left, top_right);

  // DG elements over the mesh
  const unsigned int p_order = 3;
  FE_DGQ fe(p_order);
  DoFHandler  dof_handler(volume_mesh);
  dof_handler.distribute_dofs(fe);

  // interior field phi = 1, to be integrated
  Vector field;
  field.reinit(dof_handler.n_dofs());
  VectorTools::interpolate(dof_handler, ConstantFunction(1.0), field);

  // create FE values, interface values via scratchdata
  MappingQ1 mapping;
  const QGauss quadrature(fe.degree + 1);
  const QGauss face_quadrature(fe.degree + 1);
  ScratchData sd(mapping, fe, quadrature, face_quadrature);

  for (const auto &cell : dof_handler.active_cell_iterators())
  {

sd.fe_values.reinit(cell);

// face loop
for (const auto face_index : GeometryInfo::face_indices())
{
  sd.fe_interface_values.reinit(cell, face_index);
  const std::vector> &nrml = sd.fe_interface_values.get_fe_face_values(0).get_normal_vectors();

  // if top face (and not boundary face, just to be sure)
  if (nrml[0][1] > 0 && !(cell->face(face_index)->at_boundary()))
  {
std::cout << "cell " << cell << " at " << cell->center()[0] << " " << cell->center()[1] << std::endl;

const double center_x = cell->face(face_index)->center()(0);
const double center_y = cell->face(face_index)->center()(1);

std::cout << "face " << face_index << " at " << center_x << " " << center_y << std::endl;
std::vector field_at_qpts(sd.fe_interface_values.n_quadrature_points);

// this line fails, with the error message:
// >>> The violated condition was:  cell_index == 0 || !at_boundary()
// >>> Additional information: You are on a boundary, so you can only ask for the first FEFaceValues object.
// if I c

[deal.II] Multiple material IDs

[Corbin Foucart]

Hello all,

I've read the documentation on material ID and it's not clear to me if 
there's a way to define more than one material ID. This could be helpful 
when wanting to isolate elements that correspond to the same material, but 
also when looking to isolate elements of different materials that might 
exist in a particular "layer" of the domain. 

I can see three ways to define multiple IDs per each element:

   1. Define two identical meshes and set the material IDs differently on 
   each. This seems suboptimal because I don't really want to use up the extra 
   memory
   2. Define functions that re-label the material IDs based on the current 
   task at hand, at the expense of extra computation and logic
   3. Define an external vector of structs that contain the different 
   element IDs by cell, but then it's unclear how to filter elements by 
   MaterialIDEqualTo() without first relabeling as in (2)

Is there a standard way to do this that I'm missing? Any advice is much 
appreciated!

Best,
Corbin

-- 
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/a66dcb53-8a6b-46ce-8c29-90067829a0a3n%40googlegroups.com.

Re: [deal.II] Hashing / mapping Point locations to material ID

[Corbin Foucart]

Thank you for the fast and clear answer, Wolfgang.

Just add such a hash function for Point. It's true that two nearby points 
> will 
> have different hash values, but that shouldn't stop you from using hashes 
> if 
> you query the same points over and over.

This is definitely doable. 
 

> That said, I would suggest to do the indexing not on geometric locations 
> but 
> the logic position of a cell. You can either do that with a cell iterator 
> itself, or if you want to make it work in a parallel context, the CellId 
> of a 
> cell.

I'm curious about this; if the material ID physically represents a 
location-based quantity (e.g., which vertical column of fluid as specified 
by a surface mesh), how come you recommend using the iteration order to 
specify it? Unless I'm misunderstanding what you mean by logic position of 
the cell w/r/t a cell iterator...

Best,
Corbin 

-- 
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/fdf1fd09-49c7-485a-9faa-9cc6f77c2143n%40googlegroups.com.

[deal.II] Hashing / mapping Point locations to material ID

[Corbin Foucart]

Hello all,

I'm attempting to set material ID based on x-location. I'd like to create  
std::unordered_map between x-locations of cell centers and material id, and 
then refer to that map to quickly set material_id later.

The compiler is rightly complaining that there isn't a hash method for the 
Point object, and it seems like this could indeed be dangerous, because of 
floating point comparison later between Point objects during lookup.

However, this must be a commonly done sort of thing... what's the best way 
to map a known set of Point objects to other data?

Best,
Corbin

-- 
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/fd143114-2dec-44d2-9ebe-f5395289de10n%40googlegroups.com.

Re: [deal.II] Outputting auxiliary cell-level data along with solution data

[Corbin Foucart]

Hi Praveen,

Thanks for the helpful answer. I'm new to Visit; would you mind specifying 
how to view the subset of data (e.g., faces with boundary ID 0)? 

Best,
Corbin

On Thursday, November 15, 2018 at 9:31:16 PM UTC-5 Praveen C wrote:

> You can use this
>
> https://dealii.org/developer/doxygen/deal.II/structGridOutFlags_1_1Vtu.html
>
> https://dealii.org/developer/doxygen/deal.II/structDataOutBase_1_1VtkFlags.html
>
> to set solution time.
>
> Here is an example from my code
>
> https://github.com/cpraveen/dflo/blob/master/src_mpi/output.cc#L66
>
> To visualize boundary faces, VisIt has an option to make use of material 
> id of cells and faces, which allows you to view subset of the data, e.g., 
> plot only faces with material id = 100. You must save material id into 
> vtk/vtu file and then use that in VisIt. I have never done this through 
> deal.II but maybe there is a better way to do it in deal.II which I am not 
> aware of.
>
> Best
> praveen
>
> On 16-Nov-2018, at 3:49 AM, mrjonm...@gmail.com wrote:
>
> I have a model where I'm using the boundary_id to specify boundary 
> regions where certain boundary conditions apply, and I'd like to add the 
> boundary_id numbers to the output, so I can see them spatially in Paraview 
> along with my solution output. 
>
> I'm currently trying to use a DataOutFaces object along with a class 
> derived from DataPostprocessor to get these values and output them to .vtu 
> files. I'm having trouble debugging it, and I was curious whether there is 
> a simpler way to output data that are defined for an entire cell or face. 
>
> Along with that, it would also be helpful to be able to output some global 
> scalar data (e.g. simulation time) into the output files. If there's an 
> easy way to to that, I'd appreciate it as well. In my current approach, I'm 
> just going to set a variable named "time" to the appropriate value for all 
> DoFs, but it seems like an inefficient way to do it. 
>
> Thanks,
> Jonathan
>
> -- 
> 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.
> For more options, visit https://groups.google.com/d/optout.
>
>
>

-- 
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/a2b7d3c8-78f8-40bb-b384-c1c5ce62bc6an%40googlegroups.com.

[deal.II] Surface / volume interactions for a 3D extruded mesh

[Corbin Foucart]

Hello deal.ii community,

I've really enjoyed using deal.ii over the past few weeks -- it's very 
impressive how complete the software is. I'm in the process of implementing 
an ocean equations solver which has a few unusual features:

   1. The problem is to be solved over a 3D mesh extruded downwards in 
   depth from a 2D surface mesh
   2. There is a free surface equation which makes use of the original 2D 
   mesh
   3. 3D data is communicated to the 2D surface by integrating finite 
   element fields vertically over the entire domain in depth (depth-averaged 
   fields)
   4. From a 3D scalar field, I will need to compute another 3D scalar 
   field representing the original field integrated to depth z (like a 
   hydrostatic pressure at every point in the domain)

How I was thinking of doing this:

   - (1) Mesh the 2D region with gmsh and extrude it using the approach in 
   step-49; this seems clear
   - (2,3) I was considering holding the original 2D mesh in memory and 
   creating some type of map between the top-level volume elements and the 
   surface.* Is there a better way I could go about mapping the 3D volume 
   data and depth-averaged data to the surface mesh?*
   - (4) I can re-write the hydrostatic integrals as an ode in depth and 
   march element-by-element from surface to bottom along each vertical column 
   using a DG scheme. 
  - To isolate each vertical column, I could assign a material ID to 
  each one
  - However, I need to extract the value on the bottom of each element 
  above as a boundary condition to the element below
  - I could probably do this with FEFaceValues, *but it's unclear to me 
  how to "request" values on a top or bottom face* (since element 
  orientation is general). Is there an ordering I could rely on if I know 
it 
  for a single element?
   
Any tips would be appreciated!

Best,
Corbin

-- 
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/55879d8b-7ca3-4abf-b59f-9eee8b49be25n%40googlegroups.com.

Re: [deal.II] Anisotropic diffusion implementation

[Corbin Foucart]

It seems my images were not rendered -- here they are in order as
attachments.

On Thu, Dec 17, 2020 at 10:56 PM Corbin Foucart 
wrote:

> I have a problem with direction-dependent diffusion (that is, the
> diffusivity \nu depends on the derivative direction but not spatial
> location) such that the horizontal diffusivity and vertical diffusivity are
> different, leading to a Laplacian as follows:
>
>
> I am discretizing the problem with a mixed method, which defines the
> auxiliary variable L = ∇u such that
>
>
> So in my weak form, I end up with terms like the following:
> (please excuse the abuse of index notation with repeated j three times).
> phi_ij here refers to each scalar basis function in the tensor-valued
> finite element space.
>
> It seems using the deal.ii FEValuesExtractors, I can access the divergence
> of each tensor-valued basis function, but I need to multiply each term in
> the divergence with the appropriate value of \nu by component.
>
> I suppose I could ask for the gradient of each phi and do the sum
> manually,but I was wondering if there's a way to do it directly. Is there a
> standard way of doing this in deal.ii?
>
> Any pointers would be appreciated!
>
> --
> 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/ptfmuBKUoiA/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to
> dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/65594873-33d7-44c2-834b-ba398fd9fb1bn%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/65594873-33d7-44c2-834b-ba398fd9fb1bn%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/CAE%2BMPodoB41SAR8UrJHq_cNDhCqed1cKxtQNr%3DZ-Q5UBgk4syQ%40mail.gmail.com.

[deal.II] Anisotropic diffusion implementation

[Corbin Foucart]

I have a problem with direction-dependent diffusion (that is, the 
diffusivity \nu depends on the derivative direction but not spatial 
location) such that the horizontal diffusivity and vertical diffusivity are 
different, leading to a Laplacian as follows:


I am discretizing the problem with a mixed method, which defines the 
auxiliary variable L = ∇u such that


So in my weak form, I end up with terms like the following:
(please excuse the abuse of index notation with repeated j three times). 
phi_ij here refers to each scalar basis function in the tensor-valued 
finite element space.

It seems using the deal.ii FEValuesExtractors, I can access the divergence 
of each tensor-valued basis function, but I need to multiply each term in 
the divergence with the appropriate value of \nu by component.

I suppose I could ask for the gradient of each phi and do the sum 
manually,but I was wondering if there's a way to do it directly. Is there a 
standard way of doing this in deal.ii?

Any pointers would be appreciated!

-- 
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/65594873-33d7-44c2-834b-ba398fd9fb1bn%40googlegroups.com.

Re: [deal.II] Re: Interacting with Python data // external codes

[Corbin Foucart]

Thank you for the advice and help! 

Now that I'm a bit further along in my project, I think I can ask some more 
precise questions. The code will have to be driven in Python, but hopefully 
I can find a good way to access the deal.ii data structures without writing 
to file. I'm using a time-dependent HDG scheme on an adaptive grid. From a 
python session, is it possible to access:

   1. The state of the triangulation ( this seems doable from the python 
   bindings, as in tutorial 1 in the notebooks Bruno linked)
   2. For every cell in the mesh
  1. the nodal DG coefficients u_h
  2. the solution jumps [[u_h]] on all faces of the cell (between the 
  solution on this cell and its neighbors)
  3. the right hand side (source term and time-dependent forcing), 
  ideally at the nodal points of the cell
   3. The coarsen or refine flags for the triangulation (seems possible as 
   in 1)?

I think the best way to do this would be with the Cython approach suggested 
by Alex, however, it's unclear to me the best way to pass the cell-based 
data to Python. Is it possible to define a struct on from deal.ii 
containing the data in (2) and pass an array of these structs (one per 
cell) to the python session? 

Thank you,
Corbin
 

On Tuesday, November 24, 2020 at 11:02:21 PM UTC-5 Alex Cobb wrote:

> Hi Corbin,
>
> On Wed, 25 Nov 2020 at 00:05, Bruno Turcksin  wrote:
>
>> deal.II has some limited support for python mainly for mesh manipulation. 
>> We have some python notebooks here 
>> .
>>  
>> I think what you want to do is similar to the step-62 notebook. Right now, 
>> the only way to interact with numpy is to print the data to a file and then 
>> load it (see here 
>> 
>>  
>> and here 
>> ).
>>  
>> If you want to manipulate the mesh directly in python, you need 
>> boost.python and you need to configure deal.II with 
>> -DDEAL_II_COMPONENT_PYTHON_BINDING=ON. It's sometimes a little bit tricky 
>> to enable the python binding so don't hesitate to ask any question on the 
>> mailing list if you need help.
>>
>> On Tuesday, November 24, 2020 at 12:16:18 AM UTC-5 corbin@gmail.com 
>> wrote:
>>
>>>
>>>- what is the best practice for exporting deal.ii solution data in a 
>>>way that Python / numpy can interact with it?
>>>
>>> Depending on what you want from the data, another option might be to 
> export the data as .vtu (*not* .vtk) and read it with VTKPython.  For 
> example, you can use VTKPython to grid it, and then work with the gridded 
> data with numpy.  
>
> If I/O turns out to be a bottleneck, you could also consider writing the 
> data as HDF5
> https://www.dealii.org/current/doxygen/deal.II/namespaceHDF5.html
> and then accessing the HDF5 file from Python (using h5py, for example).  I 
> haven't used deal.II's HDF5 export but HDF5 can sometimes vastly improve 
> I/O performance compared to text files or even other binary formats (e.g., 
> with blosc compression).
>
>>
>>>-  Is there a good way for external software to 'hook' into the 
>>>deal.ii pipeline? Something like:
>>>   - initialize a triangulation / grid
>>>   - run the solver
>>>   - make a call like: new_data = external_software(deal_ii_output, 
>>>   grid)
>>>   - reinitialize the grid based on new_data
>>>   - loop 
>>>
>>> Have you used Cython at all?
> https://cython.org/
> It is my favorite way to use C and C++ libraries from Python (compared to 
> binding generators, or writing extensions by hand, or cffi).  It provides a 
> very natural way to transfer data to and from C or C++ library code using 
> typed numpy arrays.  It can take some patience to get what you want from 
> the documentation, but if you are already used to Python and have some 
> familiarity with C and C++ it makes it very easy to migrate code between C 
> / C++ and Python.  So, for example, you can prototype in Python and 
> gradually shift functionality over to C or C++.
>
> What I have found easiest is to put the heavy-lifting code in a static C++ 
> library, building from the (awesome) deal.II tutorials and documentation, 
> and then wrap that C++ library into a loadable Python extension module 
> using Cython.  Then you can pass arguments from Python to your solver using 
> the extension module.
>
> In my experience, the trickiest part of this is not Cython per se, but 
> getting testing and continuous integration working with the mix of 
> languages: C++, Cython, Python, and your build system (CMake?) 
> mini-language.
>
> If you choose this route, another option for controlling your solver from 
> Python (possibly with less pain / dependencies than Boost.Python?) might be 
> pybind11 (disclaime

[deal.II] Interacting with Python data // external codes

[Corbin Foucart]

Hi everyone,

I am new to deal.ii but after reading through the documentation and the 
tutorials I'm very excited about all the functionality that the software 
offers! I'm aiming to link a finite element solver with a machine learning 
code that I have in python. To that end:

   - what is the best practice for exporting deal.ii solution data in a way 
   that Python / numpy can interact with it? 
   - Is there a good way for external software to 'hook' into the deal.ii 
   pipeline? Something like:
  - initialize a triangulation / grid
  - run the solver
  - make a call like: new_data = external_software(deal_ii_output, grid)
  - reinitialize the grid based on new_data
  - loop 
   
Apologies if this is explained elsewhere!

Corbin


-- 
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/4ffb7687-8af5-4c90-ba09-246d0c20dac1n%40googlegroups.com.