[deal.II] MPI: Evaluating several distributed functions on each node

[Konrad Simon]

Dear all,

I am facing a little challenge: Suppose I have several objects myFunction 
of a certain class MyFunction derived from a Function object on 
each process in my MPI universe. The number of MyFunction objects may vary 
from process to process.
I am trying to evaluate each object at certain points p that live on a 
distributed::Triangualtion. The evaluation points are different for 
each MyFunction object. A point may be contained in a cell owned by another 
process whose id is known, so I know where to send it to to ask a 
MyFunction object on that specific process to evaluate it for me and send 
the value back.

Now my question: How do I deal with several MyFunction objects on one 
process (I can identify these by a hash). Is it a good idea to have a class 
with a set of static variables and methods on each node to manage the 
communication and distribute necessary information within the process using 
the specific MyFunction hash? Does that make sense?

Let me stress again, the difficulty here is that I have several objects 
that may ask for values one by one or even in a (shared memory) parallel 
manner that may or may not trigger communication (sort of some_to_some 
style as implemented in Deal.II). 

Thanks in advance :-)

Best,
Konrad

-- 
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/4fb6b294-aa7f-4dc2-b823-e42247bf2302n%40googlegroups.com.

[deal.II] Triangulation from lower-dimensional object

[Konrad Simon]

Dear all,

I was recently looking for a solution to the following problem:

I have a triangulation, e.g., of a unit cube. I need to solve a 
lower-dimensional (FEM-like) problem on all faces with BCs on edges and for 
each face I need to solve a problem on all edges with BCs at vertices.

Can I somehow extract a Triangulation form a Triangulation by 
providing, for example, the boundary id on the Triangulation level?

Best,
Konrad

-- 
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/ae0eae00-a572-4036-b20e-31f45dd70892n%40googlegroups.com.

[deal.II] Elasticity Preconditioner

[Konrad Simon]

Dear all, 

A student and me are trying to deal with (parallel, static) linearized 
elasticity similar to step-8.
However, the problem is slightly different since the Lame parameters are 
functions instead of constants. 
Can anyone recommend good solvers & preconditioners for the resulting 
system? Our boundary conditions fix all degrees of freedom in the kernel.

Best,
Konrad

-- 
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/eba0cc0f-6028-4fde-9cee-307e09f8d2ban%40googlegroups.com.

[deal.II] Re: Question about deal.II flexibility

[Konrad Simon]

Dear Abdo,

I do think that Deal.II can help you with solving your problem. If you are 
just starting with Deal.II have a look at the tutorials 
. It is always 
a good idea to start with some code basis that other people have shared and 
that was thoroughly tested. As to your future problem, please have a look 
at the code gallery 
 as well. 
There are codes that other users shared dealing with visco-plastic 
materials, in particular (but not only Jean-Paul).

Cheers,
Konrad

On Saturday, June 19, 2021 at 10:23:37 PM UTC+2 aael...@ncsu.edu wrote:

> Hi
> I am Abdo, Ph.D. student at NC state university and was looking for a FE 
> package library that is flexible and can be used in modeling and analyzing 
> 2D and 3D elastic wave equation propagation within an incompressible 
> medium. Also, at some point, I might need to add the viscoelastic behavior 
> to the model. So, I am not sure if deal.II is flexible enough to add, 
> modify some feature in it or not such that I can do what I have mentioned 
> earlier. So, Can you please give me some pointers or advice regarding 
> this?, your help is much appreciated.
>
> Best regards.
> Abdo.
>

-- 
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/481afb9d-fa92-4eae-b369-0881cbd7d70fn%40googlegroups.com.

Re: [deal.II] p4est user pointer and distributed::triangulation

[Konrad Simon]

Thank you, Wolfgang! This is exactly what I am looking for. I don‘t need to 
access the pointer from outside, I am trying to extend the p4est-interface 
itself.

Best,
Konrad

On Friday, May 14, 2021 at 12:42:48 AM UTC+2 Wolfgang Bangerth wrote:

> On 5/13/21 3:05 PM, Konrad Simon wrote:
> > 
> > I am currently writing an interface to some p4est functions. The 
> structure of 
> > p4est makes it often necessary to pass data around through a user 
> pointer of 
> > whatever type (void*). When creating a new triangulation this pointer is 
> set 
> > to „this“ (the triangulation itself), see for example
> > 
> > 
> https://www.dealii.org/current/doxygen/deal.II/distributed_2tria_8cc_source.html#l2989
>  
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2Fdistributed_2tria_8cc_source.html%23l2989=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cb4cadd28da84445fd29308d91652e2c3%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637565367510208994%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000=v2G100jGfI6SbGwV35raSFjMxfnJMJurP2F2qAc61Ls%3D=0
> >
> > 
> > Is this pointer used in other interfaces somehow? p4est itself does not 
> touch 
> > it so I am wondering if I can reset it to anything I‘d like.
>
> Konrad,
> we use this user pointer in many of the functions that are called back 
> from 
> p4est. Take a look at
> RefineAndCoarsenList::refine_callback
> for example (line 761 of the file), and many other functions that you can 
> find 
> by searching for
> forest->user_pointer
> in this file.
>
> We consider the p4est forest object as internal to the p::d::T class, so I 
> would expect that you can't access it, or expect that you can use any of 
> its 
> content for anything other than what the p::d::T class does with it.
>
> 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/7fbd1adf-c6d4-4656-ab6e-7c833e5e92cen%40googlegroups.com.

[deal.II] p4est user pointer and distributed::triangulation

[Konrad Simon]

Dear all, 

I am currently writing an interface to some p4est functions. The structure 
of p4est makes it often necessary to pass data around through a user 
pointer of whatever type (void*). When creating a new triangulation this 
pointer is set to „this“ (the triangulation itself), see for example

https://www.dealii.org/current/doxygen/deal.II/distributed_2tria_8cc_source.html#l2989

Is this pointer used in other interfaces somehow? p4est itself does not 
touch it so I am wondering if I can reset it to anything I‘d like.

Any Ideas?

Best,
Konrad

-- 
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/f23950cf-7688-45d2-b986-140ccc3a58f9n%40googlegroups.com.

Re: [deal.II] interpolate FE_Nelelec

[Konrad Simon]

Hi John,

Maybe I can give some hints of how I would approach this problem (these are 
just some quick thoughts):
Write this problem as a vector Laplace problem

curl(nu*curl(A)) - grad(div(A)) = J

which you then have to write as a system of PDEs. Note, this can only be 
the same as

curl(nu*curl(A)) = J

if J is a curl itself.

If you do not want to impose conditions on the divergence then you have two 
options:

1.) You use sigma = div(A) as an auxiliary variable and impose natural 
boundary conditions (i.e., the ones that you get through integration by 
parts) on A*n and nu*curl(A)xn, where n is the outward unit normal. This 
means that sigma is a 0-form and A is a 1-form -> sigma is then discretized 
(for example) by FE_Q(n+1) and A by FE_Nedelec(n)

System:

(A1) sigma + div(A) = 0
(B1) grad(sigma) + curl(nu*curl(A)) = J

For the weak form you multiply (A1) with a FE_Q test function and integrate 
the second summand by parts. Secondly, you multiply (B1) with a Nedelec 
test function and also integrate the second summand by parts. You will see 
your natural boundary conditions pop up.

2.) You use sigma = nu*curl(A) as an auxiliary variable and impose 
essential boundary conditions on A*n and nu*curl(A)xn, where n is the 
outward unit normal. This means that sigma is interpreted as a 1-form and A 
is a 2-form -> sigma is then discretized (for example) by FE_Nedelec(n) and 
A by FE_RaviartThomas(n)

System:

(A2) (1/nu)*sigma - curl(A) = 0
(B2) curl(sigma) -  grad(div(A)) = J

For the weak form you multiply (A2) with a Nedelec test function and 
integrate the second summand by parts. Secondly, you multiply (B2)
with a Raviart-Thomas test function and also integrate the second summand 
by parts. You will NOT see your natural boundary conditions pop up since 
conditions on A*n and nu*curl(A)xn are essential in this case. You need to 
enforce them on the function spaces directly. In deal.ii you do so by using 
project_boundary_values_curl_conforming_l2 

 
or project_boundary_values_div_conforming 


Note that the additional boundary conditions make the system invertible and 
if J is a curl it will turn out that div(A)=0.

There is a field in mathematics that is called finite element exterior 
calculus (FEEC) that answers the question of stability when solving such 
problems. See this book 
. I guess Chapters 
4, 5 and 6 are most interesting for you.

> I do not understand when a geometry gets complicated. Is a toroid inside a 
> sphere, both centred at the origin, a complicated geometry? To start with, 
> I will use a relatively simple mesh. I will not use local refinement, only 
> global. I can do without refinement as well, i.e. making the mesh in gmsh.
>
> Yes, I would like to know more about orientations issues on complicated 
> geometries. Could please tell me about it? I would very much prefer to use 
> FE_Nedelec without hierarchical shape functions. The first order 3D 
> FE_Nedelec will do nicely provided the orientation issues will be 
> irrelevant to the mesh.
>
 Lowest order Nedelec and all other elements I mentioned are fine on the 
meshes you need I guess.

Hope that helps.
Best,
Konrad  

-- 
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/c2773841-4e03-42a1-9f1e-e1a41c0933c9n%40googlegroups.com.

Re: [deal.II] interpolate FE_Nelelec

[Konrad Simon]

Dear John,

As Jean-Paul pointed out the entries of the solution vector for a Nedelec 
element have a different meaning. The entries are weights that result from 
edge moments (and in case of higher order also face moments and volume 
moments), i.e, integrals on edges. Nedelec elements have a deep geometric 
meaning: they discretize quantities that you can integrate (a physicist 
would say "measure") along 1D-lines in 3D.

Btw, if you are using Nedelec in 2d be aware that it has some difficulties 
on complicated meshes. In 3D the lowest order is fine. NedelecSZ (for 2D 
lowest order and 3D all orders), as Jean-Paul pointed out, are another 
option. They are meant to by-pass certain mesh orientation issues on 
complicated geometries (I can tell you a bit about that since currently I 
am chasing some problems there).

Best,
Konrad

On Tuesday, April 20, 2021 at 9:23:23 PM UTC+2 Jean-Paul Pelteret wrote:

> Hi John,
>
> Unfortunately, you’ve fallen into the trap of confusing what the entries 
> in the solution vector mean for the different element types. For Nedelec 
> elements, they really are coefficients of a polynomial function, so you 
> can’t simply set each coefficient to 1 to visualise the shape function 
> (like one might be inclined to try for Lagrange type polynomials). If you 
> want a little piece of code that will plot the Nedelec shape functions for 
> you, then take a look over here:
> https://github.com/dealii/dealii/issues/8634#issuecomment-595715677
>
> BTW. In case you’re not aware, a newer (and probably more robust) 
> implementation of a Nedelec element is the FE_NedelecSZ 
> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> that 
> uses hierarchical shape functions to form the higher order bases. The PhD 
> thesis of S. Zaglmayr that describes the element is mentioned in the class 
> documentation.  
>
> Best,
> Jean-Paul
>
> On 20. Apr 2021, at 15:47, John Smith  wrote:
>
>
> Dear Konrad,
>
> Thank you for your help! I have recently modified the step-3 tutorial 
> program to save FE_Nedelec<2> fe(0) shape function into *.csv files. It is 
> attached to this message. The result of this program is not what I would 
> expect it to be. I have also attached a pdf file which compares my 
> expectations (marked as “Expectations” in the pdf file) with the results of 
> the modified step-3 program (marked as “Reality” in the pdf file). My 
> expectations are based on the literature I have. My question is: why the 
> shape functions of FE_Nedelec<2> fe(0) differ from my expectations? I guess 
> my expectations are all wrong. I would also appreciate if you will share 
> with me a reference to a paper or a book, so I can extend my modest library.
>
> Thanks!
>
> John
>
> P.S. I have sent you an e-mail as you have requested. Could you please 
> check your spam folder? 
>
>
>
> On Monday, April 19, 2021 at 8:47:29 AM UTC+2 Konrad Simon wrote:
>
>> Dear John,
>>
>> I had a look at the pdf you sent. I noticed some conceptual 
>> inconsistencies that are important for the discretization. Let me start 
>> like this: The curl-curl-problem that you are trying to solve is - 
>> according to your description of the gauge - actually a curl-curl + 
>> grad-div problem and hence a Laplace-problem that describes a 
>> Helmholtz-decomposition. In other words in order to get a well-posed 
>> simulation problem you need more boundary conditions. You already noticed 
>> that since your curl-curl-matrix is singular (the kernel is quite 
>> complicated and contains all longitudinal waves). 
>>
>> You have of course the option to solve (0.0.2) with a Krylov-solver but 
>> then you need to make sure that the right-hand side is in the orthogonal 
>> complement of this kernel before each iteration which is quite difficult. I 
>> would not recommend that option.
>>
>> Another point is that if you have a solution for your field A like in 
>> (0.0.4) you can not have a similar representation for T. This is 
>> mathematically not possible.
>>
>> Since you not care about the gauge let me tell you how I would tackle 
>> this: The reason  is  - to come back to my original point - you are missing 
>> a boundary condition that makes you gauge unique. Since you apply natural 
>> boundary conditions (BCs) on A you must do so as well to determine div(A). 
>> This second BC for A is applied to a different part of the system that is 
>> usually neglected in the literature and A is partly determined from this 
>> part (this part then describes the missing transversal waves). The 
>> conditions you mention on curl(A) are some what contradictory to the space 
>> in whic

Re: [deal.II] Re: interpolate FE_Nelelec

[Konrad Simon]

Dear John,

I had a look at the pdf you sent. I noticed some conceptual inconsistencies 
that are important for the discretization. Let me start like this: The 
curl-curl-problem that you are trying to solve is - according to your 
description of the gauge - actually a curl-curl + grad-div problem and 
hence a Laplace-problem that describes a Helmholtz-decomposition. In other 
words in order to get a well-posed simulation problem you need more 
boundary conditions. You already noticed that since your curl-curl-matrix 
is singular (the kernel is quite complicated and contains all longitudinal 
waves). 

You have of course the option to solve (0.0.2) with a Krylov-solver but 
then you need to make sure that the right-hand side is in the orthogonal 
complement of this kernel before each iteration which is quite difficult. I 
would not recommend that option.

Another point is that if you have a solution for your field A like in 
(0.0.4) you can not have a similar representation for T. This is 
mathematically not possible.

Since you not care about the gauge let me tell you how I would tackle this: 
The reason  is  - to come back to my original point - you are missing a 
boundary condition that makes you gauge unique. Since you apply natural 
boundary conditions (BCs) on A you must do so as well to determine div(A). 
This second BC for A is applied to a different part of the system that is 
usually neglected in the literature and A is partly determined from this 
part (this part then describes the missing transversal waves). The 
conditions you mention on curl(A) are some what contradictory to the space 
in which A lives. Nedelec elements which you must use (you can not use FE_Q 
to enforce the conditions) cannot generate your desired T_0.

There are some principles when discretizing these problems (which are not 
obvious) that you MUST adhere to (choice of finite elements, boundary 
conditions, what is the exact system etc) if you want to get a stable 
solution and these are only understood recently. I am solving very similar 
problems (with Deal.II) in a fluid context and will be happy to share my 
experiences with you. Just email me: konrad.si...@uni-hamburg.de

Best,
Konrad

On Monday, April 19, 2021 at 5:51:17 AM UTC+2 Wolfgang Bangerth wrote:

> On 4/18/21 12:24 PM, John Smith wrote:
> > 
> > Thank you for your reply. It is a bit difficult to read formulas in 
> text. So I 
> > have put a few questions I have in a pdf file. Formulas are better 
> there. It 
> > is attached to this message. May I ask you to have a look at it?
>
> John:
>
> If I understand your first question right, then you are given a vector 
> field J 
> and you are looking for a vector field T so that
> curl T = J
> I don't really have anything to offer to this kind of question. There are 
> many 
> vector fields T that satisfy this because of the large null space of curl. 
> You 
> have a number of condition you would like to "approximate" but there are 
> many 
> ways to "approximate" something. In essence, you have two goals: To 
> satisfy 
> the equation above and to approximate something. You have to come up with 
> a 
> way to weigh these two goals. For example, you could look to minimize
> F(T) = ||curl T - J||^2 + alpha ||T-T_desired||^2
> where T_desired is the right hand side of (0.0.5) and you have to 
> determine 
> what alpha should be.
>
> As for your other questions: (0.0.9) is indeed called "interpolation"
>
> The issue with the Nedelec element (as opposed to a Q(k)^dim) field is 
> that 
> for the Nedelec element,
> \vec phi(x_j)
> is not independent of
> \vec phi(x_k)
> and so you can't choose the matrix in (0.0.8) as the identity matrix. You 
> realize this in (0.0.13). The point I was trying to make is that (0.0.13) 
> cannot be exactly satisfied if T(x_j) is not a vector parallel to \hat 
> e_j, 
> which in general it will not be. You have to come up with a different 
> notion 
> of what it means to solve (0.0.13) because you cannot expect the left and 
> right hand side to be equal.
>
> 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/75104172-4893-4ae6-af8e-8ec2c12d335an%40googlegroups.com.

Re: [deal.II] 2d meshes and orientation

[Konrad Simon]

Dear Wolfgang,

On Tuesday, April 6, 2021 at 5:53:15 AM UTC+2 Wolfgang Bangerth wrote:

> It's possible this is indeed a bug. In most cases, we run the 2d meshes 
> through the mesh orienter, so we rarely see 2d meshes that are not 
> correctly 
> oriented and it wouldn't surprise me if we have the kind of bug you 
> describe 
> where things almost but not quite work correctly.


I see. Many thanks for the input. This is my impression as well. These 
things rarely happen in 2d.
One example would be simple periodic meshes which cause line orientations 
not to match (this case is caught). In more
complicated scenarios involving periodicity and rotations of faces it can 
happen though that lines do match but normals do not.


> The questions then are (i) where does this need to be fixed, and (ii) how 
> many 
> tests fail if we apply the fix? One of the things I learned (and I think 
> you 
> did as well) is that for bugs in these kinds of layers, it is quite 
> possible 
> that one ends up with 100 failing tests and then it takes a *lot* of work 
> to 
> look through all of these to figure out whether the new output is correct 
> or 
> not. But then, as mentioned above, it's also possible that only one or two 
> tests fail because we almost never see these kinds of meshes.
>

Sound worth trying. Please see also this new PR 
 which is work in progress.

-- 
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/0125f6dc-b61c-4db6-84da-b101b6759f1bn%40googlegroups.com.

[deal.II] 2d meshes and orientation

[Konrad Simon]

Dear deal.ii community,

I stumbled over something in 2d meshes that confuses me a bit. Supoose I 
have 2 squares (a left and a right one) sharing one edge. If everything is 
fine all lines and faces have standard orientation, i.e., 
line_orientation==true and face_orientation==true for all lines and faces. 
Note that lines and faces are NOT the same geometric entity here since 
their orientation is determined differently (faces have normals and lines 
have tangents).

Now, if I rotate the right square by 90 degrees I get (on the shared edge) 
non-matching face normals but all line orientations do match.

If I rotate the right square by 180 degrees I get (on the shared edge) 
non-matching face normals and non-matching line orientations.

If I rotate the right square by 270 degrees I get (on the shared edge) 
matching face normals and non-matching line orientations.

This is at least my intuition. When I query these information for each 
rotation case I see that my intuition is confirmed for lines but the face 
orientations are always true.

Is this a bug?

I am asking because it seems that 2d RaviartThomas elements (which make use 
of face normals) are broken for some edge cases (3d RT I could fix).

Best,
Konrad

-- 
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/f739c450-c37d-4d24-bf5b-2778dec5d501n%40googlegroups.com.

[deal.II] Re: Need help installing deal.ii with p4est and Trilinos

[Konrad Simon]

Dear Budhyant,, 

Hard, to say what should be the paths on your system. Dealing with all 
these large dependencies can be quite cumbersome. I would say the way you 
install Deal.II depends on what you want to do with it. Let me tell you how 
I would approach this.

1. If you just want to develop a parallel application (using all the 
packages that you mention) on your local laptop it is in my opinion usually 
enough to use the packages that are provided by the OS package manager. 
Simply use apt install if it is your own machine and you have super user 
permissions. Then cmake should find these packages in the PATH and you just 
need to tell cmake to use them through, e.g., DEAL_II_WITH_P4EST = ON etc, 
for the build of Deal.II. This is usually sufficient for development since 
for running large 3D problems a single machine is anyway not enough. So no 
need to optimize a build for speed

2. If you do not have root permissions (or if you want to do some larger 
computations on a cluster, or you don't want to change the system packages) 
installing Deal.II using spack  is quite a nice option. 
Spack will resolve all dependencies, optimize the builds for your system 
and can be customized nicely (you can also link against MPI libraries and 
other stuff provided by the target system's admin). This is the option that 
I use.

Hope that helps.
Best regards,
Konrad

On Friday, February 19, 2021 at 10:07:00 AM UTC+1 venepalli...@gmail.com 
wrote:

> Dear All:
>
> I would appreciate it if you could please help me with the installation on 
> Ubuntu 20. 
>
> I have installed p4est, petsc, Trilinos and metis in "~/src/", they are 
> working individually. 
>
> Cmake succeeds in recognizing petsc and metis, however, it fails to 
> recognize the p4est and Trilinos installations. 
>
> *cmake -DDEAL_II_WITH_MPI=ON -DP4EST_DIR=~/src/p4est/build-mpi/ 
> -DPETSC_DIR=~/src/petsc-3.14.4/ -DPETSC_ARCH=amd64 
> -DTRILINOS_DIR=~/src/Trilinos/MY_BUILD/ -DMETIS_DIR=~/src/metis-5.1.0/ 
> -DCMAKE_INSTALL_PREFIX=~/src/dealii-9.2.0/ ..*
>
> What should the path be in my case?
>
> Thank you.
>
> Regards,
> Budhyant Venepalli
>

-- 
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/af5964fe-cfb1-4be7-a679-0235cdbf2b08n%40googlegroups.com.

Re: [deal.II] Face and line orientation dependent dof permutations in FE_Q

[Konrad Simon]

Dear Wolfgang, 
 

> Fundamentally, the place to look is to search through the finite element 
> classes where they ask for face and line orientations. It is possible that 
> you 
> don't actually have to do much there: Finite element classes generally 
> work in 
> the local (cell) coordinate system, where you enumerate shape functions as 
> 0...dofs_per_cell. So when, for example, you compute the local matrix, it 
> really is local. The interesting part happens when you translate what 
> *global* 
> degrees of freedom correspond to the local ones, i.e., when you call 
> cell->get_dof_indices() (==DoFCellAccessor::get_dof_indices). This is 
> where on 
> every cell has to collect which DoFs live on it, and that has to take into 
> account the orientations of lines and faces come into play. The cell DoF 
> indices are cached, but you can see how this works in 
> dof_accessor.templates.h, starting in line 902 using the 
> process_dof_indices() 
> functions. 
>

Yes, I found that function. It requires though, that two certain tables
 `adjust_quad_dof_index_for_face_orientation_table` and 
`adjust_line_dof_index_for_face_orientation_table`
are implemented in derived classes. This is the case for FE_Q but not for 
any element derived from FE_PolyTensor.
For the DG elemements this is not problematic since all dofs are interior. 
But for H(div) conforming elements
(FE_ABF, FE_RaviartThomas, FE_BDM) and for  H(curl), which is essentially 
only FE_Nedelec, it is problematic.

The issue is even worse than just the signs: Dofs are permuted for higher 
order and signs can switch depending on the permutation.

 (NedelecSZ btw does not inherit from FE_PolyTensor and has a dynamic way 
to fix the orientation issue but can not be used adaptively.)


> What this means is that if you have two cells that disagree about the 
> relative 
> orientation of their shared face, they when you call 
> cell->get_dof_indices() 
> on these two cells, you'd end up with a different ordering of the DoFs on 
> the 
> face. My usual approach to check this would be to create a mesh with two 
> cells 
> so that one has a non-standard face orientation, assign a Q3 element to 
> the 
> DoFHandler, and check that that is really what happens. (I generally turn 
> these little programs into tests for the test suite, just so that what I 
> just 
> verified to be correct -- or not! -- is actually checked automatically in 
> the 
> future. I would very much love to see you submit such tests as well!) 
>

This is exactly what I was doing and exactly the way I tested it. :-) 

>
> So this is the process for FE_Q, where there is really not very much to 
> do, or 
> in fact for all interpolatory elements. The situation becomes a bit more 
> difficult for things such as the Nedelec element (also the Raviart-Thomas 
> one) 
> where the shape functions are tangential to the edge. Let's say you have 
> the 
> lowest order Nedelec element. Then there is just one DoF on each edge, and 
> both adjacent cells agree about this. But they also have to agree *which 
> direction the vector-valued shape function should point to*, and that's 
> trickier. I can't remember how we do this (or whether we actually do this 
> at 
> all). I *think* that the RT element gets this right in 2d, but I don't 
> remember about the 3d case. Again, it would probably be quite helpful to 
> write 
> little programs and turn them into tests. 
>
> I know or strongly suspect that there are bugs for RT and Nedelec elements 
> on 
> meshes where faces have non-standard orientation. I've been thinking for 
> years 
> how I would go about debugging this, but I've never had the time to 
> actually 
> do it. But I'll share what I think would be the tools to do this: 
>

I was already working on a static fix. First, I fixed RaviartThomas. BDM 
and ABF can be taken care of in a similar static fashion.
For this I implement the above mentioned tables puls an additional table 
for the sign switches. I needed to change the implementation
of all elements but I do not break legacy behavior (apart form the fixed 
RaviartThomas) or the interface.

I actually got help from two developers/maintainers, see PR #11414 
.
 

> 1/ Creating meshes with non-standard orientations: You can feed 
> Triangulation::create_triangulation() directly so that the faces have 
> non-standard orientations with regard to each other. (See step-14 for how 
> to 
> call this function.) In 2d, this would work as follows: Think of a mesh 
> like this: 
>
> 3---4---5 
> | | | 
> 0---1---2 
>
> You'd describe the two cells via their vertex indices as 
> 0 1 3 4 
> 1 2 4 5 
> but if you want a mismatched edge orientation, you'd just use 
> 0 1 3 4 (thinks the edge is 1->4) 
> 5 4 2 1 (thinks the edge is 4->1) 
> or any other rotation. A good test would probably just try all 2x4 
> possible 
> rotations of both cells. 
>

Yes, I submitted another pull request containing exactly this.
 

>
> 

[deal.II] Face and line orientation dependent dof permutations in FE_Q

[Konrad Simon]

Dear Deal.II community,

I am currently fixing an orientation issue for some (non-primitive) vector 
elements with the help of two other Deal.II developers/maintainers.

I have a quick question: On complicated meshes some faces can have a 
non-standard orientation or they can be rotated against the neighboring 
face. This causes inconsistencies for some vector elements.

In principle this potentially causes inconsistencies for FE_Q elements as 
well (for all dofs located at vertices, lines and faces). But that is being 
taken care of in the library, so everything seems fine.

In order to understand the issue better for other elements: where in the 
library is this being taken care of? I have difficulties to find that. 

Remark: I found a table in (fe.h/fe.cc) that derived elements need to 
implement that takes care of local dof permutations on quads (not full 
faces and hence not line dofs or vertex dofs). Now, lines can have 
different orientations but lines also permute within a cell upon face 
orientation/rotation changes.

Can anyone point me to the place in Deal.II where this is being taken care 
of for FE_Q?

Best,
Konrad

-- 
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/7c3e0787-1b1c-465e-9747-3f82ff677d95n%40googlegroups.com.

Re: [deal.II] Mean value of component of TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

>
> Then you need to create a vector with locally relevant entries as ghosts 
> and 
> copy your fully distributed vector to it. There is no other way if you 
> want to 
> use that function. 
>
> But if your goal is to fix the pressure in a Stokes problem, it doesn't 
> have 
> to be the integral mean of the pressure that is zero. It could just be the 
> arithmetic mean, which you can compute without access to ghost elements.  
>

Many  thanks, Wolfgang, that works! :-)

Konrad

-- 
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/9849962a-5b13-4fb3-bcb5-bc74e045f2a5n%40googlegroups.com.

Re: [deal.II] Re: Mean value of component of TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

> That looks like an unrelated error. Can you create a small testcase for 
> this 
> issue here? 
>

I will try to come up with an example. 

Thank you again and best regards,
Konrad

-- 
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/83b966c1-4e3a-477c-827a-9d17f13682d8n%40googlegroups.com.

Re: [deal.II] Mean value of component of TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

>
> I suspect you are passing a fully distributed vector to that function, but 
> it 
> needs to read ghost elements of the vector. Have you tried copying the 
> vector 
> into a locally_relevant vector, and passing that to the function in 
> question? 
>

Thank you, Wolfgang, that was the issue.  I am using the mean value 
function at several different places in my code. One of them is inside a 
preconditioner vmult that accepts only a fully distributed vector. Is there 
also a simple workaround? I have a patch since I only need the pressure 
which is the last component of the distributed vector in each locally owned 
cell but the patch is ugly...

-- 
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/1d0a185d-84af-4be5-8d2c-9b7979528c1cn%40googlegroups.com.

[deal.II] Re: Mean value of component of TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

I am pretty sure that this is related to the following problem:

When I try to run AMR with this FESystem (Nedelec-RaviartThomas-DGQ) with 
more than one MPI rank I get this error

"An error occurred in line <1626> of file 
 
in function 
   void 
dealii::AffineConstraints::add_line(dealii::AffineConstraints::size_type)
 
[with number = double; dealii::AffineConstraints::size_type
= unsigned int] 
The violated condition was:  
   line_n != numbers::invalid_size_type 
Additional information:  
   This exception -- which is used in many places in the library -- usually 
indicates that some condition which the author of the code thought must be 
satisfied at a
certain point in an algorithm, is not fulfilled. An example would be that 
the first part of an algorithm sorts elements of an array in ascending 
order, and a second 
part of the algorithm later encounters an element that is not larger than 
the previous one. 

There is usually not very much you can do if you encounter such an 
exception since it indicates an error in deal.II, not in your own program. 
Try to come up with the 
smallest possible program that still demonstrates the error and contact the 
deal.II mailing lists with it to obtain help.
"

which is thrown in AffineConstraints::add_line 

 
(see the second assert in line 1826) after being called by 
make_hanging_node_constraints. It seems that deal.ii is trying to setup 
constraints on dofs that are not owned by the MPI process.

I read a little bit in the discussions referring to Nedelec and refinement 
- seems a difficult issue. Did anyone run into this, too?

Best,
Konrad

-- 
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/a3b92b45-1767-42ee-87ce-3e668a331a93n%40googlegroups.com.

[deal.II] Mean value of component of TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

Dear all,

I am trying to compute the mean value of the pressure component of a 
Trilinos block vector with three blocks and 7 components 
(vorticity-velocity-pressure).

Using one MPI rank is fine but if I use more I get the error

"An error occurred in line <666> of file 
 in 
function 
   dealii::TrilinosScalar 
dealii::TrilinosWrappers::MPI::Vector::operator()(dealii::TrilinosWrappers::MPI::Vector::size_type)
 
const 
The violated condition was:  
   false 
Additional information:  
   You tried to access element 3023 of a distributed vector, but this 
element is not stored on the current processor. Note: There are 4456 
elements stored on the current processor from within the range 5544 through 
 but Trilinos vectors need not store contiguous ranges on each 
processor, and not every element in this range may in fact be stored 
locally."

Now, I know what that means. Is there anything I can do? A workaround?

Thanks in advance and best regards,
Konrad

-- 
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/15e015b4-af0f-466e-b628-f3777ab64327n%40googlegroups.com.

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Dear Jean-Paul,

Many thanks for your reply.

On Tuesday, December 29, 2020 at 9:31:49 PM UTC+1 Jean-Paul Pelteret wrote:

> Hi Konrad,
>
> I'm sorry for taking some time to reply. To be honest, the inner working 
> of the FE classes is not something that I've ever had the time or 
> opportunity to dig into, so I'm really not well positioned to answer many 
> of your questions. I'm glad that you've submitted a PR for the work that 
> you've done thus far, and I hope that you can be guided by someone who has 
> more familiarity with this part of the library.
>
 
It will be my pleasure if other people will be able to benefit from 
implementations that I had to add and to contribute it as well as I 
benefited from other peoples‘ work. Let‘s hope that the stuff I asked to 
add makes sense. 
 

>
> What I can say, though, is that the various FE classes have been 
> implemented by numerous people over time, and the logic that dictates the 
> way that the DoFs are assigned varies between each FE. We have always 
> considered this to be perfectly fine -- this is, to the average person, an 
> implementational detail, and the interface that we have to interrogate DoFs 
> (e.g. which vector component they're associated with, whether or not they 
> have support points, etc...) seems to be sufficiently rich to most users. 
> So, as the end user, you shouldn't be writing code that would care that the 
> FE_Q element orders DoFs like this: [N_{0,x} N_{0,y} N_{0,z} N_{1,x} ...] ; 
> while IIRC the FE_DGQ element orders DoFs like this: [N_{0,x} N_{1,x} 
> N_{2,x} ... N_{0,y} N_{1,y} ...] .
>
> So, returning to your original set of questions: In general, I don't 
> believe that it can be assumed that all cell vertex DoFs are enumerated 
> before line/edges DoFs, before face DoFs, before cell interior DoFs. 
> Through the FiniteElementData class, the FiniteElement base class seems to 
> have a set of n_dofs_per_[vertex,line,face,cell]_dofs() 
> 
>  
> and get_first_[line/quad/hex]_index() 
> 
>  
> functions, which does suggest that I'm wrong about this, so I'd be happy to 
> be corrected on this point.  (Question: What does a vertex/edge/face DoF 
> even mean for a DG element?) Maybe these are the exact functions that you 
> were looking for (and maybe with this , if the [vertex,edge,face,cell] 
> groupings do apply, you can map from the cell DoF->associated geometric 
> entity using these functions). But I don't that they would help 
> understanding the ordering (e.g. lexiographic or something else; element 
> base index fast or spacedim fast, etc.) within each grouping -- again, 
> that's probably a decision that was left to the designer of each FE class, 
> and I don't think that there's been a need for any consistency between the 
> classes. But again, I'm just putting my thoughts out there and anticipate 
> being corrected on these points.
>

In my humble opinion the FE should only be implemented in terms of a 
reference element. This is, as far as I can tell, the case for all elements 
in Deal.II. The mapping though is problematic on complicated domains and I 
saw that one needs to renumber or switch signs for some shape functions if 
one wants to satisfy continuity conditions for the respective spaces. This 
is partly taken care of by the DoFHandler (see renumbering of FE_Q on 
faces) class which confused me, while other parts like the sign change is 
not. I guess I got the idea. Also, the concept of (generalized) support 
points is nice but sometimes it is not really clear to me if this is 
helpful in the context of element pairings (for example in the language of 
finite element exterior calculus where this concept does not really make 
sense globally). I may be wrong and miss a point here  and would love to 
get a lesson.
 

>
> I get the feeling that there is no overlap in assignment of DoFs between 
> geometric entities. Consider each of these entities not to have closure: so 
> you're actually referring to DoFs associated with the cell interior (excl. 
> faces, edges and vertices), face interior (excl edges and vertices) and 
> edge interiors (excl. vertices). To me, that the only sensible 
> interpretation of this. 
>
 
This is partly what is not fully clear to me. Vertex/edge/face/volume dofs 
(again FEEC) depends on the functionals that define them and less on 
location (a Raviart-Thomas or Nedelec shape function is not regular enough 
for point evaluations, although viewing them as a polynomial suggests so). 


> Thanks for giving this so much thought, and for the interesting questions. 
> I've certainly got a lot to learn from the responses that may come from the 
> other developers.
>
 
As I said, I hope other people can benefit from that but I am not super 
familiar with this part of the 

[deal.II] Re: Mean Value Constraints

[Konrad Simon]

Little correction: I wrote "In the vmult() function remove the mean value 
(i.e., project the rhs on the orthogonal complement of the kernel of the 
kernel)", I meant "In the vmult() function remove the mean value after you 
multiply (i.e., project the rhs on the orthogonal complement of the kernel 
of the system matrix)"

This paper <https://epubs.siam.org/doi/abs/10.1137/S0036144503426074> by 
Bochev and Lehoucq describes the technique with a toy model (Poisson with 
purely natural BCs). It was not new at that time but to the best of my 
knowledge they were the first ones to write this down in a proper way. The 
technique I mentioned above works for other kernels in system matrices as 
well. For example in linearized elasticity one maybe needs to remove not 
the mean value but rigid motions (well not really rigid but infinitesimal 
rotations - invariance under rotations cannot be linear) from the kernel. 
Such motions do not affect stress but constitute of nontrivial 
displacements. Also in this case one can project the rhs onto the 
orthogonal complement of the kernel of the system matrix. In this case the 
kernel has dimension 3 in 2D (one for scaling+90 degree rotation and two 
independent shifts) and dimension 6 in 3D. Once you project the rhs a 
Krylov solver can deal with your singular problem. 

Cheers,
Konrad 
On Saturday, December 26, 2020 at 11:06:36 AM UTC+1 Konrad Simon wrote:

> Hi,
>
> On Saturday, December 26, 2020 at 6:43:21 AM UTC+1 smetca...@gmail.com 
> wrote:
>
>> Hi all,
>>
>> Does anyone here have any experience applying mean value constraints 
>> (specifically with periodic boundary conditions)? I'm having some trouble. 
>> As far as I can tell, there are two approaches (both of which I have been 
>> able to get working). The first approach is to just add the mean value 
>> constraint straight into the constraint matrix. I got this working just 
>> fine and was able to get an output but the problem is that this condition 
>> seems to make my matrices dense(?) or, at the very least, require a huge 
>> amount of memory hindering this approach for the types of fine spatial 
>> meshes I require for this application. 
>>
>
> If you think of the mean value constraint as an additional row it looks 
> like (for Poisson equation) u_0+u_1+...+u_n=0 and hence it couples all 
> DoFs. This is mathematically ok but not efficient since this constraint is 
> non-local. 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. 
>  
>
>> The other approach is to just apply periodic boundary conditions to the 
>> matrices, try to solve the system anyway and, because UMFPACK is awesome, 
>> somehow manage to get a solution out of all of this (which will of course 
>> only be defined up to a constant). The problem here is that the solution I 
>> am getting is of order 10^12! When I postprocess and subtract the mean 
>> value, what I get looks correct but I am losing a lot of accuracy in the 
>> process but it does work for the fine meshes that I require. Anyone have 
>> any ideas or do I just have to live with the reduced accuracy here? Thanks!
>>
>
> This is of course like solving an ill-posed problem but some direct 
> solvers are smart enough to deal with zero pivots (redundant rows).
>
> Both approaches above do work but only for small and non-parallel 
> problems. This is what you can do:
>
> Many people do not think of the fact that even iterative solvers can deal 
> with singular problems if the right hand side is in the orthogonal 
> complement of the kernel. For the mean value constraint do the following:
>
> 1.  Setup boundary conditions in constraints if necessary. 
> 2. Assemble the system matrix and rhs and apply your constraints
> 3. Create a class that wraps your system matrix. This class must provide a 
> vmult() function
> 4. Before starting the iterative solver remove the mean value, see step-57 
> <https://www.dealii.org/current/doxygen/deal.II/step_56.html#StokesProblemprocess_solution>
>  
> 4. In the vmult() function remove the mean value (i.e., project the rhs on 
> the orthogonal complement of the kernel of the kernel)
> 5. use the wrapped matrix in the solver.
>
> Hope that helps.
>
>  Best,
> Konrad
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum

[deal.II] Re: Mean Value Constraints

[Konrad Simon]

Hi,

On Saturday, December 26, 2020 at 6:43:21 AM UTC+1 smetca...@gmail.com 
wrote:

> Hi all,
>
> Does anyone here have any experience applying mean value constraints 
> (specifically with periodic boundary conditions)? I'm having some trouble. 
> As far as I can tell, there are two approaches (both of which I have been 
> able to get working). The first approach is to just add the mean value 
> constraint straight into the constraint matrix. I got this working just 
> fine and was able to get an output but the problem is that this condition 
> seems to make my matrices dense(?) or, at the very least, require a huge 
> amount of memory hindering this approach for the types of fine spatial 
> meshes I require for this application. 
>

If you think of the mean value constraint as an additional row it looks 
like (for Poisson equation) u_0+u_1+...+u_n=0 and hence it couples all 
DoFs. This is mathematically ok but not efficient since this constraint is 
non-local. 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. 
 

> The other approach is to just apply periodic boundary conditions to the 
> matrices, try to solve the system anyway and, because UMFPACK is awesome, 
> somehow manage to get a solution out of all of this (which will of course 
> only be defined up to a constant). The problem here is that the solution I 
> am getting is of order 10^12! When I postprocess and subtract the mean 
> value, what I get looks correct but I am losing a lot of accuracy in the 
> process but it does work for the fine meshes that I require. Anyone have 
> any ideas or do I just have to live with the reduced accuracy here? Thanks!
>

This is of course like solving an ill-posed problem but some direct solvers 
are smart enough to deal with zero pivots (redundant rows).

Both approaches above do work but only for small and non-parallel problems. 
This is what you can do:

Many people do not think of the fact that even iterative solvers can deal 
with singular problems if the right hand side is in the orthogonal 
complement of the kernel. For the mean value constraint do the following:

1.  Setup boundary conditions in constraints if necessary. 
2. Assemble the system matrix and rhs and apply your constraints
3. Create a class that wraps your system matrix. This class must provide a 
vmult() function
4. Before starting the iterative solver remove the mean value, see step-57 

 
4. In the vmult() function remove the mean value (i.e., project the rhs on 
the orthogonal complement of the kernel of the kernel)
5. use the wrapped matrix in the solver.

Hope that helps.

 Best,
Konrad

-- 
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/67679343-a4c9-46c3-96b9-dd037445046fn%40googlegroups.com.

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Dear Deal.II community,

Suppose I chose FE_Nedelec<3> of some order as finite element.

1. Is there a way to query for a given cell_dof_index whether 
it is a face_dof and to get the face_index?

2. If this dof is a line_dof is it by definition also a 
face_dof? I need to exclude line_dofs and treat them 
separately.

3. Is there a convention that (if I am not using FESystem) in the set of 
all cell_dofs all line_dofs come before all 
face_dofs and only then the volume_dofs?

I am asking since I am trying to come up with a patch for RaviartThomas and 
a pattern to implement it in the current structures without interfering 
with them (too much).

Best,
Konrad
On Wednesday, December 16, 2020 at 9:33:57 PM UTC+1 Konrad Simon wrote:

> Dear Jean-Paul, dear Deal.II community,
>
> I partially solved the problem of sign flipping and permuting the degrees 
> of freedom and would like to come up with a merge-request at least for the 
> RaviartThomas elements. I have some questions before I can go on and 
> guidance would be appreciated.
>
> I decided to go a similar way that is taken for FE_Q elements. There one 
> has only a permutation of dofs on non-standard or flipped or rotated faces 
> since a sign change in the orientation does not matter. A vector of size 
> n_faces_per_cell of which each contains a table mapping each face_dof and 
> each of the 8 combinations of the flags (non-standard | flipped | rotated) 
> to an offset of that face_dof. The system behind that relies on the 
> lexicographic ordering of dofs if I understand correctly. My questions:
>
> 1. Is there a similar numbering scheme for dofs on faces for RT elements. 
> There one has face moments as dofs and it is not fully clear to me what 
> lexicographic means there (if such a principle is usable there at all).
>
> 2. If the order (and sign) of face_dofs if permuted for each of the flags 
> (non-standard | flipped | rotated) then in what order are the permutations 
> applied? (They do not commute)
>
> I would like to avoid hard-coding that (I think up to order 2 would be 
> doable -> but already 54 face_dofs). 
>
> I believe the table similar 
> to adjust_quad_dof_index_for_face_orientation_table of FE_Q_Base makes 
> sense but I guess every FE inheriting from FE_PolyTensor must implement 
> such a table then. This is a lot of work and more complicated for Nedelec 
> elements.
>
> Best,
> Konrad
> On Saturday, December 12, 2020 at 8:07:19 PM UTC+1 Konrad Simon wrote:
>
>> Hi Jean-Paul,
>>
>> On Thursday, December 10, 2020 at 11:39:08 PM UTC+1 Jean-Paul Pelteret 
>> wrote:
>>
>>> HI Konrad,
>>>
>>> I have no familiarity with the H-div elements, so I could be wrong with 
>>> this suggestion...
>>>
>>> The Fe_Nedelec element suffered from a similar issue, where adjacent 
>>> cells had to agree on the edge sense/directionality. This requirement could 
>>> not be unconditionally guaranteed for that element type, hence the 
>>> following note in its introduction 
>>> <https://dealii.org/current/doxygen/deal.II/classFE__Nedelec.html>
>>>
>>> Even if this element is implemented for two and three space dimensions, 
>>> the definition of the node values relies on consistently oriented faces in 
>>> 3D. Therefore, care should be taken on complicated meshes.
>>>
>>
>> Yes, this issue is similar. I assume I will need to check that for these 
>> elements as well since I also need them.
>>
>>
>>> The more recent FE_NedelecSZ 
>>> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> 
>>> element resolves the issue by querying the vertex numbers that make up each 
>>> edge/face and applying a rule that gives each face and its bounding edges a 
>>> direction. This means that the orientation of the face/edges is always 
>>> consistent no matter which of the two cells that share the face you're on. 
>>> The details are more explicitly laid out in the FE_NedelecSZ 
>>> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> 
>>> introduction, as well as the two papers that are listed therein. 
>>>
>>
>> This is a nice alternative for the Nedelec element. Often though such 
>> elements are used in conjunction with H1 or H(div) conformal elements and 
>> need to satisfy a compatibility condition (stronger than just LBB, i.e., 
>> they must form a bounded co-chain complex that is a subcomplex to something 
>> else). I am not sure that this element does satisfy this condition together 
>> with classical RT elements etc.
>>
>>
>>> Now, I'm by no means certain tha

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Dear Jean-Paul, dear Deal.II community,

I partially solved the problem of sign flipping and permuting the degrees 
of freedom and would like to come up with a merge-request at least for the 
RaviartThomas elements. I have some questions before I can go on and 
guidance would be appreciated.

I decided to go a similar way that is taken for FE_Q elements. There one 
has only a permutation of dofs on non-standard or flipped or rotated faces 
since a sign change in the orientation does not matter. A vector of size 
n_faces_per_cell of which each contains a table mapping each face_dof and 
each of the 8 combinations of the flags (non-standard | flipped | rotated) 
to an offset of that face_dof. The system behind that relies on the 
lexicographic ordering of dofs if I understand correctly. My questions:

1. Is there a similar numbering scheme for dofs on faces for RT elements. 
There one has face moments as dofs and it is not fully clear to me what 
lexicographic means there (if such a principle is usable there at all).

2. If the order (and sign) of face_dofs if permuted for each of the flags 
(non-standard | flipped | rotated) then in what order are the permutations 
applied? (They do not commute)

I would like to avoid hard-coding that (I think up to order 2 would be 
doable -> but already 54 face_dofs). 

I believe the table similar 
to adjust_quad_dof_index_for_face_orientation_table of FE_Q_Base makes 
sense but I guess every FE inheriting from FE_PolyTensor must implement 
such a table then. This is a lot of work and more complicated for Nedelec 
elements.

Best,
Konrad
On Saturday, December 12, 2020 at 8:07:19 PM UTC+1 Konrad Simon wrote:

> Hi Jean-Paul,
>
> On Thursday, December 10, 2020 at 11:39:08 PM UTC+1 Jean-Paul Pelteret 
> wrote:
>
>> HI Konrad,
>>
>> I have no familiarity with the H-div elements, so I could be wrong with 
>> this suggestion...
>>
>> The Fe_Nedelec element suffered from a similar issue, where adjacent 
>> cells had to agree on the edge sense/directionality. This requirement could 
>> not be unconditionally guaranteed for that element type, hence the 
>> following note in its introduction 
>> <https://dealii.org/current/doxygen/deal.II/classFE__Nedelec.html>
>>
>> Even if this element is implemented for two and three space dimensions, 
>> the definition of the node values relies on consistently oriented faces in 
>> 3D. Therefore, care should be taken on complicated meshes.
>>
>
> Yes, this issue is similar. I assume I will need to check that for these 
> elements as well since I also need them.
>
>
>> The more recent FE_NedelecSZ 
>> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> 
>> element resolves the issue by querying the vertex numbers that make up each 
>> edge/face and applying a rule that gives each face and its bounding edges a 
>> direction. This means that the orientation of the face/edges is always 
>> consistent no matter which of the two cells that share the face you're on. 
>> The details are more explicitly laid out in the FE_NedelecSZ 
>> <https://dealii.org/current/doxygen/deal.II/classFE__NedelecSZ.html> 
>> introduction, as well as the two papers that are listed therein. 
>>
>
> This is a nice alternative for the Nedelec element. Often though such 
> elements are used in conjunction with H1 or H(div) conformal elements and 
> need to satisfy a compatibility condition (stronger than just LBB, i.e., 
> they must form a bounded co-chain complex that is a subcomplex to something 
> else). I am not sure that this element does satisfy this condition together 
> with classical RT elements etc.
>
>
>> Now, I'm by no means certain that underlying issue that you're seeing 
>> here is the same as that which plagued the canonical Nedelec element 
>> (although your comment #3 makes me suspect that it could be). But if it is, 
>> then perhaps what was done to fix the orientation conflict in the Zaglmayr 
>> formulation of the Nedelec element can serve as inspiration for a fix to 
>> the BDM and RT elements?
>>
>> Thanks for looking into this, and I hope that you have some success at 
>> finding a solution to the problem. Naturally, if there's anything that we 
>> can do to help please let us know. It would be very nice to have a robust 
>> implementation to these elements in the library, so anything that you might 
>> be willing to contribute would be greatly appreciated!
>>
>  
> Lowest order RaviartThomas is already fixed. I will hunt down this issue 
> as good as I can. I might need some help with local and global dof 
> numbering conventions. I will also try to move the discussion to the issue 
> page <https://github.com/dea

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Hi Jean-Paul,

On Thursday, December 10, 2020 at 11:39:08 PM UTC+1 Jean-Paul Pelteret 
wrote:

> HI Konrad,
>
> I have no familiarity with the H-div elements, so I could be wrong with 
> this suggestion...
>
> The Fe_Nedelec element suffered from a similar issue, where adjacent cells 
> had to agree on the edge sense/directionality. This requirement could not 
> be unconditionally guaranteed for that element type, hence the following note 
> in its introduction 
> 
>
> Even if this element is implemented for two and three space dimensions, 
> the definition of the node values relies on consistently oriented faces in 
> 3D. Therefore, care should be taken on complicated meshes.
>

Yes, this issue is similar. I assume I will need to check that for these 
elements as well since I also need them.


> The more recent FE_NedelecSZ 
>  
> element resolves the issue by querying the vertex numbers that make up each 
> edge/face and applying a rule that gives each face and its bounding edges a 
> direction. This means that the orientation of the face/edges is always 
> consistent no matter which of the two cells that share the face you're on. 
> The details are more explicitly laid out in the FE_NedelecSZ 
>  
> introduction, as well as the two papers that are listed therein. 
>

This is a nice alternative for the Nedelec element. Often though such 
elements are used in conjunction with H1 or H(div) conformal elements and 
need to satisfy a compatibility condition (stronger than just LBB, i.e., 
they must form a bounded co-chain complex that is a subcomplex to something 
else). I am not sure that this element does satisfy this condition together 
with classical RT elements etc.


> Now, I'm by no means certain that underlying issue that you're seeing here 
> is the same as that which plagued the canonical Nedelec element (although 
> your comment #3 makes me suspect that it could be). But if it is, then 
> perhaps what was done to fix the orientation conflict in the Zaglmayr 
> formulation of the Nedelec element can serve as inspiration for a fix to 
> the BDM and RT elements?
>
> Thanks for looking into this, and I hope that you have some success at 
> finding a solution to the problem. Naturally, if there's anything that we 
> can do to help please let us know. It would be very nice to have a robust 
> implementation to these elements in the library, so anything that you might 
> be willing to contribute would be greatly appreciated!
>
 
Lowest order RaviartThomas is already fixed. I will hunt down this issue as 
good as I can. I might need some help with local and global dof numbering 
conventions. I will also try to move the discussion to the issue page 
 on github.

Best,
Konrad

-- 
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/90968d76-09ae-4dd5-82a8-c49bc6383792n%40googlegroups.com.

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Dear David, dear all,

It seems like the issue is much more involved than I initially thought. 
Some updates since I started going into this issue 
<https://github.com/dealii/dealii/issues/7970>.

1. I could solve the problem for lowest order RaviartThomas Elements (RT). 
This is already good.

2. The fix does not work for lowest order BDM elements.

3. The fix does not work for higher order elements, neither BDM nor RT 
elements
Suspected reason: For shape functions of higher order elements (whose 
normal components) are supported on a face flipping the sign on a flipped 
face is not enough since the flipped face changes/reverts the ordering of 
the face dofs hence the translation of local to global numbering of degrees 
of freedom changes for cells with flipped faces. I am not sure but at least 
I suspect that. Is this being taken care of by Deal.II somewhere?

I will go on chasing the bug and hopefully come up with a merge request 
soon.

Best,
Konrad

On Wednesday, December 9, 2020 at 3:44:58 PM UTC+1 Konrad Simon wrote:

> Hi David, 
>
> Many thanks for the hint. After some research I believe I stumbled over this 
> issue#7970 <https://github.com/dealii/dealii/issues/7970>? Let's see what 
> I can do.
>
> Best,
> Konrad
> On Tuesday, December 8, 2020 at 10:56:44 PM UTC+1 Wells, David wrote:
>
>> Hi Konrad,
>>
>> This isn't a very satisfying answer but, to the best of my knowledge, 
>> these Hdiv elements weren't written with general unstructured 3D meshes in 
>> mind. I suspect, without looking more closely, that this is a problem on 
>> our side.
>>
>> We document this for the FE_RaviartThomasNodal class but for 
>> FE_RaviartThomas this is just a comment in the source file. Would you be 
>> willing to write a patch explaining things better for the other Hdiv 
>> classes?
>>
>> Best,
>> David
>> --
>> *From:* dea...@googlegroups.com  on behalf of 
>> Konrad Simon 
>> *Sent:* Monday, December 7, 2020 12:48 PM
>> *To:* deal.II User Group 
>> *Subject:* [deal.II] Div-conforming elements in 3D 
>>  
>> Dear Deal.ii community, 
>>
>> Sorry for spamming the mailing list again. I posted some of the below 
>> stuff already in another post but that might have been with too many 
>> unnecessary details and confusing ideas from my side.
>>
>> I am not sure if I stumbled over something and would like to report on 
>> that. 
>>
>> Goal is to project a given function onto a finite element space over a 
>> certain mesh. The focus here shall be on H(div)-conforming elements such as 
>> Raviart-Thomas and Brezzi-Douglas-Marini Elements.
>>
>> I do so with several different meshes and I noticed an inconsistency in 
>> the projection on meshes with, say, not so simple topology (e.g., holes and 
>> voids). 
>>
>> Observations:
>> 1. The problematic meshes do have faces with non-standard orientation. My 
>> first thought was this was not carried over to the mapping of the 
>> H(div)-basis but I was wrong. The mesh orientation is consistent (as far as 
>> I can tell). A cube, for example, is not problematic. Please see two 
>> example plots below.
>>
>> 2. The issue arises with both Raviart-Thomas and Brezzi-Douglas-Marini 
>> Elements.
>>
>> 3. The issue arises in 3D only. 
>>
>> Any ideas what I ran into?
>>
>> Best,
>> Konrad[image: plate_with_hole.png][image: hyper_shell.png]
>>
>> -- 
>> 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/fa304b41-f83f-4c96-a2f7-c4e9e813553fn%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/fa304b41-f83f-4c96-a2f7-c4e9e813553fn%40googlegroups.com?utm_medium=email_source=footer>
>> .
>>
>

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

Re: [deal.II] Div-conforming elements in 3D

[Konrad Simon]

Hi David, 

Many thanks for the hint. After some research I believe I stumbled over this 
issue#7970 <https://github.com/dealii/dealii/issues/7970>? Let's see what I 
can do.

Best,
Konrad
On Tuesday, December 8, 2020 at 10:56:44 PM UTC+1 Wells, David wrote:

> Hi Konrad,
>
> This isn't a very satisfying answer but, to the best of my knowledge, 
> these Hdiv elements weren't written with general unstructured 3D meshes in 
> mind. I suspect, without looking more closely, that this is a problem on 
> our side.
>
> We document this for the FE_RaviartThomasNodal class but for 
> FE_RaviartThomas this is just a comment in the source file. Would you be 
> willing to write a patch explaining things better for the other Hdiv 
> classes?
>
> Best,
> David
> --
> *From:* dea...@googlegroups.com  on behalf of 
> Konrad Simon 
> *Sent:* Monday, December 7, 2020 12:48 PM
> *To:* deal.II User Group 
> *Subject:* [deal.II] Div-conforming elements in 3D 
>  
> Dear Deal.ii community, 
>
> Sorry for spamming the mailing list again. I posted some of the below 
> stuff already in another post but that might have been with too many 
> unnecessary details and confusing ideas from my side.
>
> I am not sure if I stumbled over something and would like to report on 
> that. 
>
> Goal is to project a given function onto a finite element space over a 
> certain mesh. The focus here shall be on H(div)-conforming elements such as 
> Raviart-Thomas and Brezzi-Douglas-Marini Elements.
>
> I do so with several different meshes and I noticed an inconsistency in 
> the projection on meshes with, say, not so simple topology (e.g., holes and 
> voids). 
>
> Observations:
> 1. The problematic meshes do have faces with non-standard orientation. My 
> first thought was this was not carried over to the mapping of the 
> H(div)-basis but I was wrong. The mesh orientation is consistent (as far as 
> I can tell). A cube, for example, is not problematic. Please see two 
> example plots below.
>
> 2. The issue arises with both Raviart-Thomas and Brezzi-Douglas-Marini 
> Elements.
>
> 3. The issue arises in 3D only. 
>
> Any ideas what I ran into?
>
> Best,
> Konrad[image: plate_with_hole.png][image: hyper_shell.png]
>
> -- 
> 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/fa304b41-f83f-4c96-a2f7-c4e9e813553fn%40googlegroups.com
>  
> <https://groups.google.com/d/msgid/dealii/fa304b41-f83f-4c96-a2f7-c4e9e813553fn%40googlegroups.com?utm_medium=email_source=footer>
> .
>

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

Re: [deal.II] Boundary conditions on components in with FESystem

[Konrad Simon]

Dear Wolfgang,

I may have stumbled over something. I created a 3D hyper_shell with 6 
coarse cells which is nothing but a cube with a cubical void in the center. 
I then asked for the face orientations of each face in each cell. In each 
of the 6 cells all faces have standard orientation according to the deal.ii 
convention. Mathematically this is not possible if you want a consistent 
orientation of normals cross the mesh (right?). I think, two cells must 
have flipped two flipped face normals for consistency.

In the GeometryInfo documentation it says:

However, it turns out that a significant number of 3d meshes cannot satisfy 
this convention. This is due to the fact that the face convention for one 
cell already implies something for the neighbor, since they share a common 
face and fixing it for the first cell also fixes the normal vectors of the 
opposite faces of both cells. It is easy to construct cases of loops of 
cells for which this leads to cases where we cannot find orientations for 
all faces that are consistent with this convention.

For this reason, above convention is only what we call the *standard 
orientation*. deal.II actually allows faces in 3d to have either the 
standard direction, or its opposite, in which case the lines that make up a 
cell would have reverted orders, and the above line equivalences would not 
hold any more. You can ask a cell whether a given face has standard 
orientation by calling cell->face_orientation(face_no): if the result 
is true, then the face has standard orientation, otherwise its normal 
vector is pointing the other direction. There are not very many places in 
application programs where you need this information actually, but a few 
places in the library make use of this. Note that in 2d, the result is 
always true. More information on the topic can be found in this glossary 

 article.
Is there a simple way to tell deal.ii to flip the normal of such a face? If 
yes, does this affect the edge orientation? I am asking since this is 
important for mapping Nedelec shape functions and I do not want to solve 
the problem for faces create one for edges.

A way out could be to simply flip the signs of (mapped) Raviart-Thomas 
shape functions that are supported on a non-consistent face. I guess this 
strategy would work for any H(div)-conformal shape function. Am I correct?

I am willing to fix that (in deal.ii itself) if necessary. 

Best,
Konrad

-- 
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/fc6bebbb-7f36-4e63-9028-f3ced9ba82f8n%40googlegroups.com.

[deal.II] Boundary conditions on components in with FESystem

[Konrad Simon]

Dear all,

I have a code version in which I am using a FESystem (3d) composed of three 
elements: FE_Nedelec (w), FE_RaviartThomas (u) and DGQ (p). 

The code compiles fine but gives unreasonable results and I do not see what 
I am doing wrong or what I forgot.

I need to impose (essential) boundary conditions on the first two 
components: 

   - on w: w x n =0 (n is normal)
   - on u: u*n=0 (no normal flux)

I assume that I am doing something wrong in the setup of BCs or the system 
(I guess the solver is likely ruled since two different solvers give 
exactly the same nonsense). This is the code snippet that I suspect to be 
buggy. I hope the names are self explanatory (essentially it is a 
modification of step-32). Note that dim=3 here.

Any help would be appreciated. 

Best,
Konrad

/
  // System and dof setup
  /

  template 
  void
  BoussinesqModel::setup_nse_matrices(
const std::vector _partitioning,
const std::vector _relevant_partitioning)
  {
nse_matrix.clear();
LA::BlockSparsityPattern sp(nse_partitioning,
nse_partitioning,
nse_relevant_partitioning,
this->mpi_communicator);

Table<2, DoFTools::Coupling> coupling(2 * dim + 1, 2 * dim + 1);
for (unsigned int c = 0; c < 2 * dim + 1; ++c)
  {
for (unsigned int d = 0; d < 2 * dim + 1; ++d)
  {
if (c < dim)
  {
if (d < 2 * dim - 1)
  coupling[c][d] = DoFTools::always;
else
  coupling[c][d] = DoFTools::none;
  }
else if ((c >= dim) && (c < 2 * dim))
  {
coupling[c][d] = DoFTools::always;
  }
else if (c == 2 * dim)
  {
if ((d >= dim) && (d < 2 * dim))
  coupling[c][d] = DoFTools::always;
else
  coupling[c][d] = DoFTools::none;
  }
  }
  }

DoFTools::make_sparsity_pattern(nse_dof_handler,
coupling,
sp,
nse_constraints,
false,
Utilities::MPI::this_mpi_process(
  this->mpi_communicator));
sp.compress();

nse_matrix.reinit(sp);
  }

  template 
  void
  BoussinesqModel::setup_dofs()
  {
TimerOutput::Scope timing_section(
  this->computing_timer, "BoussinesqModel - setup dofs of systems");

/*
 * Setup dof handlers for nse and temperature
 */
std::vector nse_sub_blocks(3, 0);
nse_sub_blocks[1] = 1;
nse_sub_blocks[2] = 2;

nse_dof_handler.distribute_dofs(nse_fe);

/*
 * Reduce bandwidth for ILU preconditioning
 */
DoFRenumbering::Cuthill_McKee(nse_dof_handler);

DoFRenumbering::block_wise(nse_dof_handler);

/*
 * Count dofs
 */
std::vector nse_dofs_per_block(3);
DoFTools::count_dofs_per_block(nse_dof_handler,
   nse_dofs_per_block,
   nse_sub_blocks);
const unsigned int n_w = nse_dofs_per_block[0], n_u = 
nse_dofs_per_block[1],
   n_p = nse_dofs_per_block[2],
   n_T = temperature_dof_handler.n_dofs();

/*
 * Comma separated large numbers
 */
std::locale s = this->pcout.get_stream().getloc();
this->pcout.get_stream().imbue(std::locale(""));

/*
 * Print some mesh and dof info
 */
this->pcout << "Number of active cells: "
<< this->triangulation.n_global_active_cells() << " (on "
<< this->triangulation.n_levels() << " levels)" << std::endl
<< "Number of degrees of freedom: " << n_w + n_u + n_p + n_T
<< " (" << n_w << " + " << n_u << " + " << n_p << " + " << 
n_T
<< ")" << std::endl
<< std::endl;
this->pcout.get_stream().imbue(s);

/*
 * Setup partitioners to store what dofs and matrix entries are stored 
on
 * the local processor
 */
{
  nse_index_set = nse_dof_handler.locally_owned_dofs();

  nse_partitioning.push_back(nse_index_set.get_view(0, n_w));
  nse_partitioning.push_back(nse_index_set.get_view(n_w, n_w + n_u));
  nse_partitioning.push_back(
nse_index_set.get_view(n_w + n_u, n_w + n_u + n_p));

  DoFTools::extract_locally_relevant_dofs(nse_dof_handler,
  nse_relevant_set);

  nse_relevant_partitioning.push_back(nse_relevant_set.get_view(0, 
n_w));
  nse_relevant_partitioning.push_back(
nse_relevant_set.get_view(n_w, n_w + n_u));
  nse_relevant_partitioning.push_back(

Re: [deal.II] Tools for parallel debugging

[Konrad Simon]

Thank you, Wolfgang and Daniel. Seems like I will go with the command line 
then. I was just wondering if people here use Eclipse's PTP which sounded 
like a good graphical tool. In my case it frequently crashes or simply gets 
stuck.

Best,
Konrad

On Sunday, October 11, 2020 at 11:35:34 PM UTC+2 Wolfgang Bangerth wrote:

> On 10/11/20 3:26 PM, Daniel Arndt wrote:
> > 
> > have a look at 
> > 
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-do-i-debug-mpi-programs
>  
> > <
> https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fwiki%2FFrequently-Asked-Questions%23how-do-i-debug-mpi-programs=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cab82fabba6ff45f8e49708d86e2c6359%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637380484218698974=jECVtixHfq%2BTSXkb1xww9SYWSyGwOonjWg%2Fqug7zLPk%3D=0>.
>  
>
> > For me, the "mpiexec -np n gdb ./executable" approach normally works 
> fairly well.
> > 
>
> There is also an interactive demonstration of this in lecture 41.25:
> https://www.math.colostate.edu/~bangerth/videos.676.41.25.html
>
> 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/3ed1b810-e78e-4cf0-9712-70aec571875en%40googlegroups.com.

Re: [deal.II] Tools for parallel debugging

[Konrad Simon]

Oh, seems like I missed that. Thanks you!

Konrad

On Sunday, October 11, 2020 at 11:35:34 PM UTC+2 Wolfgang Bangerth wrote:

> On 10/11/20 3:26 PM, Daniel Arndt wrote:
> > 
> > have a look at 
> > 
> https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-do-i-debug-mpi-programs
>  
> > <
> https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fdealii%2Fdealii%2Fwiki%2FFrequently-Asked-Questions%23how-do-i-debug-mpi-programs=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cab82fabba6ff45f8e49708d86e2c6359%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637380484218698974=jECVtixHfq%2BTSXkb1xww9SYWSyGwOonjWg%2Fqug7zLPk%3D=0>.
>  
>
> > For me, the "mpiexec -np n gdb ./executable" approach normally works 
> fairly well.
> > 
>
> There is also an interactive demonstration of this in lecture 41.25:
> https://www.math.colostate.edu/~bangerth/videos.676.41.25.html
>
> 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/3af7489d-9b56-4796-af4c-498291829c2an%40googlegroups.com.

[deal.II] Tools for parallel debugging

[Konrad Simon]

Hi deal.ii community,

I have a short question about tools (which probably also concerns many 
other people). 

I am using eclipse for C++ development and I frequently use the debugger. 
So far everything fine. Now I need to use parallel debugging tools for MPI 
but my eclipse crashes many times and/or is super slow.

What tools are you using? Any recommendations or hints?

Best,
Konrad

-- 
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/d89c4b39-a069-4825-b96b-741665e55095n%40googlegroups.com.

Re: [deal.II] Re: Evaluating FE-solution on distributed mesh, semi-Lagrangian method

[Konrad Simon]

Dear all,

I finally had time to go back to my code. Quick problem summary: I am 
currently working on a semi-Lagrangian method for an advection-diffusion 
equation with a multiscale method that I designed for advection dominated 
processes. During runtime I must evaluate the (known) multiscale 
FE-solution object at a previous time step. 

Problem: The mesh is distributed and the object that will give me the 
evaluated data at a point p is possibly owned on another processor. The 
processor that owes the object is the same as the one owing the cell that 
surrounds p. To find the correct cell on the same processor that needs the 
information I am doing things similarly as the FEFieldFunction class to 
find the cell. But what if p is in a cell that is not locally owned? It 
could be in a ghost cell or an artificial one (also ghosts no not owe the 
data I need).

Now, I have have to manually MPI communicate such points to find out the 
MPI rank that owes the cell surrounding p as I found out with your help. 
Any way of doing this efficiently? I would like to avoid the brute force 
version: communicating all points to all ranks and search the owned cells. 
(Hint: Mostly (but not always) p will be in a ghost layer. )

Any way to get the actual MPI rank of the owner?

Any help would be much appreciated.

Best,
Konrad

On Thursday, May 14, 2020 at 4:10:33 PM UTC+2 Konrad Simon wrote:

> Thank you, Bruno. :-)
>
>
> On Tuesday, May 12, 2020 at 2:50:35 PM UTC+2, Bruno Turcksin wrote:
>
>> Konrad, 
>>
>> There is nothing out of the box. However, deal.II uses p4est which can 
>> use more that one layer of ghost cells. So you should take a look 
>> there to see how hard it is to change the code. If that's the way you 
>> want to go, I am sure that someone will be able to provide you with 
>> some guidance. 
>>
>> Best, 
>>
>> Bruno 
>>
>> Le mar. 12 mai 2020 à 06:04, Konrad Simon  a écrit : 
>> > 
>> > Thank you, Bruno. What is sometimes done for semi-Lagrangian methods is 
>> that one defines a halo region around the locally owned cells (like ghost 
>> cells) that contains information from previous time steps and then limits 
>> the time step so that one never leaves the halo region when asking for 
>> information from previous time steps. Is there a way in Deal.ii to control 
>> the size of the halo region. In my case I have pretty uniform grids if it 
>> helps. 
>> > 
>> > Best, 
>> > Konrad 
>> > 
>> > On Monday, May 11, 2020 at 2:49:21 PM UTC+2, Bruno Turcksin wrote: 
>> >> 
>> >> Konrad, 
>> >> 
>> >> Unfortunately, you will need to do the communication yourself. You can 
>> only evaluate the solution on cells that locally owned or on ghost cells. 
>> >> 
>> >> 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 a topic in the 
>> Google Groups "deal.II User Group" group. 
>> > To unsubscribe from this topic, visit 
>> https://groups.google.com/d/topic/dealii/9li9twVxrGE/unsubscribe. 
>>
> > To unsubscribe from this group and all its topics, send an email to 
>> dea...@googlegroups.com. 
>>
> > To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/bfc6dd81-abb3-421c-8f6c-539c234b624f%40googlegroups.com.
>>  
>>
>>
>

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

Re: [deal.II] Re: Evaluating FE-solution on distributed mesh, semi-Lagrangian method

[Konrad Simon]

Thank you, Bruno. :-)

On Tuesday, May 12, 2020 at 2:50:35 PM UTC+2, Bruno Turcksin wrote:
>
> Konrad, 
>
> There is nothing out of the box. However, deal.II uses p4est which can 
> use more that one layer of ghost cells. So you should take a look 
> there to see how hard it is to change the code. If that's the way you 
> want to go, I am sure that someone will be able to provide you with 
> some guidance. 
>
> Best, 
>
> Bruno 
>
> Le mar. 12 mai 2020 à 06:04, Konrad Simon  > a écrit : 
> > 
> > Thank you, Bruno. What is sometimes done for semi-Lagrangian methods is 
> that one defines a halo region around the locally owned cells (like ghost 
> cells) that contains information from previous time steps and then limits 
> the time step so that one never leaves the halo region when asking for 
> information from previous time steps. Is there a way in Deal.ii to control 
> the size of the halo region. In my case I have pretty uniform grids if it 
> helps. 
> > 
> > Best, 
> > Konrad 
> > 
> > On Monday, May 11, 2020 at 2:49:21 PM UTC+2, Bruno Turcksin wrote: 
> >> 
> >> Konrad, 
> >> 
> >> Unfortunately, you will need to do the communication yourself. You can 
> only evaluate the solution on cells that locally owned or on ghost cells. 
> >> 
> >> 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 a topic in the 
> Google Groups "deal.II User Group" group. 
> > To unsubscribe from this topic, visit 
> https://groups.google.com/d/topic/dealii/9li9twVxrGE/unsubscribe. 
> > To unsubscribe from this group and all its topics, send an email to 
> dea...@googlegroups.com . 
> > To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/bfc6dd81-abb3-421c-8f6c-539c234b624f%40googlegroups.com.
>  
>
>

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

[deal.II] Evaluating FE-solution on distributed mesh, semi-Lagrangian method

[Konrad Simon]

Dear all,

I am currently working on a semi-Lagrangian method for an 
advection-diffusion equation. During runtime I must evaluate the (known) 
FE-solution at a previous time step but the mesh is distributed. 

Problem: It can well be that I must evaluate the solution at a point that 
is in a cell that is not locally owed by the processor requesting that 
value or not even a ghost cell. I am aware of the FEFieldFunction class and 
that it says that an exception will be thrown if the point is found to be 
in an artificial cell. 

Now, does anyone here have experience with that? Am I doomed and will I 
have to manually MPI communicate such points to any processor and test if 
they own the relevant cell?

Any help would be much appreciated.

Best,
Konrad

-- 
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/9e19c2b8-29e2-48fd-bc99-18b912696cfd%40googlegroups.com.

Re: [deal.II] Installation with spack - no access to adol-c

[Konrad Simon]

Jean-Paul, thanks a lot. Works! :-)

Cheers,
Konrad

-- 
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/47803471-9b2b-4af4-aa9c-b77d39a68af6%40googlegroups.com.

[deal.II] Installation with spack - no access to adol-c

[Konrad Simon]

Hi all,

A colleague and myself had problems to install dealii 9.1.1 via spack.

Quick question: Did anyone have problems to install deal.ii via spack 
lately? I see that it depends on the development version of adol-c but this 
one requires a gitlab account and access to adol-c. Does deal.ii v9.1.1 
really depend on the development version of adol-c? Other versions <= 2.6.3 
are on github.

Best,
Konrad

-- 
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/a09a369b-ef5c-48d4-bdf5-1c3c681d614c%40googlegroups.com.

[deal.II] Re: Interpolating to globally refined distributed mesh

[Konrad Simon]

Dear deal.ii community,

I solved it, sorry for bugging you with it but simple mistakes can bug you 
for long... 

I simply forgot to re-distribute the dofs for the finite element after 
refining the mesh. Ooofff :-/

Best,
Konrad

-- 
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/00bb0838-ea8a-4704-8a13-94089509a3c6%40googlegroups.com.

[deal.II] Re: Interpolating to globally refined distributed mesh

[Konrad Simon]

Add-on: The problem is in the function 
RefineInterpolate::refine_and_interpolate_on_distributed_mesh()

-- 
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/51923469-509d-4c73-ac27-cc5e5b79cce1%40googlegroups.com.

[deal.II] Re: Interpolating to globally refined distributed mesh

[Konrad Simon]

Add-on: The problem is in the function 
RefineInterpolate::refine_and_interpolate_on_distributed_mesh()

-- 
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/beca16b0-92ec-46f3-90bf-a1ccfc2913fb%40googlegroups.com.

[deal.II] Interpolating to globally refined distributed mesh

[Konrad Simon]

Dear deal.ii community,

I am trying to interpolate a distributed solution vector onto a globally 
refined distributed mesh. Actually I need to globally refine more than once 
but I believe I can iterate this procedure (is there a more efficient way? 
The finite element can be vector or scalar valued, and the solution vector 
a block vector)

However, despite looking at tutorials and the documentation I get an error 
(segmentation fault) without further explanation. I did not manage to track 
down the error - probably I am using the interface not correctly. Anyone 
has experience with this?

A minimal example is attached. (I used deal.ii 9.2.0 but it should work 
with 9.1.1, too.)

Best regards,
Konrad

-- 
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/9978ebc5-4ff9-411f-b892-d8220a2146ec%40googlegroups.com.
##
#  CMake script for the step-40 tutorial program:
##

# Set the name of the project and target:
SET(TARGET "interpolate")

# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
#  FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
#  FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
#  SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
  ${TARGET}.cc
  )

# Usually, you will not need to modify anything beyond this point...

CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)

FIND_PACKAGE(deal.II 9.2.0 QUIET
  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
  )
IF(NOT ${deal.II_FOUND})
  MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to 
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()

#
# Are all dependencies fulfilled?
#
IF(NOT ((DEAL_II_WITH_PETSC AND NOT DEAL_II_PETSC_WITH_COMPLEX) OR 
DEAL_II_WITH_TRILINOS) OR NOT DEAL_II_WITH_P4EST) # keep in one line
  MESSAGE(FATAL_ERROR "
Error! This tutorial requires a deal.II library that was configured with the 
following options:
DEAL_II_WITH_PETSC = ON
DEAL_II_PETSC_WITH_COMPLEX = OFF
DEAL_II_WITH_P4EST = ON
or
DEAL_II_WITH_TRILINOS = ON
DEAL_II_WITH_P4EST = ON
However, the deal.II library found at ${DEAL_II_PATH} was configured with these 
options
DEAL_II_WITH_PETSC = ${DEAL_II_WITH_PETSC}
DEAL_II_PETSC_WITH_COMPLEX = ${DEAL_II_PETSC_WITH_COMPLEX}
DEAL_II_WITH_P4EST = ${DEAL_II_WITH_P4EST}
DEAL_II_WITH_TRILINOS = ${DEAL_II_WITH_TRILINOS}
which conflict with the requirements.
One or both of the aforementioned combinations of prerequisites are not met by 
your installation, but at least one is required for this tutorial step."
)
ENDIF()

DEAL_II_INITIALIZE_CACHED_VARIABLES()
SET(CLEAN_UP_FILES *.log *.gmv *.gnuplot *.gpl *.eps *.pov *.vtk *.ucd *.d2 
*.vtu *.pvtu)
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
#include 
#include 

#include 

namespace LA {
using namespace dealii::LinearAlgebraTrilinos;
} // namespace LA

#include 
#include 
#include 

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

#include 
#include 
#include 
#include 

#include 
#include 
#include 

#include 
#include 

using namespace dealii;

template  class MyScalarFunction : public Function {
public:
  MyScalarFunction() : Function() {}

  virtual double value(const Point ,
   const unsigned int component = 0) const override;
  virtual void value_list(const std::vector> ,
  std::vector ,
  const unsigned int component = 0) const override;
};

template 
double MyScalarFunction::value(const Point ,
const unsigned int /*component*/) const {
  double return_value = 0;

  for (unsigned int d = 0; d < dim; ++d)
return_value += (p(d) - 0.5) * (p(d) - 0.5);

  return return_value;
}

template 
void MyScalarFunction::value_list(
const std::vector> , std::vector ,
const unsigned int /*component = 0*/) const {
  Assert(points.size() == values.size(),
 ExcDimensionMismatch(points.size(), values.size()));

  for (unsigned int p = 0; p < points.size(); ++p) {
values[p] = value(points[p]);
  } // 

Re: [deal.II] Re: cmake with library and executables

[Konrad Simon]

Many thanks, Matthias!

Works!

Best,
Konrad

On Wednesday, December 25, 2019 at 2:07:41 PM UTC+1, Matthias Maier wrote:
>
>
> On Fri, Dec 20, 2019, at 13:07 CST, Konrad Simon  > wrote: 
>
> > 
> ###
>  
>
> > 
> ###
>  
>
> > ADD_CUSTOM_TARGET(debug 
> >   COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR} 
> >   COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all 
> >   COMMENT "Switch CMAKE_BUILD_TYPE to Debug" 
> >   ) 
> > 
> > ADD_CUSTOM_TARGET(release 
> >   COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release 
> ${CMAKE_SOURCE_DIR} 
> >   COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all 
> >   COMMENT "Switch CMAKE_BUILD_TYPE to Release" 
> >   ) 
> > 
> ###
>  
>
> > 
> ###
>  
>
>
> Do you get the jobserver warning while running "make release" or "make 
> debug" with "-j8"? 
>
> If so - that's because you cannot recursively call into make with the 
> jobserver feature that way. 
>
> We have fixed the example steps a while ago but probably never updated 
> the documentation regarding these two custom targets. It is best to 
> simply have: 
>
> ADD_CUSTOM_TARGET(debug 
>   COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR} 
>   COMMENT "Switch CMAKE_BUILD_TYPE to Debug" 
>   ) 
>
> ADD_CUSTOM_TARGET(release 
>   COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR} 
>   ) 
>
> and then 
>
>   $ make release 
>   $ make -j8 
>
> or 
>
>   $ make debug 
>   $ make -j8 
>
> Best, 
> Matthias 
>

-- 
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/b39bba70-1033-4a18-b618-572309514c45%40googlegroups.com.

Re: [deal.II] Re: cmake with library and executables

[Konrad Simon]

Hi Wolfgang,

On Tuesday, December 24, 2019 at 5:59:31 PM UTC+1, Wolfgang Bangerth wrote:
>
>
> Konrad, 
> your email has no question :-) Is your problem that you can't call 'make 
> -j8' 
> and your question how to make that possible? If so, what happens if you 
> try? 
> What is the error message? 
>

The message essentially says I can not use the jobserver so I can only use 
once core to compile. My library does compile but falls back to make -j1 
every time I invoke make -j8. My question is why? And how can I use all 
cores?

Thank you and merry Christmas!

Konrad

-- 
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/1d1a9cfa-5975-432c-9093-b542ee213cab%40googlegroups.com.

[deal.II] Re: cmake with library and executables

[Konrad Simon]

Hello deal.ii community,

I am posting again (sorry) with a bit more info since I did not find the 
mistake in my cmake setup. I have essentially the following project 
structure:

project/source
project/include
project/doc
project/test. 

In project/ I have this CMakeLists.txt:

#
#
# Macro to print all arguments
macro (print_all_args)
  message (STATUS "")
  
  set(list_var "${ARGN}")
  foreach(_currentItem IN LISTS list_var)
message (STATUS "Adding library source file:   ${_currentItem}")
  endforeach (_currentItem)
  
  message (STATUS "")
endmacro (print_all_args)
#
#



###
###
message(STATUS "This is CMake ${CMAKE_VERSION}")
message(STATUS "")

cmake_minimum_required(VERSION 2.8.12)

find_package(deal.II 9.1.1 QUIET
  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
  )
if (NOT ${deal.II_FOUND})
  message(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. 
***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to 
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this 
path."
)
endif ()

DEAL_II_INITIALIZE_CACHED_VARIABLES()

set(PROJECT_NAME MsFEComplex)

# The version number
set(MsFEComplex_VER_MAJOR 0)
set(MsFEComplex_VER_MINOR 1)

project(${PROJECT_NAME})
###
###



###
###
# Check for the existence of various optional folders:
if (EXISTS ${CMAKE_SOURCE_DIR}/doc/CMakeLists.txt)
  set(MsFEComplex_HAVE_DOC_DIRECTORY TRUE)
endif ()

if (EXISTS ${CMAKE_SOURCE_DIR}/test/CMakeLists.txt)
  set (MsFEComplex_HAVE_TEST_DIRECTORY TRUE)
endif ()
###
###



###
###
# Change default CMAKE_INSTAL_PREFIX to ${CMAKE_BINARY_DIR}/lib
if (CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
  set(CMAKE_INSTALL_PREFIX "${CMAKE_SOURCE_DIR}/lib" CACHE PATH "default 
install path" FORCE )
endif ()

set (MsFEComplex_LIBRARY msfecomplex)
###
###



###
###
set(MsFEComplex_INCLUDE_DIR
${CMAKE_SOURCE_DIR}/include)

set(MsFEComplex_SRC_DIR 
${CMAKE_SOURCE_DIR}/source)
add_subdirectory (${MsFEComplex_SRC_DIR})


# Only add subdirectories for doc if exists
if (MsFEComplex_HAVE_DOC_DIRECTORY)
  add_subdirectory (${CMAKE_SOURCE_DIR}/doc)
endif ()

# Only add subdirectories for doc if exists
if (MsFEComplex_HAVE_TEST_DIRECTORY)
  add_subdirectory (${CMAKE_SOURCE_DIR}/test)
  ENABLE_TESTING()
endif ()
###
###



###
###
ADD_CUSTOM_TARGET(debug
  COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR}
  COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
  COMMENT "Switch CMAKE_BUILD_TYPE to Debug"
  )

ADD_CUSTOM_TARGET(release
  COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR}
  COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
  COMMENT "Switch CMAKE_BUILD_TYPE to Release"
  )
###
###


and in project/source I have:


###
###
#
# Include directory for sources
#
include_directories(${MsFEComplex_INCLUDE_DIR})
###

[deal.II] projecting function onto TrilinosWrappers::MPI::BlockVector

[Konrad Simon]

Hi all,

I am having a little problem with projecting a function onto (parts of) FE 
spaces. I am getting the error

The violated condition was:  
   (dynamic_cast *>( 
&(dof.get_triangulation())) == nullptr) 
Additional information:  
   You are trying to use functionality in deal.II that is currently not 
implemented. In many cases, this indicates that there simply didn't appear 
much of a need for 
it, or that the author of the original code did not have the time to 
implement a particular case. If you hit this exception, it is therefore 
worth the time to look int
o the code to find out whether you may be able to implement the missing 
functionality. If you do, please consider providing a patch to the deal.II 
development sources 
(see the deal.II website on how to contribute).

I of course get what this error message suggests and I am wondering if I 
could fix this somehow. The funny thing is that when I step through the 
code in debug mode I see that exactly the cast above fails. Funnily, the 
cast dynamic_cast *>( 
&(dof.get_triangulation())) works. 

Now I am asking myself why?  Am I missing something here?

Best regards,
Konrad


This is my function:

{
TrilinosWrappers::MPI::BlockVector locally_relevant_exact_solution;
locally_relevant_exact_solution.reinit(owned_partitioning,
mpi_communicator);

{ // write sigma to Nedelec_0 space

// Quadrature used for projection
QGauss<3> quad_rule (3);

// Setup function
ExactSolutionLin_A_curl exact_sigma(parameter_filename); // Exact solution 
for first FE ---> Nedelec FE

DoFHandler<3> dof_handler_fake (triangulation);
dof_handler_fake.distribute_dofs (fe.base_element(0));

if (parameters.renumber_dofs)
{
DoFRenumbering::Cuthill_McKee (dof_handler_fake);
}

AffineConstraints constraints_fake;
constraints_fake.clear ();
DoFTools::make_hanging_node_constraints (dof_handler_fake, 
constraints_fake);
constraints_fake.close();

VectorTools::project (dof_handler_fake,
constraints_fake,
quad_rule,
exact_sigma,
locally_relevant_exact_solution.block(0));

dof_handler_fake.clear ();
}

{// write u to Raviart-Thomas_0 space

// Quadrature used for projection
QGauss<3> quad_rule (3);

// Setup function
ExactSolutionLin exact_u(parameter_filename); // Exact solution for second 
FE ---> Raviart-Thomas FE

DoFHandler<3> dof_handler_fake (triangulation);
dof_handler_fake.distribute_dofs (fe.base_element(1)); 

if (parameters.renumber_dofs)
{
DoFRenumbering::Cuthill_McKee (dof_handler_fake);
}

AffineConstraints constraints_fake;
constraints_fake.clear ();
DoFTools::make_hanging_node_constraints (dof_handler_fake, 
constraints_fake);
constraints_fake.close();

VectorTools::project (dof_handler_fake,
constraints_fake,
quad_rule,
exact_u,
locally_relevant_exact_solution.block(1));

dof_handler_fake.clear ();
}
}

-- 
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/7a66b73f-4fe3-449f-b0ba-c4b16e8fbca3%40googlegroups.com.

[deal.II] Question about postprocessing

[Konrad Simon]

Hi deal.ii community,

I also have a little question about postprocessing in different spaces. I 
am post processing two solutions of the same problem but solved in two 
different (pairs of) spaces. One quantity, for example, is called u and is 
either in H(curl) or in H(div) depending on the form of the problem. 
Another is sigma and it is either in H^1 or in H(curl), respectively.

I need to compute things that involve the divergence of u when u is in 
H(div) or the curl when u is in H(curl). Both things I do using the 
DataPostprocessor class with entries of the (matrix valued) gradient of u 
but I don't know the internals.

My problem ist that when comparing quantities that should be similar 
according to the math then I get differences that are too large 
(intuitively). 

The thing is that when I have a quantity that is in H(curl) using Nedelec 
approximation then I can nicely take the curl but divergence will be zero 
by construction. Same if u is in H(div) with Raviart-Thomas approximation 
(then the curl is zero by construction). How is the computation of the 
gradient done internally?

Best,
Konrad

-- 
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/68c66960-82d8-4a63-92ac-9c104926959a%40googlegroups.com.

[deal.II] cmake with library and executables

[Konrad Simon]

Hi deal.ii community,

Little cmake question: I set up a user project with include, test, doc and 
source directory. I am collecting some code in a library and then I link a 
few application files to the library. I followed the guidelines of the 
documentation and everything is fine. Now one little but annoying thing. I 
get

make[4]: warning: jobserver unavailable: using -j1.  Add '+' to parent make 
rule.

when typing make -j8.

Do you have an idea what I could do?

Best,
Konrad

-- 
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/8d5a7e5e-3cc1-4fb9-bd3d-63f6db4622f8%40googlegroups.com.

Re: [deal.II] Evaluate shape functions for an element on a given cell in a triangulation

[Konrad Simon]

Many thanks, Wolfgang. 

-- 
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/8ce85110-2bf7-4f76-a4cb-ac1f9bdf3781%40googlegroups.com.

[deal.II] Re: How to use CellId class

[Konrad Simon]

Hi Zhidong,

On Monday, October 21, 2019 at 1:42:30 AM UTC+2, Zhidong Brian Zhang wrote:
>
> Thank you very much for your prompt reply, Konrad!
>
> My confusion is the output of cell->id(), for example,
> 0_3:000
> 0_3:200
> 0_3:003
> 0_3:006
> 0_3:406
> 0_3:606
> 0_3:206
> 0_3:007
> 0_3:407
> 0_3:607
> 0_3:207.
>
> What type I need to use as the key in std::map? And is the map 
> parallel-distributed? Because my model requires that.
>
> You can use the CellId object as a key itself: 

std::map cell_info_map;

If you need the map distributed across nodes according to the distribution 
of cells then you can initialize the map (on each node separately) through

for (; cell!=endc; ++cell)
{
if (cell->is_locally_owned())
{
CellId current_cell_id(cell->id());

std::pair::iterator, bool > 
result;
result = cell_basis_map.insert(std::make_pair(cell->id(),
cell_info_map));

 // ClassWithCellInfo must have a copy constructor
Assert(result.second,
ExcMessage ("Insertion of cell_info_map into std::map failed. "
"Problem with copy constructor?"));
 } // end if
} // end ++cell

Best,
Konrad

-- 
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/51256993-cd1b-487d-82a4-ee08ce56c501%40googlegroups.com.

[deal.II] Evaluate shape functions for an element on a given cell in a triangulation

[Konrad Simon]

Hi deal.ii community,

I have a little implementation problem to bug you with... 

Is there a way to evaluate a given shape function on a given cell like a 
normal Function or TensorFunction? I need to do many computations on a 
coarse cell that itself is meshed.

I have an implementation but it is not very efficient: I create an FEValues 
object and first map the evaluation points back to the unit cell using an 
appropriate mapping. These I use as "quadrature points" and get the values 
on the real (coarse) cell by then reinitializing the FEValues object with 
it. Works for all shape functions but is slow and cumbersome.

I actually only need it for Q1 elements and lowest order Nedelec and 
Raviart-Thomas elements.

Anyone has an idea?

Cheers,
Konrad

-- 
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/c3226e74-c327-4f63-8be9-ba71468a3a7e%40googlegroups.com.

[deal.II] Re: How to use CellId class

[Konrad Simon]

Hi Zhidong,

I don't know what exactly you want to do but I found it quite useful that 
the CellId class has an order relation. That makes it usable in a std::map. 
Use the CellId (you can retrieve it through cell-->id() if cell is a cell 
accessor like the one you use to loop over all your cells) as a key then 
and connect to it an object that, for example, contains specific 
information about your cell(s).

Best,
Konrad  

-- 
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/829003a9-f318-4790-b957-9684aec7b91e%40googlegroups.com.

[deal.II] Re: Issue with boost serialization and spack?

[Konrad Simon]

Hi Denis,

I don't have the build folders any more so I can not post the error log. 
But the error (using spack) occurred with both boost versions. I will post 
something once I will find a solution.
At any rate, thank you for your help.

Best,
Konrad 

-- 
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/4fab9bae-8ed2-49de-918c-a270a4070c33%40googlegroups.com.

[deal.II] Re: Issue with boost serialization and spack?

[Konrad Simon]

>
>
>> I guess I have a clue why I get the error. My backend nodes run on a 
>> different architecture. The openmpi version (i.e. mpirun) is the correct 
>> one. This is made sure in the slurm script.
>>
> That's always a pain to deal with. If you can, I would get an interactive 
> session on a compute node and compile everything on the compute node. You 
> can also compile everything in batch mode or compute on the login node and 
> pass the flag for the architecture of the compute node.
>
> Thank you, that's what I am trying. I will find a solution the question is 
jsut how much time it will eat...
Konrad

-- 
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/f3c230bb-68b1-4491-ae3b-04ec9f6c74f8%40googlegroups.com.

[deal.II] Re: Issue with boost serialization and spack?

[Konrad Simon]

>
> Thank you, Denis. I use a pretty stupid (but simple) workaround: I setup 
>> and compile deal.ii myself since all dependencies are compiled and use the 
>> cmake command used by spack. That works. And I do not get the serialization 
>> error.
>>
>
>
> now that is strange. Are you sure you pick up the same boost?
> Could you post the original error from CMake error logs that shows how 
> serialization test fails?
>  
>
Cmake in my little modification seems to pick up the right boost. (for 
error log with boost see above)
 
DEAL_II_WITH_BOOST set up with external dependencies 
#BOOST_VERSION = 1.70.0 
#BOOST_DIR = 
/scratch/cen/numgeo/spack-lib/linux-debian8-x86_64/gcc-7.3.0/boost-1.70.0 
#BOOST_CXX_FLAGS = -Wno-unused-local-typedefs 
#BOOST_DEFINITIONS = BOOST_NO_AUTO_PTR 
#BOOST_USER_DEFINITIONS = BOOST_NO_AUTO_PTR 
#BOOST_INCLUDE_DIRS = 
/scratch/cen/numgeo/spack-lib/linux-debian8-x86_64/gcc-7.3.0/boost-1.70.0/include
 

#BOOST_USER_INCLUDE_DIRS = 
/scratch/cen/numgeo/spack-lib/linux-debian8-x86_64/gcc-7.3.0/boost-1.70.0/include
 

#BOOST_LIBRARIES = 
/scratch/cen/numgeo/spack-lib/linux-debian8-x86_64/gcc-7.3.0/boost-1.70.0/lib/libboost_iostreams-mt.so;/scratch/cen/numgeo/spack-lib/lin
ux-debian8-x86_64/gcc-7.3.0/boost-1.


 

>
>> However, now my code runs on the machine I installed it on. But once I 
>> use slurm to distribute the job across nodes I get "illegal instruction" 
>> erros. Frustrating.
>>
>
> With HPC I used to have access to they did not use slurm, so I can't 
> really comment here.
> If that's during running quick tests or so, it could be related to a wrong 
> MPIEXEC pickedup by deal.II config.
> I wanted to fix that in Spack https://github.com/spack/spack/pull/11142 
> but apparently this solution may not be fully functional for Slurm.
>

I guess I have a clue why I get the error. My backend nodes run on a 
different architecture. The openmpi version (i.e. mpirun) is the correct 
one. This is made sure in the slurm script.
 
Best,
Konrad

-- 
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/d1e8e796-969e-4c41-ae0c-1c5a54742dd2%40googlegroups.com.

[deal.II] Re: Issue with boost serialization and spack?

[Konrad Simon]

Thank you, Denis. I use a pretty stupid (but simple) workaround: I setup 
and compile deal.ii myself since all dependencies are compiled and use the 
cmake command used by spack. That works. And I do not get the serialization 
error.

However, now my code runs on the machine I installed it on. But once I use 
slurm to distribute the job across nodes I get "illegal instruction" erros. 
Frustrating.

Nevertheless, thanks for your hint.

Best,
Konrad

-- 
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/49a66335-e58d-4eab-aee3-55e87a54f88f%40googlegroups.com.

[deal.II] Issue with boost serialization and spack?

[Konrad Simon]

Dear deal.ii community,

I am having a little problem and I was wondering if this issue is known. I 
installed deal.ii on our cluster  and all dependencies build nicely with my 
chosen compiler (gcc v7.3). BLAS and LAPACK are being built as well as 
openmpi (the versions on the cluster are not suitable at the moment). 
However, when building deal.ii I get the output below (I only paste the 
relevant part). It also happens when I ask preck to build deal.ii with 
(externally spack built) boost v1.62 and v1.68. 

Could there be a configuration issue with boost in the spack build class of 
deal.ii? When looking at it I see that some version dependent boost patches 
are applied but I am not familiar with the details.

Anyone knows this? Is there a workaround?

Thanks and best regards,
Konrad

-- Include 
/tmp/u290231/spack-stage/dealii-9.1.1-gfk33pgt5rojhujhmuruaahuuzyzq2zm/spack-src/cmake/configure/configure_2_boost.cmake
-- Found Boost: 
/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/include
 
(found suitable version "1.62.0", minimum required is "1.59") found 
components:  iostreams serialization system thread regex chrono date_time 
atomic 
--   BOOST_VERSION: 1.62.0
--   BOOST_LIBRARIES: 
/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_iostreams-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_serialization-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_system-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_thread-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_regex-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_chrono-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_date_time-mt.so;/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/lib/libboost_atomic-mt.so
--   BOOST_INCLUDE_DIRS: 
/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/include
--   BOOST_USER_INCLUDE_DIRS: 
/scratch/cen/numgeo/spack-lib/linux-debian8-broadwell/gcc-7.3.0/boost-1.62.0/include
-- Found BOOST
-- Performing Test BOOST_IOSTREAMS_USABLE
-- Performing Test BOOST_IOSTREAMS_USABLE - Success
-- Performing Test BOOST_SERIALIZATION_USABLE
-- Performing Test BOOST_SERIALIZATION_USABLE - Failed
-- The externally provided Boost.Serialization library failed to pass a 
crucial test. 
Therefore, the bundled boost package is used. 
The configured testing project can be found at 
/tmp/u290231/spack-stage/dealii-9.1.1-gfk33pgt5rojhujhmuruaahuuzyzq2zm/spack-build/cmake/configure/TestBoostBugWorkdir
-- DEAL_II_WITH_BOOST has unmet external dependencies.
CMake Error at cmake/macros/macro_configure_feature.cmake:112 (MESSAGE):
  

  Could not find the boost library!

  The externally provided Boost.Serialization library failed to pass a
  crucial test.Please ensure that a suitable boost library is installed on
  your computer.

  If the library is not at a default location, either provide some hints for
  autodetection,

  $ BOOST_DIR="..." cmake <...>
  $ cmake -DBOOST_DIR="..." <...>

  or set the relevant variables by hand in ccmake.

  Alternatively you may choose to compile the bundled library of boost by
  setting DEAL_II_ALLOW_BUNDLED=on or DEAL_II_FORCE_BUNDLED_BOOST=on.

Call Stack (most recent call first):
  cmake/macros/macro_configure_feature.cmake:269 (FEATURE_ERROR_MESSAGE)
  cmake/configure/configure_2_boost.cmake:237 (CONFIGURE_FEATURE)
  cmake/macros/macro_verbose_include.cmake:19 (INCLUDE)
  CMakeLists.txt:124 (VERBOSE_INCLUDE)

-- 
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/e889b864-66ad-4e05-81f3-08420d5e6a30%40googlegroups.com.

[deal.II] Re: Configuration of PETSc

[Konrad Simon]

Hi Toni,

Seems like I missed that little note in the documentation. Thank you :-)

Best,
Konrad

-- 
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/3148cde3-83f1-46a4-be88-3ddfe3e45cdc%40googlegroups.com.

[deal.II] Re: deal.II for Inverse problems in non-linear elasticity

[Konrad Simon]

Hi Prashant,
 

> I am trying to solve an Inverse Cauchy problem in 3D nonlinear elasticity. 
> I have observed displacement data at partial boundary as well in partial 
> regions inside the body, and want to reconstruct the traction field. From 
> documentation of deal.II, I understood how tangent stiffness matrix can be 
> created at each iteration for the forward problem.
>
 
This is a challenging ill-posed problem. I did something similar myself in 
the past. An advice from my experience: Start with linear elasticity 
(Navier-Lame). In fact the stiffness-matrix is the tangent stiffness matrix 
of Saint-Venant Materials around the zero-displacement field. Then with 
Saint-Venant. If this works try "simple" Neo-Hookean laws (keep in mind 
that tangent stiffness is not necessarily definite etc). etc... but I guess 
that was not your question. (Sorry for the unrequested advice)

While solving the inverse problem iteratively,  I want to retrieve this 
> tangent stiffness matrix at each iteration, perform necessary algebraic 
> manipulations and update the increments in unknown displacement and 
> traction variables.  Is it possible to do so using deal.II?  If you have 
> any suggestions, please feel free to include them.
>

Look at tutorial step-44. I think you will find what you are looking for.

Hope that helps.

Best,
Konrad

-- 
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/56f9a37a-58e0-4b56-a265-41ebdd1557b0%40googlegroups.com.

[deal.II] Configuration of PETSc

[Konrad Simon]

Dear deal.ii community,

I am using deali.ii with PETSc and Trilinos. However, when I am using the 
PETSc PreconditionILU I get an error that suggests that a solver package is 
missing (with Trilinos it works). Petsc's PreconditionAMG works fine 
(although not very efficiently for my problem). 
Do I need to do any special configuration steps for PETSC? I followed the 
instructions that are documented on the deal.ii pages on "how to configure 
Petsc". https://www.dealii.org/current/external-libs/petsc.html

Best,
Konrad

This is the error:

Running using PETSc. 
Number of active cells: 262144 
Total number of cells: 24865 (on 7 levels) 
Number of degrees of freedom: 1609920 (811200+798720) 
[0]PETSC ERROR: - Error Message 
-- 
[0]PETSC ERROR: See 
http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html for 
possible LU and Cholesky solvers 
[0]PETSC ERROR: Could not locate a solver package. Perhaps you must 
./configure with --download- 
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for 
trouble shooting. 
[0]PETSC ERROR: Petsc Release Version 3.9.4, Sep, 11, 2018  
[0]PETSC ERROR: main on a x86_64 named thunder5 by u290231 Fri Sep 20 
10:00:23 2019 
[0]PETSC ERROR: Configure options --with-shared-libraries=1 --with-x=0 
--with-mpi=yes --download-hypre=yes --with-64-bit-indices 
--with-debugging=yes --with-hypre=yes 
[0]PETSC ERROR: #1 MatGetFactor() line 4328 in 
/scratch/cen/numgeo/lib/petsc-3.9.4/src/mat/interface/matrix.c 
[0]PETSC ERROR: #2 PCSetUp_ILU() line 142 in 
/scratch/cen/numgeo/lib/petsc-3.9.4/src/ksp/pc/impls/factor/ilu/ilu.c 
[0]PETSC ERROR: #3 PCSetUp() line 923 in 
/scratch/cen/numgeo/lib/petsc-3.9.4/src/ksp/pc/interface/precon.c 
- 
TimerOutput objects finalize timed values printed to the 
screen by communicating over MPI in their destructors. 
Since an exception is currently uncaught, this 
synchronization (and subsequent output) will be skipped to 
avoid a possible deadlock. 
- 
WARNING! There are options you set that were not used! 
WARNING! could be spelling mistake, etc! 
Option left: name:-p value: ../MsFEComplex/parameter.in 
ERROR: Uncaught exception in MPI_InitFinalize on proc 0. Skipping 
MPI_Finalize() to avoid a deadlock. 


 
Exception on processing:  

 
An error occurred in line <421> of file 
 
in function 
   void dealii::PETScWrappers::PreconditionILU::initialize(const 
dealii::PETScWrappers::MatrixBase&, const 
dealii::PETScWrappers::PreconditionILU::AdditionalData&) 
The violated condition was:  
   ierr == 0 
Additional information:  
deal.II encountered an error while calling a PETSc function. 
The description of the error provided by PETSc is "See 
http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html for 
possible LU and Cholesky solvers". 
The numerical value of the original error code is 92. 
 

Aborting! 
 
-- 
Primary job  terminated normally, but 1 process returned 
a non-zero exit code. Per user-direction, the job has been aborted. 
-- 
-- 
mpirun detected that one or more processes exited with non-zero status, 
thus causing 
the job to be terminated. The first process to do so was: 

 Process name: [[37149,1],0] 
 Exit code:1 
--

-- 
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/71e73e1e-0a03-417d-b086-89b62e44a2e1%40googlegroups.com.

Re: [deal.II] Initialization of iterative solver

[Konrad Simon]

>
> Yes, you just set the solution vector to your best guess before you pass 
> it to the solver. But you'll find that in practice, this rarely reduced 
> the number of iterations by any significant amount :-( 
>

Many thanks, Wolfgang. Now, I also saw that in the source code of 
solver_cg.h. I really hoped that would speed up things a bit when solving 
time dependent problems.. Hm. 

Best,
Konrad

-- 
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/c1d3ff11-38c3-4c03-b8a0-5329b3247c27%40googlegroups.com.

[deal.II] Initialization of iterative solver

[Konrad Simon]

Dear deal.ii community,

I have a quick question: Suppose I know a vector that is close to the 
solution of my linear system. Is there a method to initialize an iterative 
solver with this vector?

Best,
Konrad

-- 
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/c9a7e762-bf1b-4ce3-91c1-515610ffdb5b%40googlegroups.com.