Re: [deal.II] Finite Volume Method - Matrix/RHS assembly

[Abbas Ballout]

Nick, 

2)
When you use DG0, all of the terms in the weak from cancel to zero except 
for the face term penallty*jump(u)*jump(v).  (du/dx and dv/dy are zero 
inside an element)
That penalty term in the DG community is seen as something that adds 
stability and isn't at the core of the scheme.  

In FV after doing your Gauss integral you end up asking yourself "what is 
the value of the gradient at the face" and set that to be mu*(u2-u1)/h. You 
get to construct one based on the values on an element and it's neighbour
That's a shortcoming of FV in a way because if you want to construct 
anything of higher order you'll have to querry values form the neighbours 
of neighbours of neighbours... In 2d/3d FV and with unstructured grids, 
constructing this isn't easy. 
The second term in the DG weak form: the jump(u) *mu * average(gradv), is 
supposed to be exactly that, but it gets cancelled out when you use DG0. In 
DG grad(v) on a face is calculated based on data of the local cell, in FV 
grad(v) gets calculated based on data from the neighbouring cells. BTW, 
that term is what you get from integrating by parts. 
Back to your question and if you insist on using DG0:
So to match the penalty*jump(u)*jump(v) with what you get from an FV  
scheme (mu*(u2-u1)/h) maybe you might try to set the penalty to be 
something like (mu/h) where h is the distance between the element centers. 
(That's not it tho)
 
1) 
That penalty term is there to enhance stability but I don't know if there 
is more to it. It penalizes the jump across the elements. See if going 
higher with a constant penalty helps with the accuracy.
Or maybe scale it with the element size. 
There are other DG tutorials for the Poisson equation maybe they'll have 
something there for you too. 
(Never liked that penalty term tbh either. All my homies hate the penalty 
parameter. JK). 

Abbas
On Sunday, June 23, 2024 at 1:24:38 AM UTC+2 Wolfgang Bangerth wrote:

>
> Nik:
>
> > (1) trying to solve in 1D:
> > 
> > * I get errors of the following type:
> > *error: ‘class dealii::DoFAccessor<0, 1, 1, false>’ has no member
> > named ‘measure’
> >   411 | re() / cell->face(face_no)->measure();* / I can make it run
> > by removing the calls to these functions and replace the penalty
> > factor which would involve them by a constant number, which according
> > to the output seems to work o.k. but might need further tuning of the
> > penalty constant (here just chosen as a constant of 100)./
>
> Yes, that's because the measure of a face in 1d is hard to define -- what 
> is 
> the measure of a point in a zero-dimensional space?
>
> If you have a good idea for how the 'extent1/2' variables should be 
> defined in 
> 1d, let us know and we can patch the program so that it will also work in 
> 1d.
>
>
> >  (2) trying to solve for polynomial degree 0 (i.e. piecewise constant):
> > 
> > * /When I try to solve the default convergence_rate test case in step-74
> > with degree p = 0, I don't see a decrease in the L2 norm of the solution.
> > I tested different values for the penalty_factor for which the given
> > function would evaluate to 0 because of the p*(p+1) term. I tested just
> > removing the p*(p+1) term and also tested constants like 1000 or 0.001,
> > but in all cases, I got results similar to the table below, whre the L2
> > norm does not decrease. The table below is for the case with penalty term
> > /= 0.5 *(1./cell_extent_left+1./cell_extent_right);
> > 
> > degree = 0
> > | cycle | cells | dofs  | L2| L2...red.rate.log2 | H1| 
> > H1...red.rate.log2 | Energy|
> > | 0 | 16| 16| 3.016e-01 | -  | 4.443e+00 | - 
> >| 6.045e+00 |
> > | 1 | 64| 64| 2.273e-01 | 0.41   | 4.443e+00 | 
> 0.00 
> > | 5.680e+00 |
> > | 2 | 256   | 256   | 2.284e-01 | -0.01  | 4.443e+00 | 
> 0.00 
> > | 5.554e+00 |
> > | 3 | 1024  | 1024  | 2.368e-01 | -0.05  | 4.443e+00 | 
> 0.00 
> > | 5.497e+00 |
> > | 4 | 4096  | 4096  | 2.428e-01 | -0.04  | 4.443e+00 | 
> 0.00 
> > | 5.470e+00 |
> > | 5 | 16384 | 16384 | 2.462e-01 | -0.02  | 4.443e+00 | 
> 0.00 
> > | 5.456e+00 |
> > | 6 | 65536 | 65536 | 2.481e-01 | -0.01  | 4.443e+00 | 
> 0.00 
> > | 5.448e+00 |
> > 
> > Maybe this is well-known and it shouldn't work. However, it would be 
> nice to 
> > be able to solve with degree 0 because that would allow me to mimic a 
> finite 
> > volume type implemenation.
>

[deal.II] Re: Finite Volume Method - Matrix/RHS assembly

[Abbas Ballout]

Nick, 
You can look at step 74. It has EXACTLY what you need. 


On Wednesday, June 19, 2024 at 10:33:18 PM UTC+2 nikl...@stanford.edu wrote:

> Dear all,
>
> I'm new to deal.ii and I would like to use it to gradually build a 
> simulator for charge/mass transport inside batteries.
>
> Since the equations I'm trying to solve are conservations laws, I would 
> like to use a finite volume type of implementation. I know that deal.ii is 
> a finite element library. However, since I would like to use a fully 
> implicit approach and later hopefully octree refined grids, I think deal.ii 
> offers the things needed for this task as well as many functionalities and 
> interfaces for future tasks.
>
> To start, I've summarized a simplified problem I would like to solve in 
> the appended .pdf file. It's a 1D Laplace equation using finite volumes. 
>
> I would like to build on the step-12 or step-12b tutorials because they 
> use discontinuous elements too. I don't call the beta function of the 
> original tutorials such that a 1D version should be possible.
>
> The appended nik-step-12.cc file contains a hard-coded version of the 
> three assembly functions that results in the correct matrix and solution 
> for a very specific case described in the .pdf (section 3).
>
> My questions would therefore be the following:
>
>- How should I best write the three functions (cell, boundary and 
>interior face worker) for the matrix/RHS assembly with as less hard-coding 
>as possible for a finite volume scheme like equations (7-9) in the .pdf.
>- Related - How can I access the elements of the solution vector 
>(phi_i) of the current cell and neighbor cell that are used for example in 
>equation (7) in the .pdf.
>
> Thank you very much in advance for your help or inputs.
>
> Best,
> Nik
>
>

-- 
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/356b0a48-8ccb-4c50-bb01-6ee6bcff8a22n%40googlegroups.com.

Re: [deal.II] Project cell quadrature point data to boundary

[Abbas Ballout]

Almonaci, 

I think you'll need the compute_projection_from_face_quadrature_points_matrix 
tfunction. 
Here's how you might use it:

std::vector distributed_field;
distributed_field.reinit(dof_handler.locally_owned_dofs(), 
mpi_communicator);

FullMatrix qpoint_to_dof_matrix(fe.dofs_per_face,
face_quadrature.size());

std::vector> local_values_at_qpoints, local_fe_values;
local_values_at_qpoints.reinit(face_quadrature.size());
local_fe_values.reinit(fe.dofs_per_cell);

for (const auto cell :dof_handler.active_cell_iterators())
for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++
face)
if(face_at_boundary)
{
FETools::compute_projection_from_face_quadrature_points_matrix(fe,
face_quadrature,
face_quadrature,
cell,
face,
qpoint_to_dof_matrix);

for (unsigned int q = 0; q < face_quadrature->size(); ++q)
{
local_values_at_qpoints[q] = kirchof_stress_or_something;
}

qpoint_to_dof_matrix.vmult_add(local_fe_values,
local_values_at_qpoints);


cell->set_dof_values(local_fe_values,
distributed_field);

local_fe_values = 0;
}


The  distributed_field should have what you want in the end. Hope this 
helps. 

Abbas 

-- 
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/e83497f4-4336-487d-aebc-41875a7d36dbn%40googlegroups.com.

Re: [deal.II] Shape gradients of (f(x)*\phi_i)

[Abbas Ballout]

Yes product rule (so broke can't even pay attention D:)

This is a "In order for me to explain this I have to test it, but in order 
to test this I have to explain it" situation
so to answer the why but extremely poorly: 
It is common in the spectral element community to use Gauss-Lobatto 
quadrature and shape functions defined 
at the same Guass-Lobatto points (no more). Since the integration is going 
to be inexact, aliasing errors are introduced 
and the descrte integral evaluations for the product rule are no longer 
equivalent (I guess?). 
The split form for advection is used to remove aliasing errors in such a 
scenario. 
It is also the defacto to write discrete derivatives as matrix operators 
there.   

You can see: Split form nodal discontinuous Galerkin schemes with 
summation-by-parts property for the compressible Euler equations
Or remark 3.7 in: Analysis of the SBP-SAT Stabilization for Finite Element 
Methods Part I: Linear Problems. 
  
I am probably missing something here so I'll be back if I figure this out. 
Asking if I can query for shape_grad(f*phi, point) was more of a shot in 
the dark. 
Thanks for the input as well
Abbas

On Tuesday, January 23, 2024 at 10:24:55 PM UTC+1 Wolfgang Bangerth wrote:

> On 1/23/24 13:09, Abbas Ballout wrote:
> > 
> > Sorry I meant to say I don't want to use that chain rule.
>
> You mean the product rule? Either way, why is it that you don't want to 
> use 
> it? The formula Daniel shows is an *identity*, not some approximation. 
> This is 
> how the derivative of a product is defined.
>
> 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/823febdb-9b58-4257-92de-7b195807376an%40googlegroups.com.

Re: [deal.II] Shape gradients of (f(x)*\phi_i)

[Abbas Ballout]

Daniel, 
Sorry I meant to say I don't want to use that chain rule.
It's because I am looking at papers that compare integrating both versions 
of the weak form:
post and prior using the chain rule.   
Abbas 

On Tuesday, January 23, 2024 at 7:47:14 PM UTC d.arnd...@gmail.com wrote:

> Abbas, 
>
> \dfrac{\partial f(x) \phi_i}{\partial x} is just \dfrac{\partial 
> f(x)}{\partial x} \phi_i+f(x)\dfrac{f(x)}{\partial x} which I would use.
> There is no need to integrate by parts.
>
> Best,
> Daniel
>
> On Tue, Jan 23, 2024 at 2:37 PM Abbas Ballout  
> wrote:
>
>> I know I can query the standard shape gradient at a quadrature point 
>> with  fe_values.shape_grad(i, q).
>> Is it possible to query for something like the gradient of the shape 
>> function multiplied by another function? I need \dfrac{\partial f(x) 
>> \phi_i}{\partial x} but I don't want to integrate by parts. 
>>
>> Abbas 
>>
>> -- 
>> 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/64e24fc5-9876-4e34-97c4-ba6bc77d22a4n%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/64e24fc5-9876-4e34-97c4-ba6bc77d22a4n%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/2d8de9d1-23f2-45a6-83ac-6311a779fa01n%40googlegroups.com.

[deal.II] Shape gradients of (f(x)*\phi_i)

[Abbas Ballout]

I know I can query the standard shape gradient at a quadrature point with  
fe_values.shape_grad(i, 
q).
Is it possible to query for something like the gradient of the shape 
function multiplied by another function? I need \dfrac{\partial f(x) 
\phi_i}{\partial x} but I don't want to integrate by parts. 

Abbas 

-- 
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/64e24fc5-9876-4e34-97c4-ba6bc77d22a4n%40googlegroups.com.

Re: [deal.II] Liltle help with MUMPS

[Abbas Ballout]

Thanks W. and Blais. 

I ran all six 2-core combinations and I got the same results. I passed 
--cpu-set 0,1,2,3 to np2 and I still got the same thing. 
I monitored the CPU usage and nothing is out of the ordinary.  

The good news is that I don't get this issue when I run on the desktop 
which has more cores and RAM.  
All good then c: 


On Wednesday, November 15, 2023 at 11:42:43 PM UTC+1 blais...@gmail.com 
wrote:

> In addition to what Wolfgang wrote,
> Sometimes the OS will map the cores of a CPU in the following order:
> Physical-Logical-Physical-Logical-Physical-Logical etc.
> Consequently, if your processor supports hyperthreading through the use of 
> logical core, running with core 0-1 means you are essentially running on a 
> single physical core but using two logical cores. For intensive operations 
> that use the same instruction set, this is generally never a good idea. In 
> the present case though, since the assembly is fast but not the matrix 
> solve, I would presume it is more a question of memory lane issues than 
> anything else. Your CPU might be really memory-bandwith bound.
>
>
> On Wednesday, November 15, 2023 at 11:39:59 a.m. UTC-5 Wolfgang Bangerth 
> wrote:
>
>>
>>
>> On 11/15/23 04:09, Abbas Ballout wrote: 
>> > 
>> > If I run mpirun -np2 the assemble times are 25.8 seconds and the MUMPS 
>> > solve times are 51.7 seconds. 
>> > If I run mpirun --cpu-set 0-1 -np 2, the assemble times are 26 seconds 
>> > (unchanged) but the solve time are at 94.9 seconds! 
>> > Is this normal and expected? 
>>
>> What happens if you use `--cpu-set 0,2` or other combinations? 
>>
>> Cores have memory channels that they often share with other cores. I 
>> wonder whether, for example, cores 0 and 1 share a memory channel and 
>> consequently step on each other's toes during the direct solve (which is 
>> memory bound). Or they share SIMD execution units. It would be 
>> interesting to see what happens if you choose other sets of cores to use. 
>>
>> Best 
>> W. 
>>
>

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

Re: [deal.II] In example 40, the results are significantly different between the first-order and second-order finite element meshes

[Abbas Ballout]

ztdep, 
You got me curious 

You can output the results of the Kelly error estimate with something like 
this in the data output: 








*Vector 
estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator::estimate(dof_handler,QGauss(fe.degree + 1),std::map 
*>(),locally_relevant_solution,estimated_error_per_cell);  
  data_out.add_data_vector(estimated_error_per_cell, "kelly");*

and you'll see that the estimate is horrible with fe(1). Why? because the 
estimate is based on the jump of the gradients across cells and fe(1) isn't 
cutting it 
So if you output the gradient with:


*GradientPostprocessor gradient_postprocessor;
data_out.add_data_vector(locally_relevant_solution, 
gradient_postprocessor);*
(copy paste the postprocessor 
from 
https://www.dealii.org/current/doxygen/deal.II/classDataPostprocessorVector.html)
 

You can see in Paraview that the gradient doesn't look continuous/clean 
across elements when you use fe(1). 
So this isn't weird and unexpected. 

Apparently the gradient approx with bi-linear fe(1) is bad? 
The only way would be for me to be 100% sure would be to write these terms 
down and check for myself. 
If you do it let me know. 


On Monday, November 13, 2023 at 1:45:17 PM UTC+1 blais...@gmail.com wrote:

> Changing to Fe1 really alters the solution, so consequently the kelly 
> error estimator will also give you a different error estimation and thus a 
> different mesh adaptation.
> If your mesh adaptation is based on an error estimator, altering the 
> finite element interpolation order will alter the solution and thus alter 
> the mesh adaptation process.
>
>
> On Monday, November 6, 2023 at 9:13:31 p.m. UTC-5 ztdep...@gmail.com 
> wrote:
>
>>
>> The default setting is  'fe(2)" , we can clearly see a refined region in 
>>  the resulting mesh. While after i changed it to "fe(1)", we cann't see 
>> this region.  
>> This is what i wnat to know .
>> ztdepyahoo
>> ztdep...@gmail.com
>>
>> 
>>  
>>  Replied Message  
>> From Wolfgang Bangerth 
>> Date 11/7/2023 02:59 
>> To  
>> Subject Re: [deal.II] In example 40, the results are significantly 
>> different between the first-order and second-order finite element meshes 
>> On 11/6/23 03:12, ztdep...@gmail.com wrote:
>>
>> **
>>
>> the adaptivited mesh are very diffent as shown in attachments.
>>
>>
>> @ztdep:
>> We don't actually know what you are asking: Your post has no question. We 
>> also 
>> don't know what it is you did: The post doesn't show how you modified 
>> step-40. 
>> Finally, you do not say why you think that the result is wrong or 
>> surprising: 
>> You just show a picture, but for all we know this may be correct.
>>
>> Please spend a bit of thought in formulating questions in such a way that 
>> it 
>> is clear what you are asking, why you are asking it, and what is behind 
>> the 
>> question you are asking.
>>
>> 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 a topic in the 
>> Google Groups "deal.II User Group" group.
>> To unsubscribe from this topic, visit 
>> https://groups.google.com/d/topic/dealii/VcsXMXhl_YQ/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to 
>> dealii+un...@googlegroups.com.
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/e17928b7-f1d0-e5f2-83b9-872ba95e492a%40colostate.edu
>> .
>>
>

-- 
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/0a611a56-c142-4153-8220-f42c394bcb00n%40googlegroups.com.

[deal.II] Liltle help with MUMPS

[Abbas Ballout]

This isn't a dealii problem per-se. 

I am trying to run a number of simulations with different parameters of the 
same code with the underlying solver being MUMPS. I am using  *mpirun 
--cpu-set * to bind and isolate the different simulations to different 
cores (As I believe I should) 

*Profiling: *
The number of dofs are something like 262,144, the assemble function is 
called 25 times and the solve function is called 20 times. The solve 
function is using MUMPS nothing fancy is happening there. 
If I run mpirun -np2 the assemble times are 25.8 seconds and the MUMPS 
solve times are 51.7 seconds. 
If I run mpirun --cpu-set 0-1 -np 2, the assemble times are 26 seconds 
(unchanged) but the solve time are at 94.9 seconds! 
 
Is this normal and expected? 

*Maybe relevant details: *
I have 4 physical cores. I am on ubuntu. PETScVectors and  PETScWrappers::
SparseDirectMUMPS 


-- 
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/24f65840-62f0-4d73-b9e0-06a2770ad7f0n%40googlegroups.com.

[deal.II] Re: Difference between 2 Wasy to Initialize System Matrix in MPI World

[Abbas Ballout]

Did you try to run in Debug? If not give it a try and see if an exception 
is thrown. 

Abbas

On Thursday, September 28, 2023 at 3:15:14 AM UTC+2 hitl...@gmail.com wrote:

> Hello all, 
>
>
> Previously, I initialized my system matrix in this way (option A)
>
> DoFTools::make_flux_sparsity_pattern(dof_handler,
> dsp,
> cell_coupling,
> face_coupling);
> constraints_newton_update.condense(dsp);
> SparsityTools::distribute_sparsity_pattern(dsp,
>hp_index_set,
>mpi_communicator,
>hp_relevant_set);
> system_matrix.reinit(hp_index_set,
>  hp_index_set,
>  dsp,
>  mpi_communicator);
>
> However, the generated system matrix is singular. I was stuck by the 
> singular matrix problem for about 2 weeks.
>
>
> Then, I initialized the matrix with option B, which generates a full rank 
> matrix. Although, I don't see the difference between these 2 ways.  Can 
> someone explain this to me? Thans in advance.
>
> DoFTools::make_flux_sparsity_pattern(dof_handler,
>  dsp,
>  cell_coupling,
>  face_coupling);
> constraints_newton_update.condense(dsp);
> sparsity_pattern.copy_from(dsp);
> system_matrix.reinit(hp_index_set,
>  sparsity_pattern,
>  mpi_communicator);
>

-- 
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/7bcd69b4-708a-44cf-9d91-32e840ae4ffdn%40googlegroups.com.

Re: [deal.II] Neumann Boundary condition in MeshWorker for Automatic Differentiation

[Abbas Ballout]

Jost, 

For looping on the boundaries, you'll have to write an extra lambda amd 
pass that to the MeshWorker. 
You can look at step 12 or 74 for this.
Also following this cause I want to know what are the MeshWorker 
alternatives since I am relying heavily on it in my code

Best, 
Abbas
 
On Thursday, August 3, 2023 at 3:11:48 PM UTC+2 jost@googlemail.com 
wrote:

> Thanks for the quick response!
>
> Would you suggest me to read more tutorials on how the MeshWorker works to 
> fix this issue then or try to set up my program without a MeshWorker? I 
> thought I had set it up correctly and tried already for a while to figure 
> out what I did wrong. Also I Initially only started using the MeshWorker 
> because it seemed to be recommended in Step-72 and I thought following the 
> recommendations is probably the best and easiest way. Now I am not so sure 
> anymore about it.
>
> Best,
>
> Jost
>
> Wolfgang Bangerth schrieb am Donnerstag, 3. August 2023 um 14:49:44 UTC+2:
>
>> On 8/3/23 04:54, 'Jost Arndt' via deal.II User Group wrote: 
>> > 
>> > for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) 
>> > { 
>> > if 
>> > (cell->face(f)->at_boundary())//(face->at_boundary()) 
>> > { 
>> > std::cout << "reached line 
>> 1032" << 
>> > std::endl; 
>> > const auto _face_values = 
>> > scratch_data.reinit(cell, f); 
>> > std::cout << "reached line 
>> 1034" << 
>> > std::endl; 
>> >  /* here integration using 
>> > fe_face_value should happen*/ 
>> > 
>> > but the last cout never gets executed, but the code crashes. I get a 
>> long but 
>> > cryptic (to me?) stacktrace but also the following error message: 
>> > 
>> > 
>> > An error occurred in line <3044> of file <./source/fe/fe_values.cc> in 
>> function 
>> > dealii::FEValuesBase::FEValuesBase(unsigned int, 
>> unsigned int, 
>> > dealii::UpdateFlags, const dealii::Mapping&, const 
>> > dealii::FiniteElement&) [with int dim = 2; int spacedim 
>> >  = 2] 
>> > The violated condition was: 
>> > n_q_points > 0 
>> > Additional information: 
>> > There is nothing useful you can do with an FEValues object when 
>> using 
>> > a quadrature formula with zero quadrature points! 
>> > 
>> > From here I am quite unsure what I did wrong or how else to create the 
>> > FEFaceValues Object? I would be really thankful about any hint! 
>>
>> The error likely means that MeshWorker hasn't initialized the 
>> FEFaceValues 
>> object. Nobody knows any more how MeshWorker really works, and so we try 
>> to 
>> discourage use of that framework. But I do know that MeshWorker allows 
>> you to 
>> not only loop over the cells, but also the faces of a triangulation. That 
>> is 
>> the recommended way to deal with face integrals (rather than dealing with 
>> these face integrals as part of cell integration -- because you want to 
>> make 
>> sure that you visit each face only once, rather than twice in your 
>> approach) 
>> and if you tell MeshWorker that you also want to use a face worker, then 
>> it 
>> should (?) initialize the FEFaceValues object as well. 
>>
>> 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/90721349-725a-49e7-976c-7d989be0812bn%40googlegroups.com.

[deal.II] Re: How to use cell stiffness matrix to build global stiffness matrix

[Abbas Ballout]

Not really experienced with eigen so I am not exactly sure. 
Abbas

On Tuesday, August 1, 2023 at 6:53:35 PM UTC+2 dim...@gmail.com wrote:

> Hello Abbas,
>
> Thanks for your reply.
>
> I will take a look at the step13.
>
> Howerver, I would like to know if the cell matrix transform seems like the 
> process below:
>
> const Eigen::MatrixXd& cell_matrix,
> const Eigen::VectorXd& cell_youngs_modulus,
> const std::vector& local_dof_indices,
> Eigen::MatrixXd& tangent_matrix,
> Eigen::VectorXd& global_youngs_modulus)
>
>const int N = cell_matrix.rows();
> const int M = tangent_matrix.rows();
>
> // Loop over the local degrees of freedom
> for (int i = 0; i < N; ++i) {
> const int global_i = local_dof_indices[i]; // Global index of the 
> i-th local dof
>
> // Add the contribution of the i-th local dof to the global 
> Young's modulus vector
> global_youngs_modulus(global_i) += cell_youngs_modulus(i);
>
> for (int j = 0; j < N; ++j) {
> const int global_j = local_dof_indices[j]; // Global index of 
> the j-th local dof
>
> // Add the contribution of the (i,j)-th entry of the cell 
> matrix multiplied by the i-th and j-th Young's modulus to the global 
> tangent matrix
> tangent_matrix(global_i, global_j) += cell_matrix(i, j) * 
> cell_youngs_modulus(i) * cell_youngs_modulus(j);
> }
> }
>
> I am not quite sure if the assembly process for the global stiffness 
> matrix is correct.
>
> Best regards
> Lance
>
> On Tuesday, August 1, 2023 at 11:47:12 AM UTC+2 abbas.b...@gmail.com 
> wrote:
>
>> Lance, 
>>
>> I am not sure if I understood you correctly. 
>> Maybe looking at the local_assemble_matrix() function in step:13 might 
>> help. 
>>
>> Abbas
>>
>> On Tuesday, August 1, 2023 at 11:00:10 AM UTC+2 dim...@gmail.com wrote:
>>
>>> Hello dear group,
>>>
>>> I have one question which is about how to build global stiffness matrix 
>>> with cell stiffness matrix.
>>>
>>> In our project,the density is a vector (actually we need a mu vector but 
>>> in the code of deallii mu is a scalar)with different values,I would like to 
>>> use the element of density vector to build cell matrix as written in the 
>>> link below:
>>>
>>>  
>>> https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
>>> (starting from line 1782)
>>>
>>> And finally the cell stiffness matrix is used to build global stiffness 
>>> matrix,but I don' know how the global stiffness matrix is bulit with cell 
>>> stiffness matrix.
>>>
>>>
>>> https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
>>> (starting from line 1599)
>>> //>>>>>>>>>>>>>>>>>>>>>>>>>>
>>>  constraints.distribute_local_to_global(
>>> data.cell_matrix, 
>>> data.cell_rhs,
>>> data.local_dof_indices,
>>> tangent_matrix, 
>>> system_rhs);
>>> }
>>>
>>>
>>> //<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>>
>>>
>>> I have a idea that I use mu[0] to obtain the first cell_matrix and use 
>>> mu[1] to get the second cell stiffness matrix ,step by step then I got 
>>> final mu[m] to obtain the last cell stiffness matrix and in the end I use 
>>> some function like  distribute_local_to_global to construct the global 
>>> stiffness matrix,and system rhs and even the solutions.
>>>
>>>I checked the code in 
>>> AffineConstraints::distribute_local_to_global 
>>> <https://www.dealii.org/current/doxygen/deal.II/classAffineConstraints.html#a373fbdacd8c486e675b8d2bff8943192>.
>>>  
>>> (
>>> https://www.dealii.org/current/doxygen/deal.II/affine__constraints_8h_source.html
>>> )
>>> But the hint information guides me to read the cm.templates.h file.
>>> ewcfp I was a little confused ,I did not find the file from dealii lib. 
>>> could anyone provide any information or hint? Thanks in advance!
>>> Best regards
>>> Lance
>>> Could 
>>>
>>> I
>>>
>>

-- 
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/33214d24-1337-4a01-84a0-932e8bbd6b77n%40googlegroups.com.

Re: [deal.II] Q: Question about extracting part of a vector

[Abbas Ballout]

Najwa, 
The files I edited are in the pull request here: 
https://github.com/dealii/dealii/pull/15805.
Let me know if I can help in a better way or of there is something I can 
do. 

Abbas. 

On Tuesday, August 1, 2023 at 5:48:41 AM UTC+2 najwaa...@gmail.com wrote:

> Hello Abbas,
>
> Thank you for your response. I have actually written something, but I need 
> to test it. It would be great to have a look at your work as well, so we 
> can compare.
>
> Best regards,
> Najwa
>
> On Mon, 31 Jul 2023 at 10:48 PM Abbas Ballout  
> wrote:
>
>> Najwaa, 
>>
>> I submitted a pull request <https://github.com/dealii/dealii/pull/15805> 
>> recently about something similar 
>> <https://groups.google.com/g/dealii/c/H0kLG8RJKVc>. (I guess)
>> Maybe that would help. 
>>
>> Abbas 
>>  
>>
>> On Monday, July 31, 2023 at 7:37:59 PM UTC+2 najwaa...@gmail.com wrote:
>>
>>> Dear Wolfgang,
>>>
>>> Thank you for your answer. I will try to write a patch to the deal.II 
>>> sources that implement the missing function. This might require some time. 
>>> I will come back here if needed.
>>>
>>>
>>> Thank you all for your quick answers,
>>> Najwa
>>>
>>> On Wednesday, July 26, 2023 at 11:41:15 PM UTC+3 Wolfgang Bangerth wrote:
>>>
>>>> On 7/26/23 12:17, Najwa Alshehri wrote: 
>>>> > 
>>>> > Regarding the function "VectorFucntionFromScalarFunctionObject," I 
>>>> have 
>>>> > observed that it works for computing the L2 norm of the error. 
>>>> However, when I 
>>>> > attempted to compute the H1_seminorm of the error, I noticed that the 
>>>> gradient 
>>>> > was not implemented for this function. Please let me know if my 
>>>> understanding 
>>>> > is incorrect, as I referred to the  for 
>>>> this. 
>>>>
>>>> Yes, correct. Someone needs to implement 
>>>> VectorFucntionFromScalarFunctionObject::gradient() 
>>>> and/or 
>>>> VectorFucntionFromScalarFunctionObject::gradient_list() 
>>>> in exactly the same way as the value() and value_list() functions are 
>>>> implemented. 
>>>>
>>>>
>>>> > Could you kindly advise me on the simplest way to accomplish this? 
>>>>
>>>> The easiest way is to create a class derived from 
>>>> VectorFucntionFromScalarFunctionObject in your own project that 
>>>> implements the 
>>>> missing function. The better way is to write a patch to the deal.II 
>>>> sources 
>>>> that implements the missing function because then others can use your 
>>>> work in 
>>>> the future as well :-) 
>>>>
>>>>
>>>> > On a side note, although "VectorFucntionFromScalarFunctionObject" 
>>>> works with 
>>>> > the L2 norm of the error, it is quite computationally intensive. 
>>>>
>>>> It isn't VectorFucntionFromScalarFunctionObject that's expensive, but 
>>>> the 
>>>> FEFieldFunction. That's sort of inherently the case with this class. If 
>>>> you 
>>>> can, it should be avoided and you should try an evaluate the FE field 
>>>> you have 
>>>> at regular quadrature points, via FEValues. 
>>>>
>>>> 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 a topic in the 
>> Google Groups "deal.II User Group" group.
>> To unsubscribe from this topic, visit 
>> https://groups.google.com/d/topic/dealii/DqH2auneWaY/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to 
>> dealii+un...@googlegroups.com.
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/dd5f93f8-fb58-45c6-9ab2-91ef6c5bdf56n%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/dd5f93f8-fb58-45c6-9ab2-91ef6c5bdf56n%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/bf142002-8af6-492b-a54b-d2fad4093ecen%40googlegroups.com.

[deal.II] Re: How to use cell stiffness matrix to build global stiffness matrix

[Abbas Ballout]

Lance, 

I am not sure if I understood you correctly. 
Maybe looking at the local_assemble_matrix() function in step:13 might 
help. 

Abbas

On Tuesday, August 1, 2023 at 11:00:10 AM UTC+2 dim...@gmail.com wrote:

> Hello dear group,
>
> I have one question which is about how to build global stiffness matrix 
> with cell stiffness matrix.
>
> In our project,the density is a vector (actually we need a mu vector but 
> in the code of deallii mu is a scalar)with different values,I would like to 
> use the element of density vector to build cell matrix as written in the 
> link below:
>
>  
> https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
> (starting from line 1782)
>
> And finally the cell stiffness matrix is used to build global stiffness 
> matrix,but I don' know how the global stiffness matrix is bulit with cell 
> stiffness matrix.
>
>
> https://www.dealii.org/current/doxygen/code-gallery/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc
> (starting from line 1599)
> //>>>>>>>>>>>>>>>>>>>>>>>>>>
>  constraints.distribute_local_to_global(
> data.cell_matrix, 
> data.cell_rhs,
> data.local_dof_indices,
> tangent_matrix, 
> system_rhs);
> }
>
>
> //<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>
>
> I have a idea that I use mu[0] to obtain the first cell_matrix and use 
> mu[1] to get the second cell stiffness matrix ,step by step then I got 
> final mu[m] to obtain the last cell stiffness matrix and in the end I use 
> some function like  distribute_local_to_global to construct the global 
> stiffness matrix,and system rhs and even the solutions.
>
>I checked the code in 
> AffineConstraints::distribute_local_to_global 
> <https://www.dealii.org/current/doxygen/deal.II/classAffineConstraints.html#a373fbdacd8c486e675b8d2bff8943192>.
>  
> (
> https://www.dealii.org/current/doxygen/deal.II/affine__constraints_8h_source.html
> )
> But the hint information guides me to read the cm.templates.h file.
> ewcfp I was a little confused ,I did not find the file from dealii lib. 
> could anyone provide any information or hint? Thanks in advance!
> Best regards
> Lance
> Could 
>
> I
>

-- 
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/1ab39ea3-2fbf-413e-924d-aa8c601766den%40googlegroups.com.

Re: [deal.II] Q: Question about extracting part of a vector

[Abbas Ballout]

Najwaa, 

I submitted a pull request <https://github.com/dealii/dealii/pull/15805> 
recently about something similar 
<https://groups.google.com/g/dealii/c/H0kLG8RJKVc>. (I guess)
Maybe that would help. 

Abbas 
 

On Monday, July 31, 2023 at 7:37:59 PM UTC+2 najwaa...@gmail.com wrote:

> Dear Wolfgang,
>
> Thank you for your answer. I will try to write a patch to the deal.II 
> sources that implement the missing function. This might require some time. 
> I will come back here if needed.
>
>
> Thank you all for your quick answers,
> Najwa
>
> On Wednesday, July 26, 2023 at 11:41:15 PM UTC+3 Wolfgang Bangerth wrote:
>
>> On 7/26/23 12:17, Najwa Alshehri wrote: 
>> > 
>> > Regarding the function "VectorFucntionFromScalarFunctionObject," I have 
>> > observed that it works for computing the L2 norm of the error. However, 
>> when I 
>> > attempted to compute the H1_seminorm of the error, I noticed that the 
>> gradient 
>> > was not implemented for this function. Please let me know if my 
>> understanding 
>> > is incorrect, as I referred to the  for this. 
>>
>> Yes, correct. Someone needs to implement 
>> VectorFucntionFromScalarFunctionObject::gradient() 
>> and/or 
>> VectorFucntionFromScalarFunctionObject::gradient_list() 
>> in exactly the same way as the value() and value_list() functions are 
>> implemented. 
>>
>>
>> > Could you kindly advise me on the simplest way to accomplish this? 
>>
>> The easiest way is to create a class derived from 
>> VectorFucntionFromScalarFunctionObject in your own project that 
>> implements the 
>> missing function. The better way is to write a patch to the deal.II 
>> sources 
>> that implements the missing function because then others can use your 
>> work in 
>> the future as well :-) 
>>
>>
>> > On a side note, although "VectorFucntionFromScalarFunctionObject" works 
>> with 
>> > the L2 norm of the error, it is quite computationally intensive. 
>>
>> It isn't VectorFucntionFromScalarFunctionObject that's expensive, but the 
>> FEFieldFunction. That's sort of inherently the case with this class. If 
>> you 
>> can, it should be avoided and you should try an evaluate the FE field you 
>> have 
>> at regular quadrature points, via FEValues. 
>>
>> 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/dd5f93f8-fb58-45c6-9ab2-91ef6c5bdf56n%40googlegroups.com.

Re: [deal.II] Integrate difference isn't picking on the

[Abbas Ballout]

Lucas 
It's not the linker but I am guessing it's when I edit something deep 
within the library.h files. 
Since everything is linked to it many files are being rebuilt   
I didn't know about other linker so thanks for the tip. I am sure I will do 
that in the future. 

Anyone, 
This is my first pull request ever and I don't know exactly how to approach 
this. 
Details: 
I extended the VectorFunctionFromTensorFunction. to return the 
gradient of a function, 
I ended up editing the files: function.h and function.templates.h. but I 
ended up copying from other files. 

I am getting the right H1 convergence in my example code so I know what I 
am doing works. 
How do I proceed? 

Best,
Abbas  
 

On Saturday, July 22, 2023 at 6:53:24 PM UTC+2 lucasm...@gmail.com wrote:

> I haven't tried this with deal.II yet, but I've recently started using the 
> mold linker for programs I write which have many compilation units (link to 
> GitHub here: https://github.com/rui314/mold). If you're only changing one 
> line then probably only one compilation unit is being recompiled, so the 
> bulk of the wait time is probably being taken by the linker. 
>
> Note that, to get a significant speedup, you'll have to recompile with all 
> of your cores so the linking can be parallelized. I've found that using all 
> cores to compile sometimes causes make to get killed because it uses too 
> much memory, but if most things are compiled and you're just on the linking 
> stage you don't really run into this. 
>
> Alternatively, I think that some other folks use gold linker or the lld 
> linker (in case you can't get mold to work). 
>
> - Lucas
>
>
> On Sat, Jul 22, 2023, 11:37 AM Abbas Ballout  wrote:
>
>> If I change a single line in deal I'll have to wait a good amount of time 
>> for it to compile. 
>> Do you all go through this? All I have is a 4 cores and 16 GBs of RAM on 
>> my Laptop.  
>>
>> On Monday, July 17, 2023 at 9:48:00 AM UTC+2 Abbas Ballout wrote:
>>
>>> YES. Would take me a while but I ll do it. 
>>> Be prepared for many questions. 
>>>
>>> On Monday, July 17, 2023 at 2:44:01 AM UTC+2 Wolfgang Bangerth wrote:
>>>
>>>> On 7/16/23 07:22, Abbas Ballout wrote: 
>>>> > I get this error when I call the 
>>>> exact_solution_vector_function.gradient(p1): 
>>>> > !ERROR Message: 
>>>> > /An error occurred in line <135> of file 
>>>> > 
>>>> 
>>>>  
>>>> in function 
>>>> > dealii::Tensor<1, dim, Number> dealii::Function>>> > RangeNumberType>::gradient(const dealii::Point&, unsigned int) 
>>>> const 
>>>> > [with int dim = 2; RangeNumberType = double] 
>>>> > The violated condition was: 
>>>> > false 
>>>> > Additional information: 
>>>> > You (or a place in the library) are trying to call a function 
>>>> that is 
>>>> > declared as a virtual function in a base class but that has not 
>>>> been 
>>>> > overridden in your derived class./ 
>>>> > / 
>>>> > / 
>>>> > Calling value  works fine but not gradient even though I have 
>>>> overridden both 
>>>> > functions. 
>>>>
>>>> You do in your own class, but VectorFunctionFromTensorFunction does 
>>>> not. 
>>>> That's not by design, but probably simply because nobody has ever had a 
>>>> need 
>>>> for this. Would you like to try and write a patch that adds 
>>>> VectorFunctionFromTensorFunction::gradient()? You could model it on 
>>>> VectorFunctionFromTensorFunction::value(). 
>>>>
>>>> Best 
>>>> Wolfgang 
>>>>
>>>> -- 
>>>>  
>>>>
>>>> Wolfgang Bangerth email: bang...@colostate.edu 
>>>> www: http://www.math.colostate.edu/~bangerth/ 
>>>>
>>>>
>>>> -- 
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see 
>> https://groups.google.com/d/forum/dealii?hl=en
>> --- 
>> You received this message because you are subscribed to the Google Groups 
>> "deal.II User Group" group.
>> To unsubscribe from this group and stop receiving emails from it, send an 
>> email to dealii+un...@googlegroups.com.
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/15f0b1e7-d5ea-49ed-b58f-ea10fae84d52n%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/15f0b1e7-d5ea-49ed-b58f-ea10fae84d52n%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/8ba758f4-06dc-4664-a715-72a7e60ff7bdn%40googlegroups.com.

Re: [deal.II] Integrate difference isn't picking on the

[Abbas Ballout]

If I change a single line in deal I'll have to wait a good amount of time 
for it to compile. 
Do you all go through this? All I have is a 4 cores and 16 GBs of RAM on my 
Laptop.  

On Monday, July 17, 2023 at 9:48:00 AM UTC+2 Abbas Ballout wrote:

> YES. Would take me a while but I ll do it. 
> Be prepared for many questions. 
>
> On Monday, July 17, 2023 at 2:44:01 AM UTC+2 Wolfgang Bangerth wrote:
>
>> On 7/16/23 07:22, Abbas Ballout wrote: 
>> > I get this error when I call the 
>> exact_solution_vector_function.gradient(p1): 
>> > !ERROR Message: 
>> > /An error occurred in line <135> of file 
>> > 
>> 
>>  
>> in function 
>> > dealii::Tensor<1, dim, Number> dealii::Function> > RangeNumberType>::gradient(const dealii::Point&, unsigned int) 
>> const 
>> > [with int dim = 2; RangeNumberType = double] 
>> > The violated condition was: 
>> > false 
>> > Additional information: 
>> > You (or a place in the library) are trying to call a function that 
>> is 
>> > declared as a virtual function in a base class but that has not 
>> been 
>> > overridden in your derived class./ 
>> > / 
>> > / 
>> > Calling value  works fine but not gradient even though I have 
>> overridden both 
>> > functions. 
>>
>> You do in your own class, but VectorFunctionFromTensorFunction does not. 
>> That's not by design, but probably simply because nobody has ever had a 
>> need 
>> for this. Would you like to try and write a patch that adds 
>> VectorFunctionFromTensorFunction::gradient()? You could model it on 
>> VectorFunctionFromTensorFunction::value(). 
>>
>> Best 
>> Wolfgang 
>>
>> -- 
>>  
>> Wolfgang Bangerth email: bang...@colostate.edu 
>> www: http://www.math.colostate.edu/~bangerth/ 
>>
>>
>>

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

Re: [deal.II] How to set refinement flag on the neighbouring cell when it is owned by another processor?

[Abbas Ballout]

I ended up passing the assemble_own_interior_faces_both flag to mesh-worker 
mesh loop. 
I will probably have to use these Consensus algorithms sooner or later 
anyway. 

Thank you. 
Abbas 

On Thursday, July 13, 2023 at 3:17:59 AM UTC+2 Wolfgang Bangerth wrote:

> On 7/11/23 08:23, Abbas Ballout wrote:
> > I don't want to just set a refinement flag for my cell but I also want 
> to set 
> > the refinement
> > flag for the neighbouring cell too and that is what 
> ncellset_refine_flag(); 
> > doesThis doesn't work if the neighbouring cell is owned by another core
>
> Correct. That's because you can't set information on cells you don't own, 
> in 
> particular because it isn't clear what should happen if different 
> non-owning 
> processes (for which this cell is a ghost cell) disagree in their 
> decisions on 
> what they want to set.
>
> If you want to do this kind of thing, you probably want to collect a list 
> of 
> (cell_id, value) that you want to send to each of the neighboring 
> processes, 
> then send it via Utilities::MPI::isend(). Then you'd receive this kind of 
> list 
> on each process via Utilities::MPI::irecv(), unpack it, and then set the 
> value 
> on each of the cells in question.
>
> If you want a higher-level interface, you can use the functions in 
> Utilities::MPI::ConsensusAlgorithms.
>
> Best
> Wolfgang
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>

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

Re: [deal.II] Integrate difference isn't picking on the

[Abbas Ballout]

YES. Would take me a while but I ll do it. 
Be prepared for many questions. 

On Monday, July 17, 2023 at 2:44:01 AM UTC+2 Wolfgang Bangerth wrote:

> On 7/16/23 07:22, Abbas Ballout wrote:
> > I get this error when I call the 
> exact_solution_vector_function.gradient(p1):
> > !ERROR Message:
> > /An error occurred in line <135> of file 
> > 
> 
>  
> in function
> > dealii::Tensor<1, dim, Number> dealii::Function > RangeNumberType>::gradient(const dealii::Point&, unsigned int) 
> const 
> > [with int dim = 2; RangeNumberType = double]
> > The violated condition was:
> > false
> > Additional information:
> > You (or a place in the library) are trying to call a function that is
> > declared as a virtual function in a base class but that has not been
> > overridden in your derived class./
> > /
> > /
> > Calling value  works fine but not gradient even though I have overridden 
> both 
> > functions.
>
> You do in your own class, but VectorFunctionFromTensorFunction does not. 
> That's not by design, but probably simply because nobody has ever had a 
> need 
> for this. Would you like to try and write a patch that adds 
> VectorFunctionFromTensorFunction::gradient()? You could model it on
> VectorFunctionFromTensorFunction::value().
>
> Best
> Wolfgang
>
> -- 
> 
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>

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

Re: [deal.II] Integrate difference isn't picking on the

[Abbas Ballout]

In hindsight I could've written something simpler. 
Here it is:


#include 
#include 

using namespace dealii;

template 
class Exact : public TensorFunction<1, dim>
{
public:
Exact() : TensorFunction<1, dim>(dim)
{}
virtual Tensor<1, dim> value(const Point ) const override;
virtual Tensor<2, dim> gradient(const Point ) const override;
};

template 
Tensor<1, dim> Exact::value(const Point ) const
{
return Tensor<1, dim>();
}

template 
Tensor<2, dim> Exact::gradient(const Point ) const
{
return Tensor<2, dim>();
}

int main()
{
const int dim = 2;
Exact exact_solution;
VectorFunctionFromTensorFunction 
exact_solution_vector_function(exact_solution, 
0, dim);

Point p1;
exact_solution_vector_function.value(p1);
exact_solution_vector_function.gradient(p1);

return 0;
}

I get this error when I call the exact_solution_vector_function.gradient
(p1):
!ERROR Message:







*An error occurred in line <135> of file 

 
in functiondealii::Tensor<1, dim, Number> dealii::Function::gradient(const dealii::Point&, unsigned int) const 
[with int dim = 2; RangeNumberType = double]The violated condition was: 
falseAdditional information: You (or a place in the library) are trying 
to call a function that isdeclared as a virtual function in a base 
class but that has not beenoverridden in your derived class.*

Calling value  works fine but not  gradient even though I have overridden 
both functions. 
On Friday, July 14, 2023 at 8:55:34 PM UTC+2 Wolfgang Bangerth wrote:

> On 7/14/23 09:34, Abbas Ballout wrote:
> > // THis throws exception
> > VectorTools::integrate_difference(dof_handler,
> > solution,
> > exact_solution_vector_function,
> > difference_per_cell,
> > quadrature_formula,
> > VectorTools::H1_seminorm);;
>
> Well, what's the error message? A good rule of thumb is that 50% debugging 
> consists of carefully reading what the error message says.
>
> 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/c1157723-d4a3-4c67-96f7-8eb98aa57fd5n%40googlegroups.com.

[deal.II] Integrate difference isn't picking on the

[Abbas Ballout]

I am overriding the value and the gradient member functions in the 
TensorFunction class then similar to step 7 I am using compute_global_error to 
compute the L2: and H1 norms. I have to convert my TensorFunction to a 
vector valued function in order to do this using 
VectorFunctionFromTensorFunction.
I am getting an exception thrown when I ask for the H1 norm saying that the 
gradient operator is not overridden. 

Here is a sketch of the code:


class Exact : public TensorFunction<1, dim>
{
public:
Exact() : TensorFunction<1, dim>(dim)
{
}

virtual Tensor<1, dim> value(const Point ) const override;
virtual Tensor<2, dim> gradient(const Point ) const override;
};


Tensor<1, dim> Exact::value(const Point ) const
{
return Tensor<1, dim>;
}


Tensor<2, dim> Exact::gradient(const Point ) const
{
return Tensor<2, dim>;
}

int main()
{ 
//blablabla tri, dofs, fe..

 // declate then convert tensor function 
Exact exact_solution;
VectorFunctionFromTensorFunction exact_solution_vector_function(
exact_solution, 0, dim);

//L2 works fine
VectorTools::integrate_difference(dof_handler,
solution,
exact_solution_vector_function,
difference_per_cell,
quadrature_formula,
VectorTools::L2_norm);
// THis throws exception 
VectorTools::integrate_difference(dof_handler,
solution,
exact_solution_vector_function,
difference_per_cell,
quadrature_formula,
VectorTools::H1_seminorm);;


}

Doing something as simple as:
Point p1;
exact_solution_vector_function.value(p1); 
exact_solution_vector_function.gradient(p1); 
is raising the exception. 


I attached the code with the cmake below too.


-- 
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/e7a9106b-ca2a-45b0-9d15-e74832198da4n%40googlegroups.com.
<>

Re: [deal.II] Deformation gradient at face quadrature point for current traction

[Abbas Ballout]

Krish,

I could do this following step 18 and how history variables are constructed 
there.
It was straightforward in case you are interested.   



On Wednesday, June 28, 2023 at 7:55:19 PM UTC+2 Wolfgang Bangerth wrote:

> On 6/26/23 00:49, Krish wrote:
> > 
> > I am following step-44 to solve an elasticity problem. I am trying to 
> > implement current traction (as opposed to the dead load used in step-44) 
> > on a part of the boundary. For that I need to use the deformation 
> > gradient and its determinant at the f_q_point(s) (face quadrature 
> > points). I am interested to know what might be a suitable way to get 
> these.
>
> If you have the deformation field as a solution vector, then computing 
> its gradient at quadrature points on faces is easily done using 
> FEFaceValues. Take a look at the computation of "fluxes" here:
>
> https://dealii.org/developer/doxygen/deal.II/step_4.html#PostprocessingWhattodowiththesolution
> In your case, the gradients are going to be tensors, and you will 
> compute them using code such as
>
> std::vector> 
> solution_gradients(face_quadrature_formula.size());
>
> FEValuesExtractors::Vector deformation(...);
> ...
>
> fe_face_values[deformation].get_function_gradients(solution, 
> solution_gradients);
>
> There is a free function called 'determinant()' that can be used to 
> compute the determinant of a dim x dim tensor.
>
> 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/56e8a831-5704-4a41-b9f9-57a1e832ad4an%40googlegroups.com.

Re: [deal.II] Is it normal for the results to vary this much when running in parallel.?

[Abbas Ballout]

This isn't a problem anymore.
At first I thought that what was happening with step 18 and with distort 
random were caused by the same issue.
I am applying BCs weakly so it's okay.  
Thanks for the help. 

On Wednesday, June 28, 2023 at 4:17:55 PM UTC+2 Wolfgang Bangerth wrote:

> On 6/24/23 09:13, Abbas Ballout wrote:
> > **
> > 
> > Bellow I have plotted the rhs before and after applying boundary 
> constraints 
> > when running with 1, 2 and 3 cores.
> > The rhs looks different at the boundary with MPI 3 as opposed to MPIs 1 
> and 2.
>
> Right. The question is whether it is *wrong*.
> 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/0655e9e3-7aaa-4454-8ac6-415b35f20564n%40googlegroups.com.

Re: [deal.II] Is it normal for the results to vary this much when running in parallel.?

[Abbas Ballout]

   
   - About distort_random

"I suspect that when you call distort_random(), the mesh is different 
whether
you run on one or four processes. Have you confirmed that the meshes are
exactly identical?" 
You were right about this one. After importing a GMSH unstructured mesh the 
results are perfect and independent of 
the MPI processes used.


   -   About step 18

I plotted the stress in Paraview along a line for both MPI runs 1 and 3. 
The results matched. 
The MatrixTools::apply_boundary_values(boundary_values, system_matrix, tmp, 
system_rhs, false) 
function appears to be the culprit . 

With MPI 1 
rhs.l2() before calling apply_boundary_values 3032.65
rhs.l2() after calling apply_boundary_values 1.86145e+10  

with MPI 3 
rhs.l2() before calling apply_boundary_values 3032.65
rhs.l2() after calling apply_boundary_values 3.67161e+10

Is it because apply_boundary_values is doing something special to ghost 
nodes? 

Best,
Abbas 






On Friday, June 16, 2023 at 11:44:43 PM UTC+2 Wolfgang Bangerth wrote:

> On 6/12/23 02:57, Abbas Ballout wrote:
> > **
> > 
> > I am running step 18.
> > 
> > This is the output I getting for a single quasi-tatic step with mpirun 1 
> :
> > 
> > /Cycle 0:
> > Number of active cells:   3712 (by partition: 3712)
> > Number of degrees of freedom: 17226 (by partition: 17226)
> > Assembling system... norm of rhs is 1.88062e+10
> > Solver converged in 103 iterations.
> > Updating quadrature point data...
> >   Cycle 1:
> > Number of active cells:   12805 (by partition: 12805)
> > Number of degrees of freedom: 51708 (by partition: 51708)
> > Assembling system... norm of rhs is 1.86145e+10
> > Solver converged in 120 iterations.
> > Updating quadrature point data...
> > Moving mesh../.
> > 
> >  And this is the output I get when mpirun 3:
> > 
> > /Timestep 1 at time 1
> >   Cycle 0:
> > Number of active cells:   3712 (by partition: 1360+1286+1066)
> > Number of degrees of freedom: 17226 (by partition: 6651+5922+4653)
> > Assembling system... norm of rhs is 1.88062e+10
> > Solver converged in 131 iterations.
> > Updating quadrature point data...
> >   Cycle 1:
> > Number of active cells:   12805 (by partition: 4565+4425+3815)
> > Number of degrees of freedom: 51708 (by partition: 19983+17250+14475)
> > Assembling system... norm of rhs is 3.67161e+10
> > Solver converged in 126 iterations.
> > Updating quadrature point data...
> > Moving mesh.../
> > 
> > The L2 norm in cycle 1 is different between runs with mpi 1 and mpi 3. 
> Is this 
> > normal?
>
> This is an interesting observation. It strikes me as a bug.
>
> Have you compared the solution you get after running either 1 or 3 MPI 
> processes? If you visualize the result, or compute some other indicator of 
> the 
> solution, in which ways do they differ?
>
> 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/3b758adf-f6b4-43e9-9dfb-27d94f3d8040n%40googlegroups.com.

[deal.II] Re: Is it normal for the results to vary this much when running in parallel.?

[Abbas Ballout]

Additionally; 

I took step 40, changed the rhs term to be a constant 1 and added a  
flux_calculation() 

function which loops on the boundary and dots the normal with the gradient 
of the solution (akin to calculating forces).
I also added a GridTools::distort_random(0.3, triangulation, true, 500 or 
something);after mesh generation to get an unstructured mesh.  
I run it for a single cycle and there seem to be a discrepancy in the flux 
calculation for an unstructured grid when running with different cores. 
This doesn't happen with a structured grid. 

The output for the flux calculation on a structured mesh with MPI 1 is:
Running with PETSc on 1 MPI rank(s)...
Cycle 0:
   Number of active cells:   64
   Number of degrees of freedom: 289
   Solved in 9 iterations.
Sigma_flux -0.994394

and for MPI 4 is:
Running with PETSc on 4 MPI rank(s)...
Cycle 0:
   Number of active cells:   64
   Number of degrees of freedom: 289
   Solved in 9 iterations.
Sigma_flux -0.994394

But after activating distort_random the output with MPI 1 is:
Running with PETSc on 1 MPI rank(s)...
Cycle 0:
   Number of active cells:   64
   Number of degrees of freedom: 289
   Solved in 8 iterations.
Sigma_flux -0.993997

but with MPI 4 the output is:
Running with PETSc on 4 MPI rank(s)...
Cycle 0:
   Number of active cells:   64
   Number of degrees of freedom: 289
   Solved in 9 iterations.
Sigma_flux -0.994323

The file is also attached bellow. 

On Monday, June 12, 2023 at 10:57:03 AM UTC+2 Abbas Ballout wrote:

> I am running step 18.  
>
> This is the output I getting for a single quasi-tatic step with mpirun 1 : 
>
>  
>
>
>
>
>
>
>
>
>
>
>
> * Cycle 0:Number of active cells:   3712 (by partition: 3712)
> Number of degrees of freedom: 17226 (by partition: 17226)Assembling 
> system... norm of rhs is 1.88062e+10Solver converged in 103 
> iterations.Updating quadrature point data...  Cycle 1:Number of 
> active cells:   12805 (by partition: 12805)Number of degrees of 
> freedom: 51708 (by partition: 51708)Assembling system... norm of rhs is 
> 1.86145e+10Solver converged in 120 iterations.Updating quadrature 
> point data...Moving mesh..*.
>
>  And this is the output I get when mpirun 3: 
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> *Timestep 1 at time 1  Cycle 0:Number of active cells:   3712 (by 
> partition: 1360+1286+1066)Number of degrees of freedom: 17226 (by 
> partition: 6651+5922+4653)Assembling system... norm of rhs is 
> 1.88062e+10Solver converged in 131 iterations.Updating quadrature 
> point data...  Cycle 1:Number of active cells:   12805 (by 
> partition: 4565+4425+3815)Number of degrees of freedom: 51708 (by 
> partition: 19983+17250+14475)Assembling system... norm of rhs is 
> 3.67161e+10Solver converged in 126 iterations.Updating quadrature 
> point data...Moving mesh...*
>
> The L2 norm in cycle 1 is different between runs with mpi 1 and mpi 3. Is 
> this normal? 
>
> I am experiencing the same problem in my code and the results seem to be 
> slightly mpi dependent.  
>
> The code is attached bellow. It's step 18 the only difference is that I 
> only am running a single step 
>

-- 
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/76e2b9fd-853f-4fc2-a8da-02f86ded0c56n%40googlegroups.com.
/* -
 *
 * Copyright (C) 2009 - 2022 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -

 *
 * Author: Wolfgang Bangerth, Texas A University, 2009, 2010
 * Timo Heister, University of Goettingen, 2009, 2010
 */


#include 
#include 
#include 
#include 

#include 

namespace LA
{
#if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
  !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
  using namespace dealii::LinearAlgebraPETSc;

[deal.II] Is it normal for the results to vary this much when running in parallel.?

[Abbas Ballout]

I am running step 18.  

This is the output I getting for a single quasi-tatic step with mpirun 1 : 

 











* Cycle 0:Number of active cells:   3712 (by partition: 3712)
Number of degrees of freedom: 17226 (by partition: 17226)Assembling 
system... norm of rhs is 1.88062e+10Solver converged in 103 
iterations.Updating quadrature point data...  Cycle 1:Number of 
active cells:   12805 (by partition: 12805)Number of degrees of 
freedom: 51708 (by partition: 51708)Assembling system... norm of rhs is 
1.86145e+10Solver converged in 120 iterations.Updating quadrature 
point data...Moving mesh..*.

 And this is the output I get when mpirun 3: 














*Timestep 1 at time 1  Cycle 0:Number of active cells:   3712 (by 
partition: 1360+1286+1066)Number of degrees of freedom: 17226 (by 
partition: 6651+5922+4653)Assembling system... norm of rhs is 
1.88062e+10Solver converged in 131 iterations.Updating quadrature 
point data...  Cycle 1:Number of active cells:   12805 (by 
partition: 4565+4425+3815)Number of degrees of freedom: 51708 (by 
partition: 19983+17250+14475)Assembling system... norm of rhs is 
3.67161e+10Solver converged in 126 iterations.Updating quadrature 
point data...Moving mesh...*

The L2 norm in cycle 1 is different between runs with mpi 1 and mpi 3. Is 
this normal? 

I am experiencing the same problem in my code and the results seem to be 
slightly mpi dependent.  

The code is attached bellow. It's step 18 the only difference is that I 
only am running a single step 

-- 
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/91bc8ba9-605b-416e-a3e2-fa3e606c41f9n%40googlegroups.com.


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

#include 

#include 

#include 
#include 
#include 

namespace Step18
{
  using namespace dealii;


  template 
  struct PointHistory
  {
SymmetricTensor<2, dim> old_stress;
  };



  template 
  SymmetricTensor<4, dim> get_stress_strain_tensor(const double lambda,
   const double mu)
  {
SymmetricTensor<4, dim> tmp;
for (unsigned int i = 0; i < dim; ++i)
  for (unsigned int j = 0; j < dim; ++j)
for (unsigned int k = 0; k < dim; ++k)
  for (unsigned int l = 0; l < dim; ++l)
tmp[i][j][k][l] = (((i == k) && (j == l) ? mu : 0.0) +
   ((i == l) && (j == k) ? mu : 0.0) +
   ((i == j) && (k == l) ? lambda : 0.0));
return tmp;
  }





  template 
  inline SymmetricTensor<2, dim> get_strain(const FEValues _values,
const unsigned int   shape_func,
const unsigned int   q_point)
  {
SymmetricTensor<2, dim> tmp;

for (unsigned int i = 0; i < dim; ++i)
  tmp[i][i] = fe_values.shape_grad_component(shape_func, q_point, i)[i];

for (unsigned int i = 0; i < dim; ++i)
  for (unsigned int j = i + 1; j < dim; ++j)
tmp[i][j] =
  (fe_values.shape_grad_component(shape_func, q_point, i)[j] +
   fe_values.shape_grad_component(shape_func, q_point, j)[i]) /
  2;

return tmp;
  }


  template 
  inline SymmetricTensor<2, dim>
  get_strain(const std::vector> )
  {
Assert(grad.size() == dim, ExcInternalError());

SymmetricTensor<2, dim> strain;
for (unsigned int i = 0; i < dim; ++i)
  strain[i][i] = grad[i][i];

for (unsigned int i = 0; i < dim; ++i)
  for (unsigned int j = i + 1; j < dim; ++j)
strain[i][j] = (grad[i][j] + grad[j][i]) / 2;

return strain;
  }


  Tensor<2, 2> get_rotation_matrix(const std::vector> _u)
  {
const double curl = (grad_u[1][0] - grad_u[0][1]);

const double angle = std::atan(curl);

return Physics::Transformations::Rotations::rotation_matrix_2d(-angle);
  }


  Tensor<2, 3> get_rotation_matrix(const std::vector> _u)
  {
const Tensor<1, 3> curl({grad_u[2][1] - grad_u[1][2],
 grad_u[0][2] - grad_u[2][0],
 grad_u[1][0] - grad_u[0][1]});

const double tan_angle = std::sqrt(curl * curl);
const double angle = std::atan(tan_angle);

if (std::abs(angle) < 1e-9)
  {
static const double 

Re: [deal.II] Re: How can I visit a face exactly once when running with MPI?

[Abbas Ballout]

Thanks Wolfgang, Time, and Bruno for your comments. I hate to say thanks 
too much so sometimes I don't. 

On Wednesday, June 7, 2023 at 4:40:51 PM UTC+2 Wolfgang Bangerth wrote:

> On 6/7/23 01:21, Abbas Ballout wrote:
> > I set the flag when I create the mesh, then I loop on the faces inside a 
> mesh 
> > worker mesh loop and do something special at a flagged face.
> > This has been working without any issues for sometime now. Is this error 
> prone?
>
> This works on a single process. But it will not work if you are running in 
> parallel because clearing the flag for a face on one process does not 
> clear 
> the flag for that face on another process.
>
> 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/61ac0f95-4da4-4d85-9920-7ad643264d2bn%40googlegroups.com.

Re: [deal.II] Re: How can I visit a face exactly once when running with MPI?

[Abbas Ballout]

What if I am setting the user flag with something like this: 

for (const auto  : triangulation.active_cell_iterators())
for (const unsigned int face : cell->face_indices())
if ((cell->face(face)->center is between x1 and x2 )
cell->face(face)->set_user_flag();

I set the flag when I create the mesh, then I loop on the faces inside a 
mesh worker mesh loop and do something special at a flagged face.
This has been working without any issues for sometime now. Is this error 
prone?  
 
On Thursday, April 27, 2023 at 11:29:34 PM UTC+2 Timo Heister wrote:

> I would suggest that you use MeshWorker::mesh_loop() to do this (or
> look at the complicated implementation of it).
>
> On Thu, Apr 27, 2023 at 11:18 AM Bruno Turcksin
>  wrote:
> >
> > Right, it won't clear it everywhere. You would need to use MPI to clear 
> the flag on the other processor. If you would do that you would create a 
> huge amount of messages and your code would be extremely slow. Best, Bruno 
> ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍
> > ZjQcmQRYFpfptBannerStart
> > This Message Is From an External Sender
> > Use caution when opening links or attachments if you do not recognize 
> the sender.
> >
> > ZjQcmQRYFpfptBannerEnd
> > Right, it won't clear it everywhere. You would need to use MPI to clear 
> the flag on the other processor. If you would do that you would create a 
> huge amount of messages and your code would be extremely slow.
> >
> > Best,
> >
> > Bruno
> >
> > Le jeu. 27 avr. 2023 à 10:09, Abbas Ballout  a 
> écrit :
> >>
> >> This works!! Thanks.
> >> I am guessing something like:
> >>
> >> if(cell_neighbour->is_ghost()==true)
> >> cell->face(neighbour_face_number)->clear_user_flag();
> >>
> >> won't clear the user flag everywhere either, because that didn't work.
> >> Right?
> >>
> >> On Thursday, April 27, 2023 at 2:48:23 PM UTC+2 bruno.t...@gmail.com 
> wrote:
> >>>
> >>> Hello,
> >>>
> >>> On Thursday, April 27, 2023 at 7:56:50 AM UTC-4 abbas.b...@gmail.com 
> wrote:
> >>>
> >>> This works with a single MPI but not for more. Is it because the 
> cell->face(face)->clear_user_flag()
> >>> doesn't clear the user flag in the other distributed meshes?
> >>>
> >>>
> >>> That's correct. The flag is only cleared locally not on the other 
> processor.
> >>>
> >>> I can't directly loop over faces because there is no 
> "face->is_locally_owned()"
> >>> I know that meshworker is an option. Is it the only option?
> >>>
> >>> You can check if the cell that shares the face (use this) is a ghost 
> cell (use this). If that's the case, you can ask for the subdomain id of 
> the locally owned cell and the ghost cell (use this). Now you can decide 
> that only the processor with the lower rank can "touch" the face and the 
> other processor should just skip it.
> >>>
> >>> 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/s07h9mCZhq4/unsubscribe.
> >> To unsubscribe from this group and all its topics, send an email to 
> dealii+un...@googlegroups.com.
> >> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/d7b3b433-e421-4b2c-a837-add3bd1d6f2en%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+un...@googlegroups.com.
> > To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/CAGVt9ePNjo-9oOPy5YdfZ-RntzJDe7eu5Zdeyy6vFL8W9Jod_A%40mail.gmail.com
> .
>
>
>
> -- 
> Timo Heister
> http://www.math.clemson.edu/~heister/
>

-- 
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/7d0c1908-a484-4c2c-868f-abb0138fcc76n%40googlegroups.com.

[deal.II] Re: Is it possible to output data in vtk at quadrature points for me to view in something like Paraview?

[Abbas Ballout]

I figured this one out. 
I just wrote the output at quadrature points to a text file that I read 
that in Paraview. 

On Wednesday, May 31, 2023 at 5:08:38 PM UTC+2 Abbas Ballout wrote:

> This one is quick: 
> I want to visualise my data at q_points. 
>
> I know that paraview will interpolate the data but still. 
>
>

-- 
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/9f94cdfb-14fa-4b25-90db-e3a87f227e40n%40googlegroups.com.

[deal.II] Is it possible to output data in vtk at quadrature points for me to view in something like Paraview?

[Abbas Ballout]

This one is quick: 
I want to visualise my data at q_points. 

I know that paraview will interpolate the data but still. 

-- 
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/19534a7a-8562-4efe-91c8-f0ace4ca9eeen%40googlegroups.com.

Re: [deal.II] How are both the gradient components calculated in deal?

[Abbas Ballout]

Thanks for the reply. 

I was under the impression that get_function_gradients calculates the 
gradient based on the \sum_j U_j that are on the face. 



On Monday, May 22, 2023 at 7:42:34 PM UTC+2 Wolfgang Bangerth wrote:

> On 5/22/23 06:21, Abbas Ballout wrote:
> > How does get_function_gradients calculateboth component gradient at a 
> face?
>
> Abbas,
> since
> u_h(x) = \sum_j U_j \varphi_j(x)
> you get the gradient as
> \nabla u_h(x) = \sum_j U_j \nabla\varphi_j(x)
> and the only other thing you need to know is that \varphi_j is defined by 
> mapping shape functions from the reference cell to the real cell.
>
> But I suspect that you have a reason to ask that question. What is it 
> really 
> that you are after?
>
> 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/2f76aeac-0a44-489c-9927-4d35348bb003n%40googlegroups.com.

[deal.II] How are both the gradient components calculated in deal?

[Abbas Ballout]

Excuse this basic question but how is the gradient at a quadrature point on 
a face calculated in deal? 
Take this code snipped for example: 

QGauss<1> face_quadrature_formula(fe.degree + 1);
FEFaceValues<2> fe_face_values(fe,
face_quadrature_formula,
update_values | update_quadrature_points |
update_normal_vectors | update_gradients |
update_JxW_values);

for (const auto  : dof_handler.active_cell_iterators())
for (const auto  : cell->face_iterators())
{
fe_face_values.reinit(cell, face);
const auto _points = fe_face_values.get_quadrature_points();
const unsigned int n_q_points = q_points.size();

std::vector> grads(n_q_points);
fe_face_values.get_function_gradients(solution, grads);

for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)

{
std::cout << "At q_point\t" << q_points[q_point] << "\n";
std::cout << "Solution grads\t" << grads[q_point] << "\n";

std::cout << "\n";
}
}

How does get_function_gradients calculate both component gradient at a face? 

-- 
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/c05938d9-6844-4114-a4d1-e913a470ce3dn%40googlegroups.com.

Re: [deal.II] Simple Query about an example in the code gallary

[Abbas Ballout]

Sorry I deleted this before it had any comments on it. 
I thought it would disappear. Won't happen again. 
 
Abbas 

On Wednesday, April 26, 2023 at 7:18:59 PM UTC+2 Wolfgang Bangerth wrote:

> On 4/26/23 03:34, Abbas wrote:
> > Why is it that in step 40 the locally relevant solution is initialised 
> with 
> > both locally_owned_dofs and locally_relevant_dofs, while in the 
> > Distributed_LDG_Method example in the code gallery the locally relevant 
> > solution is just initialised with just the locally_relevant_dofs only?
>
> The code gallery program uses this function [2/5]:
>
> https://dealii.org/developer/doxygen/deal.II/classTrilinosWrappers_1_1MPI_1_1Vector.html#af74a39fa88174f63767d510416c09e6f
> The tutorial, on the other hand, uses [3/5] immediately below.
>
> Both are reasonable. I would gladly accept a patch that switches the use 
> from 
> [2/5] to [3/5] because that version is clearly more obvious.
>
> 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/a40acae0-391b-47a3-bceb-1e4c5bc40dbdn%40googlegroups.com.

[deal.II] Re: How can I visit a face exactly once when running with MPI?

[Abbas Ballout]

This works!! Thanks. 
I am guessing something like: 

if(cell_neighbour->is_ghost()==true)
cell->face(neighbour_face_number)->clear_user_flag();

won't clear the user flag everywhere either, because that didn't work. 
Right?  

On Thursday, April 27, 2023 at 2:48:23 PM UTC+2 bruno.t...@gmail.com wrote:

> Hello,
>
> On Thursday, April 27, 2023 at 7:56:50 AM UTC-4 abbas.b...@gmail.com 
> wrote: 
>
> This works with a single MPI but not for more. Is it because the cell->
> face(face)->clear_user_flag()
> doesn't clear the user flag in the other distributed meshes? 
>
>
> That's correct. The flag is only cleared locally not on the other 
> processor.
>
> I can't directly loop over faces because there is no 
> "face->is_locally_owned()" 
> I know that meshworker is an option. Is it the only option? 
>
> You can check if the cell that shares the face (use this 
> )
>  
> is a ghost cell (use this 
> ).
>  
> If that's the case, you can ask for the subdomain id of the locally owned 
> cell and the ghost cell (use this 
> ).
>  
> Now you can decide that only the processor with the lower rank can "touch" 
> the face and the other processor should just skip it.
>
> Best,
>
> Bruno
>  
>

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

[deal.II] How can I visit a face exactly once when running with MPI?

[Abbas]

I want to postporcess some data across an interface with feInterface. 
I have a fully::distributed::triangulation. 
I am doing is something that "roughly" looks like this: 

template 
void TSL_prllel::interface_gradients(Tensor<1, dim> &
interface_gradients)
{
interface_gradients = 0;
QGauss face_quadrature_formula(degree + 1);

FEInterfaceValues fe_iv(fe,
face_quadrature_formula,
+many_flags);

Tensor<1, dim> local_interface_gradients;

for (; cell != endc; ++cell)
if (cell->is_locally_owned())
for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face
)
if (cell->face(face)->user_flag_set()) //If I am at the interface
{ 
// I want to visit the face only once so I clear the flag
cell->face(face)->clear_user_flag();

fe_iv.get_jump_in_function_values(ghosted_solution, gradu);
for (unsigned int point = 0; point < n_face_q_points; ++point)
{

local_interface_gradients += gradu[point]xJXW ;

}
}
interface_stress = Utilities::MPI::sum(local_interface_gradients, 
mpi_communicator);
}

I basically set a flag at the faces I want to postprocess at then when I 
visit the face the first time, I clear the flag so the face gets visited 
only once. This works with a single MPI but not for more. Is it because the 
cell->face(face)->clear_user_flag()
doesn't clear the user flag in the other distributed meshes? 
I can't directly loop over faces because there is no 
"face->is_locally_owned()" 
I know that meshworker is an option. Is it the only option? 

-- 
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/b6e47744-49ad-46e7-9e72-14474f35f9f1n%40googlegroups.com.

[deal.II] Simple Query about an example in the code gallary

[Abbas]

Why is it that in step 40 the locally relevant solution is initialised with 
both locally_owned_dofs and locally_relevant_dofs, while in the 
Distributed_LDG_Method 
example in the code gallery the locally relevant solution is just 
initialised with just the  locally_relevant_dofs only? 

-- 
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/0878caf8-5c27-40ed-ad97-341c15821476n%40googlegroups.com.

[deal.II] Re: Error building with cmake on vscode

[Abbas]

Do you have the cmake extension installed in vscode? 
If not give it a try

Abbas

On Saturday, April 22, 2023 at 1:44:41 PM UTC+2 vinay...@gmail.com wrote:

> Hi,
>
> I am trying to build step-1 using cmake on vscode. But vscode throws me 
> this error:
>
> `
> CMake Error at CMakeLists.txt:30 (message):
>
> *** Could not locate a (sufficiently recent) version of deal.II.  ***
>
>
>
> You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake
>
> or set an environment variable "DEAL_II_DIR" that contains this path.
>
> `
> I have installed dealii using spack. I am able to build using the terminal 
> though. What changes do i need to make for building using vscode?
>
> Thanks 
> Vinayak
>
>

-- 
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/ade9af1e-1787-4840-b9db-c49e45dc6cc3n%40googlegroups.com.

[deal.II] Re: How can I get the jump in my solution if my solution vector is distributed?

[Abbas]

It works now!!
Thanks

Abbas

On Thursday, April 20, 2023 at 8:54:44 PM UTC+2 bruno.t...@gmail.com wrote:

> Hello,
>
> You are using a non-ghosted vector. You need to created a ghosted vector 
> and they assign the non-ghosted vector to the ghosted one.
>
> Best,
>
> Bruno
>
> On Thursday, April 20, 2023 at 7:07:09 AM UTC-4 Abbas wrote:
>
>> I am trying to querry for the jump in  my solution with 
>> "get_jump_in_function_values" 
>>  . My solution vector is PETSc distributed. When I run a code that's like 
>> the one below:
>>
>> template 
>> void output_jumps::interface_jump()
>> {
>>
>> QGauss face_quadrature_formula(degree + 1);
>>
>> FEInterfaceValues fe_iv(fe,
>> face_quadrature_formula,
>> UpdateFlags(update_values |
>> update_gradients |
>> update_quadrature_points |
>> update_normal_vectors |
>> update_JxW_values));
>>
>> const unsigned int n_face_q_points = face_quadrature_formula.size();
>>
>> std::vector> gradu(n_face_q_points);
>>
>> typename DoFHandler::active_cell_iterator cell = dof_handler.
>> begin_active(), endc = dof_handler.end();
>>
>>
>> const FEValuesExtractors::Vector displacements(0);
>>
>>
>> for (; cell != endc; ++cell)
>> if (cell->is_locally_owned())
>> for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++
>> face)
>> if (cell->face(face)->user_flag_set())
>> {
>> // I want to visit the face only once
>> cell->face(face)->clear_user_flag();
>>
>> auto cellnb = cell->neighbor(face);
>> int nface_number = cell->neighbor_face_no(face);
>>
>> fe_iv.reinit(cell, face, numbers::invalid_unsigned_int, cellnb, 
>> nface_number, numbers::invalid_unsigned_int);
>>
>> // This line breaks running code in parallel
>> fe_iv[displacements].get_jump_in_function_values(current_solution, jumpu
>> );
>> }
>> }
>>
>> I get the error: 
>> "You tried to access element 512 of a distributed vector, but only
>> elements in range [0,511] are stored locally and can be accessed."
>>
>> How can I access the ghost entries to calculate the jump in my solution?
>>
>>

-- 
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/162e0355-30c2-4b6e-a487-38cfd5e82042n%40googlegroups.com.

[deal.II] How can I get the jump in my solution if my solution vector is distributed?

[Abbas]

I am trying to querry for the jump in  my solution with 
"get_jump_in_function_values" 
 . My solution vector is PETSc distributed. When I run a code that's like 
the one below:

template 
void output_jumps::interface_jump()
{

QGauss face_quadrature_formula(degree + 1);

FEInterfaceValues fe_iv(fe,
face_quadrature_formula,
UpdateFlags(update_values |
update_gradients |
update_quadrature_points |
update_normal_vectors |
update_JxW_values));

const unsigned int n_face_q_points = face_quadrature_formula.size();

std::vector> gradu(n_face_q_points);

typename DoFHandler::active_cell_iterator cell = dof_handler.
begin_active(), endc = dof_handler.end();


const FEValuesExtractors::Vector displacements(0);


for (; cell != endc; ++cell)
if (cell->is_locally_owned())
for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face
)
if (cell->face(face)->user_flag_set())
{
// I want to visit the face only once
cell->face(face)->clear_user_flag();

auto cellnb = cell->neighbor(face);
int nface_number = cell->neighbor_face_no(face);

fe_iv.reinit(cell, face, numbers::invalid_unsigned_int, cellnb, nface_number, 
numbers::invalid_unsigned_int);

// This line breaks running code in parallel
fe_iv[displacements].get_jump_in_function_values(current_solution, jumpu);
}
}

I get the error: 
"You tried to access element 512 of a distributed vector, but only
elements in range [0,511] are stored locally and can be accessed."

How can I access the ghost entries to calculate the jump in my solution?

-- 
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/81bb0095-2088-4546-a69c-a88fe3e82c38n%40googlegroups.com.

Re: [deal.II] Why does putting the "TimerOuput" before the "ConditionalOStream" leads to a seg fault?

[Abbas]

Goes without saying thanks Prof 

On Wednesday, April 19, 2023 at 4:50:03 PM UTC+2 Wolfgang Bangerth wrote:

> On 4/19/23 08:33, Abbas B. wrote:
> > **
> > 
> > I attached a simple .cc file that illustrates the tittle.
> > Switching lines 18 and 19 leads to a segmentation fault and I would like 
> to 
> > know why out of curiosity.
>
> Abbas:
> This is a C++ question. Member variables are initialized in the order in 
> which 
> they are declared. If the initialization of one variable depends on 
> another, 
> you need to make sure that they are initialized in the correct order, by 
> declaring them in the correct order.
>
> 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/17470a0d-c4e6-4f8c-a533-d7c8e00103c4n%40googlegroups.com.

[deal.II] Why does putting the "TimerOuput" before the "ConditionalOStream" leads to a seg fault?

[Abbas B.]

I attached a simple .cc file that illustrates the tittle. 
Switching lines 18 and 19 leads to a segmentation fault and I would like to 
know why out of curiosity. 
Would appreciate it. Thanks 

-- 
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/86dd356e-58a4-4e61-96e0-e9a0ecfd7236n%40googlegroups.com.


#include 
using namespace dealii;

template 
class LaplaceProblem
{
public:
LaplaceProblem();

void run();

private:
void setup_system();

MPI_Comm mpi_communicator;
TimerOutput computing_timer;
ConditionalOStream pcout;

};

template 
LaplaceProblem::LaplaceProblem()
: mpi_communicator(MPI_COMM_WORLD), pcout(std::cout,
  (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)),
  computing_timer(mpi_communicator,
  pcout,
  TimerOutput::never,
  TimerOutput::wall_times)
{
}

template 
void LaplaceProblem::setup_system()
{
TimerOutput::Scope t(computing_timer, "setup");
}
template 
void LaplaceProblem::run()
{
setup_system();
computing_timer.print_summary();
computing_timer.reset();
}

int main(int argc, char *argv[])
{

Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
LaplaceProblem<2> laplace_problem_2d;
laplace_problem_2d.run();
}

[deal.II] How does the function "evaluate_vector_field()" for example work?

[Abbas]

Hello, 

Not exactly an FEM or dealii question but here goes. 
I am not interested in knowing how the function "evaluate_vector_field" 
works per se, but I want to know how I can figure it out for myself from 
looking at the source code. 

The function "evaluate_vector_field" is a good example. It is not exactly 
documented how it works. I grep searched the include and the source files 
for the function but I can't find the hard code that does the work.

I can ask "How does data post-processor calculate the gradients at the cell 
vertices?" but I want try figure out things like that on my own before 
coming here. 

So, where is the hard code associated with the function 
"evaluate_vector_field()". I can't find it. 

Thanks of course, 
Abbas  

-- 
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/25c133d0-376b-4f7a-b560-57334556a92fn%40googlegroups.com.

Re: [deal.II] VectorTools::interpolate_boundary_values doesn't work for DG

[Abbas]

Thank you Prof. Bangerth. 
I am always greatfull for your input. 
On Friday, March 10, 2023 at 10:27:49 PM UTC+1 Wolfgang Bangerth wrote:

> On 3/8/23 07:49, Abbas wrote:
> > 
> > I am solving a nonlinear problem with DG and my BCs are in the weak 
> form. It 
> > would be helpful if I start with a solution vector that is closer to the 
> > prescribed solution so I won't have to do iterations otherwise would 
> have been 
> > unnecessary.
> > For a purely linear problem, I need to do two non_linear iterations to 
> > converge instead of 1.
> > If I can't use VectorTools::interpolate_boundary_values what are my 
> options? 
> > Does VectorTools::project_boundary_values work in that case?
>
> No, it also needs to know which DoFs are located on the boundary, but for 
> DG 
> elements no DoFs are logically on the boundary.
>
> I don't see a good solution to your problem that uses what's already in 
> the 
> library. I imagine that using code similar to what's in 
> interpolate_boundary_values() but using a different way of finding out 
> which 
> DoF's support point is physically (though not logically) located on the 
> boundary, this all could be implemented in 200 lines of code. But there is 
> nothing simple the obvious already provides you for this case.
>
> 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/4334127b-3826-404a-b4e4-8180b0def2b5n%40googlegroups.com.

Re: [deal.II] VectorTools::interpolate_boundary_values doesn't work for DG

[Abbas]

I would still want to initialise my solution vector to the values which are 
prescribed at the boundary before solving.   
Context: 
I am solving a nonlinear problem with DG and my BCs are in the weak form. 
It would be helpful if I start with a solution vector that is closer to the 
prescribed solution so I won't have to do iterations otherwise would have 
been unnecessary.
For a purely linear problem, I need to do two non_linear iterations to 
converge instead of 1. 
If I can't use VectorTools::interpolate_boundary_values what are my 
options? Does VectorTools::project_boundary_values work in that case? 

On Friday, January 6, 2023 at 7:29:17 PM UTC+1 Wolfgang Bangerth wrote:

>
> > It's cool how dealii is consitent with the maths.
> > Maybe @Wolfgang would agree that this should have thrown an exception in 
> Debug?
>
> In principle yes, but I don't know how to write such an assertion given 
> the 
> generality of the function (for the hp and non-hp case, and when selecting 
> only a subset of vector components). I will see if I can at least add 
> assertions for some simple cases.
>
> 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/582227f6-425d-4a06-acbc-5db23dcea79bn%40googlegroups.com.

Re: [deal.II] Queries on calculation on average stress.

[Abbas]

Daniel, 

Side question: 
If I have to output gradients of a variable, the straightforward approach 
would be to evaluate the gradients over quadrature points and then average. 
This gives me an average value of my gradient per cell.  
Do you know how does data DataPostProcessor evaluates the value of the 
gradient at the mesh nodes? How does it it work under the hood? Do you mind 
explaining?  

Deepika, 

This is a strain post processor. Create a function that does this: 
template 
class StrainPostprocessor : public DataPostprocessorTensor
{
public:
StrainPostprocessor()
: DataPostprocessorTensor("strain",
update_gradients)
{
}

virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector 
_data, std::vector> _quantities) const 
override
{
AssertDimension(input_data.solution_gradients.size(),
computed_quantities.size());

for (unsigned int p = 0; p < input_data.solution_gradients.size(); ++p)
{
AssertDimension(computed_quantities[p].size(),
(Tensor<2, dim>::n_independent_components));

for (unsigned int d = 0; d < dim; ++d)
for (unsigned int e = 0; e < dim; ++e)
computed_quantities[p][Tensor<2, dim>::component_to_unrolled_index(
TableIndices<2>(d, e))] = (input_data.solution_gradients[p][d][e] +
input_data.solution_gradients[p][e][d]) /
2;
}
}
};

Then in data output you can do this: 
void CG_Elasticity_Linear::output_results(unsigned int 
refinement_level, unsigned int cycle) const
{

std::vector
interpretation(dim,
DataComponentInterpretation::component_is_part_of_vector);
std::vector solution_names(dim, "u");

DataOut data_out;

data_out.add_data_vector(dof_handler, solution, solution_names, 
interpretation);
const StrainPostprocessor strain;
data_out.add_data_vector(dof_handler, solution, strain);
data_out.build_patches();

std::ofstream output("solution" + std::to_string(refinement_level) + "_" + 
std::to_string(cycle) + ".vtu");

data_out.write_vtu(output);

}
For averages you'll have to loop in the DataPostProcessor derived class and 
average your quantity.  


Best,
Abbas

On Thursday, March 2, 2023 at 12:42:46 PM UTC+1 d.arnd...@gmail.com wrote:

> Deepika,
>
> > How can I calculate the average stress and strain values for the 
> following geometry (Where element size is not uniform) using the Data 
> PostProcedure?
>
> Can you describe what you are trying to do in a formula? DataPostProcessor 
> typically gives you a result per grid point in the output step.
>
> > 2. Is it possible to make a uniform mesh so the size of elements in both 
> (inclusion and matrix) will be the same?
>
> If you can draw a mesh that you would like to use, you should also be able 
> to create that in deal.II to use. Of course, it will be very difficult to 
> achieve the same size for all elements since you will need more cells to 
> describe the inclusions anyway.
>
> Best,
> Daniel 
>
> On Wed, Mar 1, 2023 at 6:38 PM Deepika Kushwah  
> wrote:
>
>> Hello Everyone,
>>
>> I have the following queries :
>> 1. How can I calculate the average stress and strain values for the 
>> following geometry (Where element size is not uniform) using the Data 
>> PostProcedure?
>>
>> [image: image.png]
>> 2. Is it possible to make a uniform mesh so the size of elements in both 
>> (inclusion and matrix) will be the same?
>>
>> Thanks & Regards
>> Deepika
>>
>> 
>> **
>> This e-mail is for the sole use of the intended recipient(s) and may
>> contain confidential and privileged information. If you are not the
>> intended recipient, please contact the sender by reply e-mail and destroy
>> all copies and the original message. Any unauthorized review, use,
>> disclosure, dissemination, forwarding, printing or copying of this email
>> is strictly prohibited and appropriate legal action will be taken. 
>> 
>>  
>>
>> -- 
>> 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/CAKTdrpp2Ckswjqi8RsU9DDwm6AkT-d05Cin-D33bWmfLv20Ytw%40mail.gmail.com
>>  
>> <https://groups.google.com/d/msgid/dealii/CAKTdrpp2Ckswjqi8RsU9DDwm6AkT-d05Cin-D33bWmfLv20Ytw%40mail.gmail

Re: [deal.II] SIPG is misbehaving at the corners which separate a Dirichlet and a Neumann BC .

[Abbas]

I don't think you were harsh nor did take it that way. That was very 
considerate of you nonetheless. 
It was a light sleezy comment on my end. I should be more careful. 
Thanks again. 
On Wednesday, January 11, 2023 at 3:47:52 PM UTC+2 luca@gmail.com wrote:

> I didn’t mean to sound too harsh, sorry for that. I could spot the mistake 
> right away because it is a mistake I do all the time. 
>
> No embarrassment is needed, and no question is ever a stupid question. 
> Don’t worry! I’ve been in your shoes so many times I cannot count them 
> anymore… 
>
> :) 
>
> Luca.
>
> > On 11 Jan 2023, at 14:33, Abbas  wrote:
> > 
> > Well, this is embarrassing.
> > Thank you very much 
> > 
> > On Wednesday, January 11, 2023 at 3:11:11 PM UTC+2 luca@gmail.com 
> wrote:
> > You have an issue in your boundary worker. 
> > 
> > 
> > The boundary worker is called passing a cell, a face_no, scratch, and 
> copy. Yet, in your code, you loop again on the faces of the cell: 
> > 
> > 
> > for (const auto  : cell->face_iterators()) 
> > { 
> > if ( ((face->boundary_id() == 0) || (face->boundary_id() == 1))) 
> > 
> > 
> > this should become 
> > 
> > if((cell->face(face_no)->boundary_id() == 0) || 
> (cell->face(face_no)->boundary_id() == 1) { 
> > 
> > … 
> > 
> > } 
> > 
> > otherwise, you are going to assemble the face terms twice on the four 
> corners. 
> > 
> > Best, 
> > Luca. 
> > 
> > > On 11 Jan 2023, at 12:05, Abbas  wrote: 
> > > 
> > > Hello, 
> > > 
> > > Setting: 
> > > 1) Unit square domain, 
> > > 2) Simple scalar Laplacian with no source term. 
> > > 3) 0 and 1 Dirichlet on opposite ends of the square and zero Neumann 
> on the other two edges. 
> > > 
> > > I am trying to solve the above problem using SIPG by adopting step 74. 
> The expected solution should be perfectly linear with a constant gradient 
> as the problem reduces to a pseudo 1-D problem. 
> > > 
> > > That is the case if I use CG, but with DG I do not get a perfectly 
> linear solution and that is causing me problems. 
> > > 
> > > I have attached below a single image file that includes: 
> > > 1) The solution. 
> > > 2) The solution gradient with weird behaviour at the corners 
> > > 3) Line plot of the solution across the lower wall. 
> > > 
> > > and the code I am using which is based on step 74. 
> > > 
> > > Things I tried/looked into: 
> > > - I tried to adapt the code to CG and I get the expected results. I 
> can't straight up impose BCs strongly with DG here. 
> > > - I am sure this isn't a Paraview issue because I tried integrating 
> across the boundary within deal and I am getting mesh dependent results. 
> > > - The weak form is correct. The Neuman BC shows up as a boundary 
> integral in the rhs which has no contribution. (Neumann zero at these 
> walls). 
> > > - I looked into the distributed LDG method within the dealii gallary 
> and adopted it to these BCs and I didn't have this problem. 
> > > 
> > > The only thing I haven't tried yet is to do a simple problem by hand 
> and compare the dealii and my outputs. Before doing that, any ideas on how 
> I can debug this? Is this something associated with imposing BCs weakly 
> with SIPG and will happen no matter what? 
> > > 
> > > best, 
> > > Abbas 
> > > 
> > > 
> > > 
> > > -- 
> > > 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/df72814c-0cdc-4215-819c-46f5a8d05e84n%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+un...@googlegroups.com.
> > To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/35d379b6-4297-457e-90d3-14173b5c234dn%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/2c593f0e-9dad-4a50-a233-761e4f8b896an%40googlegroups.com.

Re: [deal.II] SIPG is misbehaving at the corners which separate a Dirichlet and a Neumann BC .

[Abbas]

Well, this is embarrassing.
Thank you very much   

On Wednesday, January 11, 2023 at 3:11:11 PM UTC+2 luca@gmail.com wrote:

> You have an issue in your boundary worker. 
>
>
> The boundary worker is called passing a cell, a face_no, scratch, and 
> copy. Yet, in your code, you loop again on the faces of the cell:
>
>
> for (const auto  : cell->face_iterators())
> {
> if ( ((face->boundary_id() == 0) || (face->boundary_id() == 1)))
>
>
> this should become
>
> if((cell->face(face_no)->boundary_id() == 0) || 
> (cell->face(face_no)->boundary_id() == 1) {
>
> …
>
> }
>
> otherwise, you are going to assemble the face terms twice on the four 
> corners.
>
> Best,
> Luca.
>
> > On 11 Jan 2023, at 12:05, Abbas  wrote:
> > 
> > Hello, 
> > 
> > Setting: 
> > 1) Unit square domain, 
> > 2) Simple scalar Laplacian with no source term. 
> > 3) 0 and 1 Dirichlet on opposite ends of the square and zero Neumann on 
> the other two edges. 
> > 
> > I am trying to solve the above problem using SIPG by adopting step 74. 
> The expected solution should be perfectly linear with a constant gradient 
> as the problem reduces to a pseudo 1-D problem. 
> > 
> > That is the case if I use CG, but with DG I do not get a perfectly 
> linear solution and that is causing me problems. 
> > 
> > I have attached below a single image file that includes: 
> > 1) The solution.
> > 2) The solution gradient with weird behaviour at the corners 
> > 3) Line plot of the solution across the lower wall. 
> > 
> > and the code I am using which is based on step 74. 
> > 
> > Things I tried/looked into: 
> > - I tried to adapt the code to CG and I get the expected results. I 
> can't straight up impose BCs strongly with DG here. 
> > - I am sure this isn't a Paraview issue because I tried integrating 
> across the boundary within deal and I am getting mesh dependent results. 
> > - The weak form is correct. The Neuman BC shows up as a boundary 
> integral in the rhs which has no contribution. (Neumann zero at these 
> walls). 
> > - I looked into the distributed LDG method within the dealii gallary and 
> adopted it to these BCs and I didn't have this problem. 
> > 
> > The only thing I haven't tried yet is to do a simple problem by hand and 
> compare the dealii and my outputs. Before doing that, any ideas on how I 
> can debug this? Is this something associated with imposing BCs weakly with 
> SIPG and will happen no matter what? 
> > 
> > best,
> > Abbas 
> > 
> > 
> > 
> > -- 
> > 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/df72814c-0cdc-4215-819c-46f5a8d05e84n%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/35d379b6-4297-457e-90d3-14173b5c234dn%40googlegroups.com.

Re: [deal.II] collect2 linking error when executing make

[Abbas]

Simon, I believe you shouldn't pass the "-DDEAL_II_WITH_HDF5=OFF" to cmake 
to begin with.  

On Thursday, January 5, 2023 at 6:32:27 PM UTC+2 Simon wrote:

> deal.II was installed using spack.
>
> Do I have to recompile the whole library with 
> $ cmake -DDEAL_II_WITH_HDF5=OFF?
>
>
>
> Am Do., 5. Jan. 2023 um 17:15 Uhr schrieb Timo Heister  >:
>
>> How did you install deal.II on your new machine? Note that the HDF5
>> configuration is part of deal.II and not your user project.
>>
>> On Wed, Jan 4, 2023 at 4:48 AM Simon Wiesheier
>>  wrote:
>> >
>> > I called cmake as follows $ cmake -DDEAL_II_WITH_HDF5=OFF . and get the 
>> additional lines //No help, variable specified on the command line. 
>> DEAL_II_WITH_HDF5: UNINITIALIZED=OFF in my CMakeCache. txt file. However, 
>> after $ make ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍
>> > ZjQcmQRYFpfptBannerStart
>> > This Message Is From an External Sender
>> > Use caution when opening links or attachments if you do not recognize 
>> the sender.
>> >
>> > ZjQcmQRYFpfptBannerEnd
>> > I called cmake as follows
>> > $ cmake -DDEAL_II_WITH_HDF5=OFF .
>> >
>> > and get the additional lines
>> >
>> > //No help, variable specified on the command line.
>> > DEAL_II_WITH_HDF5:UNINITIALIZED=OFF
>> >
>> > in my CMakeCache.txt file.
>> >
>> > However, after
>> > $ make
>> >
>> > the error message is unchanged.
>> > Best
>> > Simon
>> >
>> >
>> > Am Di., 3. Jan. 2023 um 20:55 Uhr schrieb Simon Wiesheier <
>> simon.w...@gmail.com>:
>> >>
>> >> I started from an empty source directory with only a CMakeLists.txt 
>> file.
>> >>
>> >> I do not need hdf5.
>> >> At least I was not aware that I had it on my old computer.
>> >>
>> >> Is there a way to work around the issue?
>> >>
>> >> Best
>> >> Simon
>> >>
>> >>
>> >> Timo Heister  schrieb am Di., 3. Jan. 2023, 20:19:
>> >>>
>> >>> Simon,
>> >>>
>> >>> Make sure that you do not transfer any temporary files
>> >>> (CMakeCache.txt, CMakeFiles/ etc.) but that you start from a fresh
>> >>> source directory. The error message hints at the fact that your new
>> >>> system does not have hdf5 installed.
>> >>>
>> >>> On Mon, Jan 2, 2023 at 4:21 AM Simon  wrote:
>> >>> >
>> >>> > Dear all, First of all, I wish you a happy and prosperous new year 
>> 2023! I transferred my dealii project to a different computer and encounter 
>> an linking error (command make) after setting up the project (command cmake 
>> .). ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍ ‍
>> >>> > ZjQcmQRYFpfptBannerStart
>> >>> > This Message Is From an External Sender
>> >>> > Use caution when opening links or attachments if you do not 
>> recognize the sender.
>> >>> >
>> >>> > ZjQcmQRYFpfptBannerEnd
>> >>> > Dear all,
>> >>> >
>> >>> > First of all, I wish you a happy and prosperous new year 2023!
>> >>> >
>> >>> > I transferred my dealii project to a different computer and
>> >>> > encounter an linking error (command make) after setting up
>> >>> > the project (command cmake .).
>> >>> > The error message is attached.
>> >>> >
>> >>> > What is the problem?
>> >>> >
>> >>> > Best
>> >>> > Simon
>> >>> >
>> >>> > --
>> >>> > The deal.II project is located at http://www.dealii.org/
>> >>> > For mailing list/forum options, see 
>> https://groups.google.com/d/forum/dealii?hl=en
>> >>> > ---
>> >>> > You received this message because you are subscribed to the Google 
>> Groups "deal.II User Group" group.
>> >>> > To unsubscribe from this group and stop receiving emails from it, 
>> send an email to dealii+un...@googlegroups.com.
>> >>> > To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/ca0ce1f2-cba0-42e1-8ee4-2178806540ebn%40googlegroups.com
>> .
>> >>>
>> >>>
>> >>>
>> >>> --
>> >>> Timo Heister
>> >>> http://www.math.clemson.edu/~heister/
>> >>>
>> >>> --
>> >>> 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/CAMRj59HYLvENRpM7%3Dwy4v%3D8pvbkm7viGjMa4Apr3ytoZKX%3DktQ%40mail.gmail.com
>> .
>> >
>> > --
>> > The deal.II project is located at http://www.dealii.org/
>> > For mailing list/forum options, see 
>> https://groups.google.com/d/forum/dealii?hl=en
>> > ---
>> > You received this message because you are subscribed to the Google 
>> Groups "deal.II User Group" group.
>> > To unsubscribe from this group and stop receiving emails from it, send 
>> an email to dealii+un...@googlegroups.com.
>> > To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/CAM50jEsJVDPE-6Y1tuviCXU%2BmgPD55%3D1pX%3DyW8JTP5eghM7ZWA%40mail.gmail.com
>> .
>>
>>
>>
>> -- 
>> Timo Heister
>> 

Re: [deal.II] VectorTools::interpolate_boundary_values doesn't work for DG

[Abbas]

Thank you Daniel. 

I was using penalties before taking this route, because I had artifacts at 
the boundary faces present at a corner between a Dirichlet and a Neuman BC 
edge. A problem I didn't have with CG.  My first intuition was to apply BCs 
strongly to fix this, but I don't think it will even if I were able to 
interpolate. 

It's cool how dealii is consitent with the maths.
Maybe @Wolfgang would agree that this should have thrown an exception in 
Debug?

On Wednesday, January 4, 2023 at 9:31:34 PM UTC+2 d.arnd...@gmail.com wrote:

> Abbas,
>
> Yes, FE_DGQ only defines dofs in the interior of cells. Thus, trying to 
> interpolate the boundary has no effect. You would normally use penalty 
> terms to prescribe boundary values for DG methods.
>
> Best,
> Daniel
>
> On Wed, Jan 4, 2023 at 6:04 PM Abbas  wrote:
>
>> Hello, 
>>
>> The VectorTools::interpolate_boundary_values function works for 
>> interpolating fe_Q functions but not fe_DGQ . 
>> Was this meant to be the case? 
>> I attached a minimal code that illustrates this.
>>
>>
>> -- 
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see 
>> https://groups.google.com/d/forum/dealii?hl=en
>> --- 
>> You received this message because you are subscribed to the Google Groups 
>> "deal.II User Group" group.
>> To unsubscribe from this group and stop receiving emails from it, send an 
>> email to dealii+un...@googlegroups.com.
>> To view this discussion on the web visit 
>> https://groups.google.com/d/msgid/dealii/0a8d1da6-e564-45de-bc3d-1f5a9d5dfad9n%40googlegroups.com
>>  
>> <https://groups.google.com/d/msgid/dealii/0a8d1da6-e564-45de-bc3d-1f5a9d5dfad9n%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/a8f4659b-c4e8-4a3d-81fc-ed802c7bd3bcn%40googlegroups.com.

[deal.II] VectorTools::interpolate_boundary_values doesn't work for DG

[Abbas]

Hello, 

The VectorTools::interpolate_boundary_values function works for 
interpolating fe_Q functions but not fe_DGQ . 
Was this meant to be the case? 
I attached a minimal code that illustrates this.
   

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

#include 
#include 

#include 
#include 

#include 
#include 
#include 

#include 
using namespace dealii;
int main()

{
// create an feq obkect and interpolate a value of 1 at boundary 0
{
Triangulation<2> tria;
GridGenerator::hyper_cube(tria, 0, 1, true);
tria.refine_global(4);
FE_Q<2> fe(1);
DoFHandler<2> dof(tria);

Vector solution;

dof.distribute_dofs(fe);
solution.reinit(dof.n_dofs());

std::map boundary_values;
VectorTools::interpolate_boundary_values(dof,
 0,
 Functions::ConstantFunction<2>(1),
 boundary_values);
for (auto _value : boundary_values)
solution(boundary_value.first) = boundary_value.second;

DataOut<2> data_out;

data_out.attach_dof_handler(dof);
data_out.add_data_vector(solution, "solution");
data_out.build_patches();

const std::string filename = "solution-feq.vtu";
std::ofstream output(filename);
data_out.write_vtu(output);
}

// Same as above but
// we create an feGQq obkect and interpolate
{
Triangulation<2> tria;
GridGenerator::hyper_cube(tria, 0, 1, true);
tria.refine_global(4);
FE_DGQ<2> fe(1); // only difference between this and the above code 
DoFHandler<2> dof(tria);

Vector solution;

dof.distribute_dofs(fe);
solution.reinit(dof.n_dofs());

std::map boundary_values;
VectorTools::interpolate_boundary_values(dof,
 0,
 Functions::ConstantFunction<2>(1),
 boundary_values);
for (auto _value : boundary_values)
solution(boundary_value.first) = boundary_value.second;

DataOut<2> data_out;

data_out.attach_dof_handler(dof);
data_out.add_data_vector(solution, "solution");
data_out.build_patches();

const std::string filename = "solution-fedgq.vtu";
std::ofstream output(filename);
data_out.write_vtu(output);
}
}

Re: [deal.II] Re: Error in Reading Gmsh File in Dealii

[Abbas]

If you get stuck you can set the material ID manually after importing it 
with something like:

Point p ;
for (const auto  : triangulation.active_cell_iterators())
  {p=cell->center() ; 
   if(sqrt(p[0]^2+p[1]^2)set_material_id(1) ;
else
cell->set_material_id(2) ;
   }

It's a cheat and might be a bad idea if your geometry gets complicated 
though. 
On Thursday, December 22, 2022 at 5:23:35 PM UTC+2 deepika...@iitgoa.ac.in 
wrote:

> Thank you for your response.
> I have changed the 2d algorithm from tool -> option-> mesh. Now it is 
> showing that the mesh curvilinear is off. Still I am getting the same error 
> while reading the gmsh file in dealii. Please find the attached updated 
> mesh and the error message.
>
> [image: image.png]
> [image: image.png]
>
> Thanks
> Deepika
>
> On Thu, Dec 22, 2022 at 7:41 PM Bruno Turcksin  
> wrote:
>
>> Deepika,
>>
>> Can you try using a linear mesh instead of curvilinear one? Does that 
>> work?
>>
>> Best,
>>
>> Bruno
>>
>> On Thursday, December 22, 2022 at 1:47:32 AM UTC-5 
>> deepika...@iitgoa.ac.in wrote:
>>
>>> Hello Everyone,
>>>
>>> I have successfully read a gmsh file for one material in my code. Now I 
>>> am trying to read the gmsh file which has two different materials, but it 
>>> is showing the following error (see below). I am not able to understand the 
>>> meaning of the error message. Please help.
>>>
>>> [image: image.png]
>>> The gmsh geometry is also attached for reference.
>>> [image: image.png]
>>>
>>> Thanks & Regards,
>>> Deepika
>>>
>>> 
>>> **
>>> This e-mail is for the sole use of the intended recipient(s) and may
>>> contain confidential and privileged information. If you are not the
>>> intended recipient, please contact the sender by reply e-mail and destroy
>>> all copies and the original message. Any unauthorized review, use,
>>> disclosure, dissemination, forwarding, printing or copying of this email
>>> is strictly prohibited and appropriate legal action will be taken. 
>>> 
>>>  
>>>
>> -- 
>> 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/819f4a60-a64f-414c-b1a2-cc980f9c6b4bn%40googlegroups.com
>>  
>> 
>> .
>>
>
> **
> This e-mail is for the sole use of the intended recipient(s) and may
> contain confidential and privileged information. If you are not the
> intended recipient, please contact the sender by reply e-mail and destroy
> all copies and the original message. Any unauthorized review, use,
> disclosure, dissemination, forwarding, printing or copying of this email
> is strictly prohibited and appropriate legal action will be taken. 
> 
>  
>

-- 
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/0fa36010-55b0-4097-b0c9-cf0f8d3b17c5n%40googlegroups.com.

[deal.II] Re: MatrixFree support for Nedelec element in deal.II

[Abbas]

Nedelec elements aren't cheap either. If I recall correctly you'd need 108 
dofs for the 6 EM components per element to achieve first order 
convergence. (I hope someone would correct me on this if I am wrong) 
With AMR and hanging nodes this can get expensive fairly quickly.  The only 
thing that worries me about DG is if they could handle the div free 
condition well, and if they are robust for high frequency problems. 
Just my 2 cents. Hope it helps in some way even if you don't chose this 
path. 

best, 
Abbas  

On Monday, December 5, 2022 at 4:45:12 PM UTC+2 yy.wayne wrote:

> Thanks Abbas,
>
> I considered about DG method, because DG allows discontinuity and 
> implements boundary condition easier than FE_Q. 
> However DG requires more DoF than continuous finite element, especially in 
> frequency-domain. 
>
> Best,
> Wayne
>
> 在2022年12月5日星期一 UTC+8 22:35:12 写道：
>
>> Maybe you can look into using DG for discritizing the problem? You're 
>> already using penalty terms so might as well. 
>> One thing I am afraid of is that Nedelec elements inherently satisfy the 
>> div free condition, while DG might need careful treatment or your 
>> preconditioner would fail at high frequency. 
>> Just my 2 cents. Maybe you know someone else who can help you with this.  
>>  
>>
>>
>> On Monday, December 5, 2022 at 4:29:35 PM UTC+2 yy.wayne wrote:
>>
>>> Hi deal.II community,
>>>
>>> I'm working on MatrixFree electromagnetic problems with deal.II. For 2d 
>>> TE waves FE_Q element works good, but for cases like in step-81 and in 3d, 
>>> Nedelec element in the appropriate one for EM.
>>>
>>> I posted a discussion about this 3 months ago. May I ask what's the 
>>> current state for MatriFree support of Nedelec elements in deal.II? 
>>> Implementing MatrixFree Nedelec is quite difficult for me. Currently I plan 
>>> to solve 3d EM problem with FE_Q elements with some penalty terms, but of 
>>> course that's not a nice policy because of spurious solution and 
>>> singularity when solving EM numerically with nodal elements.
>>>
>>

-- 
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/887683e8-eda3-4032-afc3-3be956ce7b26n%40googlegroups.com.

Re: [deal.II] How to create the magnitude of my solution in a new Vector?

[Abbas]

I figured it out. 
I need to use history variables as in step 18.  

On Tuesday, December 13, 2022 at 9:18:27 AM UTC+2 Abbas wrote:

> EXACTLY. 
> However, it is not clear to me how can I store that value, so that in the 
> next iteration I would be able to compare it with the current value at the 
> q-point, and take the maximum between the two ,and then store that for the 
> next iteration and so on.   
>
>
> On Tuesday, December 13, 2022 at 2:24:40 AM UTC+2 Wolfgang Bangerth wrote:
>
>> On 12/12/22 05:41, Abbas wrote: 
>> > 
>> > I am solving a vector problem similar to step-8. I have 
>> > aVector solution with u_x and u_y at each dof. I would like to 
>> have 
>> > another Vector magnitude which stores the magnitude 
>> \sqrt(u_x^2+u_y^2) 
>> > at each dof. Now it is clear that solution has double the size of 
>> magnitude. I 
>> > need to extract the values of the magnitude vector for my weak form 
>> within the 
>> > assembly loop not for post-processing. 
>> > 
>> > Why do I need this? Well I don't actually want the magnitude exactly, 
>> but what 
>> > I need is the maximum value of my magnitude when I iterate in time. 
>>
>> That sounds very much like you want to solve a nonlinear system -- for 
>> which 
>> step-15 shows you how to access this sort of information. 
>>
>> The point to note regarding step-15, though, is that it never actually 
>> creates 
>> a vector that stores what in your case is the magnitude. That's because 
>> the 
>> magnitude you show above is not actually a finite element field even 
>> though 
>> u_x and u_y are: Even though the latter are bilinear functions (assuming 
>> you 
>> are using a Q1 element), the magnitude is not a binlinear function. As a 
>> consequence, you don't want to compute it at each node and then 
>> interpolate it 
>> at quadrature points; rather, you want to compute it at each quadrature 
>> point 
>> from the original solutions u_x and u_y. 
>>
>> 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/480b6b41-c894-4e05-b4e3-a91480d21caan%40googlegroups.com.

Re: [deal.II] How to create the magnitude of my solution in a new Vector?

[Abbas]

EXACTLY. 
However, it is not clear to me how can I store that value, so that in the 
next iteration I would be able to compare it with the current value at the 
q-point, and take the maximum between the two ,and then store that for the 
next iteration and so on.   


On Tuesday, December 13, 2022 at 2:24:40 AM UTC+2 Wolfgang Bangerth wrote:

> On 12/12/22 05:41, Abbas wrote:
> > 
> > I am solving a vector problem similar to step-8. I have 
> > aVector solution with u_x and u_y at each dof. I would like to 
> have 
> > another Vector magnitude which stores the magnitude 
> \sqrt(u_x^2+u_y^2) 
> > at each dof. Now it is clear that solution has double the size of 
> magnitude. I 
> > need to extract the values of the magnitude vector for my weak form 
> within the 
> > assembly loop not for post-processing.
> > 
> > Why do I need this? Well I don't actually want the magnitude exactly, 
> but what 
> > I need is the maximum value of my magnitude when I iterate in time.
>
> That sounds very much like you want to solve a nonlinear system -- for 
> which 
> step-15 shows you how to access this sort of information.
>
> The point to note regarding step-15, though, is that it never actually 
> creates 
> a vector that stores what in your case is the magnitude. That's because 
> the 
> magnitude you show above is not actually a finite element field even 
> though 
> u_x and u_y are: Even though the latter are bilinear functions (assuming 
> you 
> are using a Q1 element), the magnitude is not a binlinear function. As a 
> consequence, you don't want to compute it at each node and then 
> interpolate it 
> at quadrature points; rather, you want to compute it at each quadrature 
> point 
> from the original solutions u_x and u_y.
>
> 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/4b56f1ca-2263-4aaf-a025-e1e2c71ad142n%40googlegroups.com.

[deal.II] How to create the magnitude of my solution in a new Vector?

[Abbas]

Hello, 

This might be a simple question but I am stuck. 

I am solving a vector problem similar to step-8. I have a 
Vector solution with u_x and u_y at each dof. I would like to have 
another Vector magnitude which stores the magnitude 
\sqrt(u_x^2+u_y^2) at each dof. Now it is clear that solution has double 
the size of magnitude. I need to extract the values of the magnitude vector 
for my weak form within the assembly loop not for post-processing. 

Why do I need this? Well I don't actually want the magnitude exactly, but 
what I need is the maximum value of my magnitude when I iterate in time.  

 

-- 
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/da49a96a-bb94-49ae-92ca-3e009bc76ccbn%40googlegroups.com.

[deal.II] Re: How would you document a step/code/module?

[Abbas]

Thank you Marc for your input, 
Maybe I am just complicating things. 
If I ever make a good example/module in the future, I'll make sure to share 
it to the code gallery. 

On Monday, December 5, 2022 at 1:11:14 AM UTC+2 mafe...@gmail.com wrote:

> Hello Abbas,
>
> if you are willing the share your code with the deal.II community, you can 
> submit your code to the deal.II code gallery. The documentation for each of 
> the example programs use the same doxygen formatting as deal.II.
> https://www.dealii.org/code-gallery.html
> https://github.com/dealii/code-gallery
>
> Doxygen in general is a powerful tool to document your code, but other 
> tools exist as well as for example Jupyter as you pointed out. I feel like 
> it boils down to personal taste which one you would like to use. But kudos 
> to you in thinking about documenting and sharing your code!
>
> Marc
>
> On Sunday, December 4, 2022 at 6:26:41 AM UTC-7 Abbas wrote:
>
>> Hello, 
>>
>> Let's say you made 1000+ lines of code with dealii, and you want to 
>> document and share it in a format similar to that of the dealii steps, how 
>> would you go about doing this?
>>
>> MAYBE in the future people will share single "modules" with some 
>> "tutorial" like documentation. This might be something some of us are 
>> thinking about or considering. I guess? 
>>
>> I have been suggested to use a jupyter notebook but I feel like that's a 
>> huge waste of white space. Having a dedicated web link is also very 
>> attractive so I was considering having a website hosted on github pages. I 
>> tried that using Gatsby and it, but I feel like this is unnecessarily 
>> over-complicated.
>>
>> Has anyone done something similar to this? If not, any ideas?
>>
>> Thank you for reading this far.  
>>
>

-- 
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/2dde0e1c-8394-457d-8ec1-0a2976354630n%40googlegroups.com.

[deal.II] Re: MatrixFree support for Nedelec element in deal.II

[Abbas]

Nedelec aren't cheap either. You need an fe(n+1) to achieve n accuracy. 
It's 64 dofs for all the 6 components of the electric field E in 3D to 
achieve 1st order accuracy. 
With adaptive mesh refinement, Nedelec will get a lot more expensive than 
this. 
By all means, do what seems fit. I am just trying to be a bit informative. 
Maybe you get one of those questions in the future.   

best,
Abbas 

On Monday, December 5, 2022 at 4:45:12 PM UTC+2 yy.wayne wrote:

> Thanks Abbas,
>
> I considered about DG method, because DG allows discontinuity and 
> implements boundary condition easier than FE_Q. 
> However DG requires more DoF than continuous finite element, especially in 
> frequency-domain. 
>
> Best,
> Wayne
>
> 在2022年12月5日星期一 UTC+8 22:35:12 写道：
>
>> Maybe you can look into using DG for discritizing the problem? You're 
>> already using penalty terms so might as well. 
>> One thing I am afraid of is that Nedelec elements inherently satisfy the 
>> div free condition, while DG might need careful treatment or your 
>> preconditioner would fail at high frequency. 
>> Just my 2 cents. Maybe you know someone else who can help you with this.  
>>  
>>
>>
>> On Monday, December 5, 2022 at 4:29:35 PM UTC+2 yy.wayne wrote:
>>
>>> Hi deal.II community,
>>>
>>> I'm working on MatrixFree electromagnetic problems with deal.II. For 2d 
>>> TE waves FE_Q element works good, but for cases like in step-81 and in 3d, 
>>> Nedelec element in the appropriate one for EM.
>>>
>>> I posted a discussion about this 3 months ago. May I ask what's the 
>>> current state for MatriFree support of Nedelec elements in deal.II? 
>>> Implementing MatrixFree Nedelec is quite difficult for me. Currently I plan 
>>> to solve 3d EM problem with FE_Q elements with some penalty terms, but of 
>>> course that's not a nice policy because of spurious solution and 
>>> singularity when solving EM numerically with nodal elements.
>>>
>>

-- 
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/1f4bd525-930b-4966-a0f9-f88b0a5edefbn%40googlegroups.com.

[deal.II] Re: MatrixFree support for Nedelec element in deal.II

[Abbas]

Maybe you can look into using DG for discritizing the problem? You're 
already using penalty terms so might as well. 
One thing I am afraid of is that Nedelec elements inherently satisfy the 
div free condition, while DG might need careful treatment or your 
preconditioner would fail at high frequency. 
Just my 2 cents. Maybe you know someone else who can help you with this.   


On Monday, December 5, 2022 at 4:29:35 PM UTC+2 yy.wayne wrote:

> Hi deal.II community,
>
> I'm working on MatrixFree electromagnetic problems with deal.II. For 2d TE 
> waves FE_Q element works good, but for cases like in step-81 and in 3d, 
> Nedelec element in the appropriate one for EM.
>
> I posted a discussion about this 3 months ago. May I ask what's the 
> current state for MatriFree support of Nedelec elements in deal.II? 
> Implementing MatrixFree Nedelec is quite difficult for me. Currently I plan 
> to solve 3d EM problem with FE_Q elements with some penalty terms, but of 
> course that's not a nice policy because of spurious solution and 
> singularity when solving EM numerically with nodal elements.
>

-- 
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/9ca08a1c-b9ab-46cb-8e90-1f65fad47b76n%40googlegroups.com.

Re: [deal.II] How do I set initial solution vector to boundary values in DG?

[Abbas]

Thank you Profs. 

Prof. Bangerth, I probably wasn't clear in my question. I know how to apply 
BCs for an interior penalty problem weakly. 
But, for a non-linear problem like step 15, one initialises the solution 
vector and sets it to the value of the BCs at the boundary before solving 
for the Newton update. This is done explicitly in function 
'set_boundary_values()' in step 15 with the function  '
VectorTools::interpolate_boundary_values 
<https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#ab2562d41bb26f362043f9719a8cd9b87>
'. 
I cannot use this function with DG elements apparently and I am wondering 
if there are any alternatives.  

Prof. Hiester, don't you think that subtracting small and manipulating 
large numbers in the residual will cause some problems? 
I'll give it a try for sure though. 

Thanks again 


On Monday, December 5, 2022 at 3:01:32 AM UTC+2 Timo Heister wrote:

> I would consider the boundary conditions part of the nonlinear
> iteration (in contrast to step-15, where they are strongly enforced
> and as such are only needed in step 0). Every step you can evaluate
> your nonlinear residual which contains a residual in the boundary
> condition and that difference needs to be applied in every step
> (assuming you solve for an update in each Newton step).
>
> On Sat, Dec 3, 2022 at 10:56 AM Wolfgang Bangerth
>  wrote:
> >
> > On 12/3/22 07: 11, Abbas wrote: > > Something like this doesn't work 
> with DG. > Loosely speaking, my initial approach would be to solve a system 
> with just the > BC terms but I am not sure. Do I have other options? In DG 
> methods,
> > ZjQcmQRYFpfptBannerStart
> > This Message Is From an External Sender
> > Use caution when opening links or attachments if you do not recognize 
> the sender.
> >
> > ZjQcmQRYFpfptBannerEnd
> >
> > On 12/3/22 07:11, Abbas wrote:
> > >
> > > Something like this doesn't work with DG.
> > > Loosely speaking, my initial approach would be to solve a system with 
> just the
> > > BC terms but I am not sure. Do I have other options?
> >
> > In DG methods, you impose boundary values weakly, using the same 
> strategy with
> > which you impose continuity between cells weakly. You might want to look 
> at
> > papers on DG methods: pretty much every single one will show the jump 
> terms
> > corresponding to Dirichlet boundary values.
> >
> > Best
> > W.
> >
> > --
> > 
> > Wolfgang Bangerth email: bang...@colostate.edu
> > www: 
> https://urldefense.com/v3/__http://www.math.colostate.edu/*bangerth/__;fg!!PTd7Sdtyuw!Vf-IQx7HQQ5bIq8xBQu1BzCdpZiz9Q3UC6emGlg1Z50xBT-oA9UeUBtyamyADtrTV0zLEZpa0RBQ_lV6jh5B8Uo$
> >
> >
> > --
> > The deal.II project is located at 
> https://urldefense.com/v3/__http://www.dealii.org/__;!!PTd7Sdtyuw!Vf-IQx7HQQ5bIq8xBQu1BzCdpZiz9Q3UC6emGlg1Z50xBT-oA9UeUBtyamyADtrTV0zLEZpa0RBQ_lV6CDq7NmU$
> > For mailing list/forum options, see 
> https://urldefense.com/v3/__https://groups.google.com/d/forum/dealii?hl=en__;!!PTd7Sdtyuw!Vf-IQx7HQQ5bIq8xBQu1BzCdpZiz9Q3UC6emGlg1Z50xBT-oA9UeUBtyamyADtrTV0zLEZpa0RBQ_lV6qJBbcDU$
> > ---
> > 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://urldefense.com/v3/__https://groups.google.com/d/msgid/dealii/c1d413d2-3c55-0523-19fd-9f399d1781da*40colostate.edu__;JQ!!PTd7Sdtyuw!Vf-IQx7HQQ5bIq8xBQu1BzCdpZiz9Q3UC6emGlg1Z50xBT-oA9UeUBtyamyADtrTV0zLEZpa0RBQ_lV6CppyxvM$
>  
> .
>
>
>
> -- 
> Timo Heister
> http://www.math.clemson.edu/~heister/
>

-- 
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/e6b575d7-9bfd-4137-b30b-3e3d2a335144n%40googlegroups.com.

[deal.II] How would you document a step/code/module?

[Abbas]

Hello, 

Let's say you made 1000+ lines of code with dealii, and you want to 
document and share it in a format similar to that of the dealii steps, how 
would you go about doing this?

MAYBE in the future people will share single "modules" with some "tutorial" 
like documentation. This might be something some of us are thinking about 
or considering. I guess? 

I have been suggested to use a jupyter notebook but I feel like that's a 
huge waste of white space. Having a dedicated web link is also very 
attractive so I was considering having a website hosted on github pages. I 
tried that using Gatsby and it, but I feel like this is unnecessarily 
over-complicated.

Has anyone done something similar to this? If not, any ideas?

Thank you for reading this far.  

-- 
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/00ee1344-6c78-4115-946e-2e741632b6c6n%40googlegroups.com.

Re: [deal.II] How can I integrate along a surface/edge within the domain? Can I tag the surface/edge with an ID?

[Abbas Ballout]

Thank you.This works.
 For anyone seeing this in the future, change "int face" to "int f" in the
for loop.

-- 
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/CAHZ_B5f28FCY%3DCrRnDrbX6Bu5bfdXd-OUGNfhUKYGOdqiXoJvg%40mail.gmail.com.

[deal.II] How do I set initial solution vector to boundary values in DG?

[Abbas]

Hello, 

I am solving a non-linear equation with DG. 
As per step 15, we set the boundary values initially using the following 
function: 

template 
void MinimalSurfaceProblem::set_boundary_values()
{
  std::map boundary_values;
  VectorTools::interpolate_boundary_values 

(dof_handler, 0, BoundaryValueshttp://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/9e893d7a-2629-4378-aa8e-ce457f0c0e51n%40googlegroups.com.

[deal.II] How can I integrate along a surface/edge within the domain? Can I tag the surface/edge with an ID?

[Abbas]

Hello, 

As per the screenshot below, I created and imported a mesh using gmsh, with 
two material IDs. I want to integrate on the edge that separate the two 
materials for my weak form. I tried tagging the internal edge with an ID in 
gmsh to loop on it, but of course dealii throws an error saying that I 
cannot tag an internal edge with an ID. 

How can I integrate on an arbitrary edge like this one? 

-- 
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/2f2f050b-ca85-43b4-a88b-9829070a846bn%40googlegroups.com.

[deal.II] Re: Doubt Related to DataPostprocessor class

[Abbas]

Change "data_out.write_vtk 

 (output);" 
to "data_out.write_vt 

u (output);" 
https://groups.google.com/g/dealii/c/T1ckYL0Le_U/m/54r575TiAwAJ


On Wednesday, November 23, 2022 at 11:16:42 AM UTC+2 
deepika...@iitgoa.ac.in wrote:

> Hello Everyone,
>
> I am trying to calculate the gradient of displacement for step 8, for 
> that, I am using the Datapostprocessor Tensor class (same example code 
> which is given in the documentation here - 
> https://www.dealii.org/current/doxygen/deal.II/classDataPostprocessorTensor.html)
>  
> but I am getting the following error:
>
>
>  66%] Linking CXX executable step-8
> [ 66%] Built target step-8
> [100%] Run step-8 with Debug configuration
> Cycle 0:
>Number of active cells:   256
>Number of degrees of freedom: 578
>
>
> 
> Exception on processing: 
>
> 
> An error occurred in line <5286> of file <./source/base/data_out_base.cc> 
> in function
> void dealii::DataOutBase::write_vtk(const 
> std::vector >&, const 
> std::vector >&, const 
> std::vector std::__cxx11::basic_string, 
> std::allocator >, 
> dealii::DataComponentInterpretation::DataComponentInterpretation> >&, const 
> dealii::DataOutBase::VtkFlags&, std::ostream&) [with int dim = 2; int 
> spacedim = 2; std::ostream = std::basic_ostream]
> The violated condition was: 
> std::get<3>(nonscalar_data_range) != 
> DataComponentInterpretation::component_is_part_of_tensor
>
>
>
>
>
>
>
>
>
> *Additional information: You are trying to use functionality in 
> deal.II that is currently notimplemented. In many cases, this indicates 
> that there simply didn'tappear much of a need for it, or that the 
> author of the original codedid not have the time to implement a 
> particular case. If you hit thisexception, it is therefore worth the 
> time to look into the code tofind out whether you may be able to 
> implement the missingfunctionality. If you do, please consider 
> providing a patch to thedeal.II development sources (see the deal.II 
> website on how tocontribute).*
> 
>
> Aborting!
> 
> make[3]: *** [CMakeFiles/run.dir/build.make:71: CMakeFiles/run] Error 1
> make[2]: *** [CMakeFiles/Makefile2:116: CMakeFiles/run.dir/all] Error 2
> make[1]: *** [CMakeFiles/Makefile2:123: CMakeFiles/run.dir/rule] Error 2
> make: *** [Makefile:137: run] Error 2
>
> How to resolve this? 
>
>
> Thanks & Regards,
> Deepika
>
> **
> This e-mail is for the sole use of the intended recipient(s) and may
> contain confidential and privileged information. If you are not the
> intended recipient, please contact the sender by reply e-mail and destroy
> all copies and the original message. Any unauthorized review, use,
> disclosure, dissemination, forwarding, printing or copying of this email
> is strictly prohibited and appropriate legal action will be taken. 
> 
>  
>

-- 
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/145b97cc-c131-45bc-940c-636fb3774f47n%40googlegroups.com.

Re: [deal.II] Gradient Postprocessor not working

[Abbas]

That was a quick fix and it works.  
Just changed "data_out.write_vtk 
<https://www.dealii.org/current/doxygen/deal.II/group__Exceptions.html#gacad99726038e4fca7f605fdffb3317e4>
 
(output);" to "data_out.write_vt 
<https://www.dealii.org/current/doxygen/deal.II/group__Exceptions.html#gacad99726038e4fca7f605fdffb3317e4>
u (output);" 
Thanks 

Doxygen fixes are done with a PR too? Can I do it?  It will take me forever 
tho. 


On Wednesday, November 9, 2022 at 10:35:32 PM UTC+2 Wolfgang Bangerth wrote:

> On 11/8/22 09:11, Abbas wrote:
> > *** Caution: EXTERNAL Sender ***
> > 
> > So I copied what's in the link 
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassDataPostprocessorTensor.html=05%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cc98d01e4aa464ed9fb3808dac1a3d7e0%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C638035207831485386%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C2000%7C%7C%7C=dGblC2ovakBl3ham3jsHZX2VHuCqOefmWOp%2Fg%2BLja2I%3D=0>
>  verbatim 
>
> > to step 8 but it says that I am "trying to use functionality in deal.II 
> > that is currently not implemented"
>
> The problem doesn't actually have anything to do with the postprocessor, 
> but with the fact that you're outputting your data in VTK format (using 
> data_out.write_vtk(...)). That function does not support tensor output, 
> either because we didn't implement it or because the VTK file format 
> doesn't support tensors -- I don't recall. In any case, switch to 
> write_vtu() instead and everything should work.
>
> 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/7fa490a6-1c8f-4cf7-9155-e930d14320afn%40googlegroups.com.

[deal.II] Gradient Postprocessor not working

[Abbas]

So I copied what's in the link 

 verbatim 
to step 8 but it says that I am "trying to use functionality in deal.II 
that is currently not implemented" 
I attached the code below. 
I'll use the default approach for outputting gradients of a vector; is 
there a dealii step that does this?  

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/a4f172ee-d4b4-4fc9-8347-21ddf33e018bn%40googlegroups.com.
/* -
 *
 * Copyright (C) 2000 - 2021 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -

 *
 * Author: Wolfgang Bangerth, University of Heidelberg, 2000
 */


#include 
#include 
#include 
#include 

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

#include 
#include 
#include 

#include 
#include 

#include 

#include 
#include 
#include 
#include 

#include 
#include 

#include 
#include 

namespace Step8
{
  using namespace dealii;

  template 
class GradientPostprocessor : public DataPostprocessorTensor
{
public:
  GradientPostprocessor ()
:
DataPostprocessorTensor ("grad_u",
  update_gradients)
  {}
 
  virtual
  void
  evaluate_vector_field
  (const DataPostprocessorInputs::Vector _data,
   std::vector > _quantities) const override
  {
// ensure that there really are as many output slots
// as there are points at which DataOut provides the
// gradients:
AssertDimension (input_data.solution_gradients.size(),
 computed_quantities.size());
 
for (unsigned int p=0; p::n_independent_components));
for (unsigned int d=0; d::component_to_unrolled_index(TableIndices<2>(d,e))]
  = input_data.solution_gradients[p][d][e];
  }
  }
};
  template 
  class ElasticProblem
  {
  public:
ElasticProblem();
void run();

  private:
void setup_system();
void assemble_system();
void solve();
void refine_grid();
void output_results(const unsigned int cycle) const;

Triangulation triangulation;
DoFHandlerdof_handler;

FESystem fe;

AffineConstraints constraints;

SparsityPattern  sparsity_pattern;
SparseMatrix system_matrix;

Vector solution;
Vector system_rhs;
  };



  template 
  void right_hand_side(const std::vector> ,
   std::vector> &  values)
  {
AssertDimension(values.size(), points.size());
Assert(dim >= 2, ExcNotImplemented());

Point point_1, point_2;
point_1(0) = 0.5;
point_2(0) = -0.5;

for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
  {
if (((points[point_n] - point_1).norm_square() < 0.2 * 0.2) ||
((points[point_n] - point_2).norm_square() < 0.2 * 0.2))
  values[point_n][0] = 1.0;
else
  values[point_n][0] = 0.0;

if (points[point_n].norm_square() < 0.2 * 0.2)
  values[point_n][1] = 1.0;
else
  values[point_n][1] = 0.0;
  }
  }






  template 
  ElasticProblem::ElasticProblem()
: dof_handler(triangulation)
, fe(FE_Q(1), dim)
  {}



  template 
  void ElasticProblem::setup_system()
  {
dof_handler.distribute_dofs(fe);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());

constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(dof_handler,
 0,
 Functions::ZeroFunction(dim),
 constraints);
constraints.close();

DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler,
dsp,
constraints,
/*keep_constrained_dofs = */ false);
sparsity_pattern.copy_from(dsp);

system_matrix.reinit(sparsity_pattern);
  }



  template 
  

Re: [deal.II] fe_interface_values but of a tensor quantity

[Abbas]

Thanks for the suggestions. 
For the moment I will try to use the identity {uu}={u]{u}+0.25[u][u] and 
then I'll see what happens. 

On Friday, September 9, 2022 at 7:47:29 PM UTC+3 Wolfgang Bangerth wrote:

>
> Abbas,
> I have no specific suggestions to make this easy. The two options I 
> would explore are:
> * Build a tensor-valued (discontinuous) finite element, a DoFHandler for 
> it, and then interpolate the function u*u onto that finite element field 
> -- that is, build a finite element field representation of u*u. You can 
> then evaluate that via FEInterfaceValues.
> * Take a look at how FEInterfaceValues computes jump(u) and average(u) 
> (namely, by evaluating FEValues objects on both sides of the interface), 
> and write code that does something similar for u*u.
>
> The latter approach can be made more general, namely if you wanted to 
> evaluate the jump of average of an arbitrary function f(u). I suspect it 
> wouldn't be terribly hard to generalize the implementation of jump(u) 
> and average(u) to that case, and that might be interesting to others as 
> well if it were put into deal.II itself. If you'd like to explore this 
> and write a patch, let me know and I'll walk you through how one would 
> approach that!
>
> Best
> Wolfgang
>
>
> On 9/9/22 01:08, Abbas wrote:
> > *** Caution: EXTERNAL Sender ***
> > 
> > Yes
> > 
> > On Friday, September 9, 2022 at 12:06:23 AM UTC+3 Wolfgang Bangerth 
> wrote:
> > 
> > On 9/8/22 04:40, Abbas wrote:
> > >
> > > I want the average of (uu) which is = u1*u1/2 +u2*u2/2 but I
> > don't know
> > > how to get it.
> > 
> > Do I understand right that u1 and u2 are vectors, and when you say
> > u1 * u1
> > you mean the outer product of a vector with itself, forming a tensor?
> > 
> > Best
> > W.
> > 
> > -- 
> > 
> > 
> > Wolfgang Bangerth email: bang...@colostate.edu
> > www: http://www.math.colostate.edu/~bangerth/
> > <http://www.math.colostate.edu/~bangerth/>
> > 
> > -- 
> > The deal.II project is located at http://www.dealii.org/ 
> > <
> https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F=05%7C01%7CWolfgang.Bangerth%40colostate.edu%7C0060db9985d343a0ac5c08da923214dd%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637983041060585432%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=g9HAQL%2FlGBi7QI7%2BX9Pgv72TkrXMuHcgzMUpJ1yVW5o%3D=0
> >
> > For mailing list/forum options, see 
> > https://groups.google.com/d/forum/dealii?hl=en 
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den=05%7C01%7CWolfgang.Bangerth%40colostate.edu%7C0060db9985d343a0ac5c08da923214dd%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637983041060585432%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=f0cuUMKC9ftoV6Crbmjr9GIiHXa%2BVkpFsgllUwlHrP4%3D=0
> >
> > ---
> > 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 
> > <mailto:dealii+un...@googlegroups.com>.
> > To view this discussion on the web visit 
> > 
> https://groups.google.com/d/msgid/dealii/b4f5197e-c700-422d-8a04-c776b834cfa2n%40googlegroups.com
>  
> > <
> https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2Fb4f5197e-c700-422d-8a04-c776b834cfa2n%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter=05%7C01%7CWolfgang.Bangerth%40colostate.edu%7C0060db9985d343a0ac5c08da923214dd%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637983041060585432%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=m%2Fu8%2BxOxny0rygj6T1Q2GmJOiY0a0AH%2BfM1DHEXRtf0%3D=0
> >.
>
>
> -- 
> 
> 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/7c075865-7561-46d3-bea3-5a287297e7f2n%40googlegroups.com.

Re: [deal.II] fe_interface_values but of a tensor quantity

[Abbas]

Yes

On Friday, September 9, 2022 at 12:06:23 AM UTC+3 Wolfgang Bangerth wrote:

> On 9/8/22 04:40, Abbas wrote:
> > 
> > I want the average of (uu) which is = u1*u1/2 +u2*u2/2 but I don't know 
> > how to get it.
>
> Do I understand right that u1 and u2 are vectors, and when you say
> u1 * u1
> you mean the outer product of a vector with itself, forming a tensor?
>
> 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/b4f5197e-c700-422d-8a04-c776b834cfa2n%40googlegroups.com.

[deal.II] fe_interface_values but of a tensor quantity

[Abbas]

I am using DG elements and if I want to query the average of my vector 
solution across the interface I can do something like :
  std::vector> old_velocity_averages(n_q_points);
  
fe_iv[velocities].get_average_of_function_values(solution_np0,old_velocity_averages);
which gives me the average of u=u1/2+u2/2

I want the average of (uu) which is = u1*u1/2 +u2*u2/2 but I don't know how 
to get it. 
 
I am stuck and will take any hints. Should I inherit the fe_interface class 
and implement my own function?   

-- 
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/c4a0a408-1d39-4819-b1d8-304f9ec17116n%40googlegroups.com.

[deal.II] On the discretization of the time harmonic Maxwell equation

[Abbas]

This is not a dealii question per se but again, if it's someone who can 
answer this question it's going to be someone from here. 

The literature mostly refers to two approaches todiscretizing the time 
harmonic Maxwell equation. The first one is through the use of edge curl 
conforming Nedelec elements and the other being through the use of interior 
penalty DG. 

Any comments on the use of one approach over the other? 

-- 
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/ace82c1c-42b4-4e76-ade7-c591fc898940n%40googlegroups.com.

[deal.II] The include fe/fe_values header file seems redundant

[Abbas]

Hello,

Just a tiny remark, but it seems that commenting "#include 
"  out of the header yields no errors. I tried that 
in step-3 and the step worked fine. Is the FEvalues class imported by 
another library indirectly? 

Best,
Abbas

-- 
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/48bbd092-407c-47fd-a2d6-8045b0a25ea8n%40googlegroups.com.

[deal.II] Re: SparseDirectUMFPACK - member function "initialize" takes very long

[Abbas]

As for the second question. 
So the number of Dofs isn't the only predictor of how long an LU would 
take. 
Two things to note:
1) In 2D a single node is connected to 4 elements while in 3D a single node 
is connected to 8 elements which means that an unknown will show up in more 
equations which means that your matrix is less sparse now. This makes your 
LU more expensive to compute. 
2) Node numbering plays a role. I am not sure to what extent and I am not 
sure if it is something that is done automatically by direct solvers. 
Step-2 with it's Cuthill_McKee 
<https://www.dealii.org/current/doxygen/deal.II/namespaceDoFRenumbering.html#ab938a690bf4e2adff191fe969b0f21d3>
 is 
something you might want to have a look at .

Best,
Abbas  
On Wednesday, June 2, 2021 at 10:26:29 AM UTC-4 Simon wrote:

> Dear all,
>
> I use the SparseDirectUMFPACK solver for doing some post-processing steps.
> I also measure some cpu_times in my program in order to compare different 
> methods regarding efficiency. For my question only the following three 
> lines of code from my program are relevant:
> ---
> timer.enter_section("Initialize Direct Solver");
> sparse_direct_solver.initialize(global_matrix);
> timer.leave_subsection("Initialize Direct Solver");
> ---
> So I only measure the time which is neeed for executing the member 
> function "initialize":
> ->Solving a scalar problem with *dim=2* and *263425 DoFs* the above call 
> takes about 1.3 seconds, whereas solving the linear system takes only 0.8 
> seconds.
> ->Solving the same scalar problem with *dim=3* and only *72369 DoFs* the 
> above call takes about 4.4 seconds and solving the linear system about 3.2 
> seconds.
>
> -My first question is what this function actually does since it takes 
> longer than assembling and solving the system together? 
> The documentation says: "This function does nothing. It is only here to 
> provide a interface consistent with other sparse direct solvers."
> By having a look in the function body I´ve seen that this function 
> acutally does nothing, i.e. the function body is empty.
>
> -My second question is why for my dim=3 case the function call takes much 
> longer as in the dim=2 case? In the latter I have much more DoFs.
> My guess was that the number of DoFs is the essential quantity but this is 
> obviously not the case.
>
> Any input would be greatly appreciated!
>
> Best
> Simon
>
>

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

[deal.II] Re: Preconditioner for the time harmonic Maxwell equation

[Abbas]

Thank you Sebestian, 

In case someone comes across this thread in the future 
It seems that the problem that is being solved in the "Large scale 3D 
geomag... <https://library.seg.org/doi/10.1190/geo2015-0013.1>" paper is 
different than mine. 
The paper's: curl x curl x E + sigma E = J   , sigma > 0 
Mine is: curl x curl x E  - k^2 E  = 0 
and because of that it looks Hypre's AMS 
<https://hypre.readthedocs.io/en/latest/solvers-ams.html> solver can only 
work with problems as given in the paper
Am I correct here? 

best,
Abbas 

On Wednesday, May 19, 2021 at 12:29:13 PM UTC-4 sebastian...@gmail.com 
wrote:

>
> Hi Abbas,
>
> as you already noticed, standard preconditioners do not yield very good 
> convergence for Maxwell's equations, for more details see Why it is 
> Difficult to Solve Helmholtz Problems with Classical Iterative Methods 
> <https://link.springer.com/chapter/10.1007/978-3-642-22061-6_10>. Since 
> the scalar Helmholtz problem is closely related to the Maxwell's equations, 
> the same arguments from the Helmholtz problem also hold true for the 
> Maxwell's equations. So not even your plan B would work (see Geometric 
> Multigrid Methods for Maxwell’s Equations 
> <https://julianroth.org/res/bachelorarbeit.pdf>, Chapter 6.2).
>
> Therefore you need to apply specially suited preconditioners for the 
> Maxwell's equations. For the problem you are considering, I would recommend 
> the block-preconditioner, which is presented in this paper: Large-scale 
> 3D geoelectromagnetic modeling using parallel adaptive high-order finite 
> element method <https://library.seg.org/doi/10.1190/geo2015-0013.1>. As a 
> starting point for your implementation of that block-preconditioner you can 
> use step-22 <https://www.dealii.org/current/doxygen/deal.II/step_22.html>
> . 
>
> An other approach to solve the Maxwell's equations is via domain 
> decomposition method (for example see: A quasi-optimal domain 
> decomposition algorithm for the time-harmonic Maxwell's equations 
> <https://www.sciencedirect.com/science/article/pii/S0021999115001965>).
>
> I hope this helps you, best regards
> Sebastian
>   
> Abbas schrieb am Mittwoch, 19. Mai 2021 um 15:45:35 UTC+2:
>
>> This isn't a dealii problem and I am sorry if this isn't the right place 
>> to post to, but if there is someone who could give me some hints on how I 
>> would approach this then it would be you guys. 
>>
>> I have been trying to solve the time harmonic maxwell equation 
>> discritized with the use of Nedelec elements. Details are in This stack 
>> overflow post 
>> <https://scicomp.stackexchange.com/questions/37447/how-to-create-a-good-preconditioner-for-a-system-of-linear-equations-that-is-cre>
>> . 
>>
>> I have an issue with the linear system, since using any black box 
>> Trilinos preconditioner yields no convergence or worse convergence than a 
>> simple Jaccobi or an identity preconditioner and even in that case the 
>> restarted GMRES solver needs more than 400 
>> vectors to converge and it gets worse with a refined mesh. My system is 
>> broken into 2x2 block matrices much like step-29  by the way. 
>>
>> Creating a multigrid solver like step-16 might solve this issue but I 
>> would rather use that as plan B. 
>>
>> Any hint is appreciated.  
>>
>> Thanx :3 
>> Abbas 
>>
>>

-- 
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/75d667f1-3b8c-4939-aa61-fa5d199d7648n%40googlegroups.com.

[deal.II] Preconditioner for the time harmonic Maxwell equation

[Abbas Ballout]

This isn't a dealii problem and I am sorry if this isn't the right place to 
post to, but if there is someone who could give me some hints on how I 
would approach this then it would be you guys. 

I have been trying to solve the time harmonic maxwell equation discritized 
with the use of Nedelec elements. Details are in This stack overflow post 
<https://scicomp.stackexchange.com/questions/37447/how-to-create-a-good-preconditioner-for-a-system-of-linear-equations-that-is-cre>
. 

I have an issue with the linear system, since using any black box Trilinos 
preconditioner yields no convergence or worse convergence than a simple 
Jaccobi or an identity preconditioner and even in that case the restarted 
GMRES solver needs more than 400 
vectors to converge and it gets worse with a refined mesh. My system is 
broken into 2x2 block matrices much like step-29  by the way. 

Creating a multigrid solver like step-16 might solve this issue but I would 
rather use that as plan B. 

Any hint is appreciated.  

Thanx :3 
Abbas 

-- 
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/6be07afd-4418-4917-b68f-39cd7c9d1b4dn%40googlegroups.com.

Re: [deal.II] using VectorTools::project_boundary_values_curl_conforming_l2 with FEsystem fe(FE_Nedelec<3>(0), 2)

[Abbas Ballout]

I am running my code in Debug using the latest 9.3 release. 
I have attached a minimal code and what it does it that it prints the 
constraints.

   - If BCs are applied for both parts of the system ,my output for the 
   constraint is:


   1. 0 = 0
   2. 1 = 0
   3. 2 = 0
   4. 3 = 0
   5. 4 = 1.0
   6. 5 = 0
   7. 6 = 1.0
   8. 7 = 0


   - If BCs are applied for the first part of the system, my output for the 
   constraint is:


   1. 0 = 0
   2. 2 = 0
   3. 4 = 1.0
   4. 6 = 1.0


   - If BCs are applied for he first partof the system my output for the 
   constraint is:


   1. 1 = 0
   2. 3 = 0
   3. 5 = 0
   4. 7 = 0




On Monday, May 3, 2021 at 8:22:58 PM UTC-4 Wolfgang Bangerth wrote:

> On 5/3/21 7:43 AM, Abbas Ballout wrote:
> > *The right hand side of the second equation isn't being updated when I 
> apply 
> > BCs even though the system matrix is.
> > *
> > 1) Hints on what I might be missing in my code?
> > 2) Any suggested workarounds in case that there is an actual bug?
>
> That's too far downstream. You have a suspicion that a function isn't 
> working 
> as expected (which is totally plausible), so let's check that right after 
> that 
> function is called. Dealing with the right hand side happens much later in 
> your program, and many other things could have gone wrong in the meantime. 
> I 
> would just print the constraints right after they were computed and check 
> that 
> first.
>
> 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/14582a5f-b419-481b-8129-dee84f4601d2n%40googlegroups.com.
#include 

#include 

#include 
#include 
#include 

#include 
#include 

#include 

#include 

//You will want to set this 
#include "/home/abbas/deal.ii-candi/dealii/dealii/tests/tests.h"

using namespace dealii; 

  class BC_class: public Function<3>
  {
  public:
BC_class(): Function<3>(3)  
{}

void vector_value_list(const std::vector> ,
  std::vector> &  values) const override;
  };

  void BC_class::vector_value_list(const std::vector> ,
  std::vector> &  values) const 
 {
for (unsigned int i = 0; i < points.size(); ++i)
  {
values[i][0] = 1.0;
values[i][1] = 0.0;
values[i][2] = 0.0;
  }
  }
  

int main()
{ Triangulation<3> triangulation ; 
  
  GridGenerator::hyper_cube(triangulation); 
  triangulation.begin_active()->face(4)->set_boundary_id(1);  //set the 4th boundary id to 1 to apply dirichlet BC

  DoFHandler<3>  dof_handler(triangulation);
  FESystem<3>fe(FE_Nedelec<3>(0),2);  
  dof_handler.distribute_dofs(fe);


  AffineConstraints constraints;
  constraints.clear();
  BC_class new_BC; 

  //Apply BC on the first part of the system
  VectorTools::project_boundary_values_curl_conforming_l2(
dof_handler,
0,  //the first 3 compenents
new_BC,
1,  //face ID1
constraints,
StaticMappingQ1<3>::mapping);
//Apply BC on the second part of the system
  VectorTools::project_boundary_values_curl_conforming_l2(
dof_handler,
3,  //the second 3 components 
new_BC,
1,   //face ID1
constraints,
StaticMappingQ1<3>::mapping);

 constraints.close();

initlog();
constraints.print(deallog.get_file_stream());
  
}

[deal.II] using VectorTools::project_boundary_values_curl_conforming_l2 with FEsystem fe(FE_Nedelec<3>(0), 2)

[Abbas Ballout]

Some context: 
I am trying to model a waveguide which would require me to solve the double 
curl equation curlxcurlxE-cE=0 where c is some constant(real or imaginary) 
subject to nxE=nxf on different parts of the boundary where f can be real 
or imaginary also. 

The problem requires the use of Nedelec elements but you already know that. 
Since the problem is complex I am breaking it down into two with FEsystem. 
I am trying to approach this by solving this simple system before 
complicating things:
curlxcurlxE_r-cE_r=0
nxE_r=1   on an "inlet"
nxE_r=0   everywhere else 
and 
curlxcurlxE_i-cE_i=0
nxE_i=1   on an "inlet"
nxE_i=0   everywhere else 

These two systems are uncoupled and exactly the same. The issue has been 
that my BCs aren't being applied for the second part of the equation. I am 
using VectorTools::project_boundary_values_curl_conforming_l2(). The first 
part of the equation gives a solution (inaccurate but that's not my issue 
here) with the BCs applied. 

I print out the sparsity pattern and there seems to be a difference between 
applying BCs and not applying them which makes me conclude that 
project_boundary_values_curl_conforming_l2 is effecting the matrix. 
I printed out the global right hand side for a smaller problem and the 
second half of it, which is the half that is related to the second problem, 
is just zeros so my guess is that 
constraints.distribute_local_to_global(cell_matrix, cell_rhs, 
local_dof_indices, system_matrix, system_rhs) isn't applying the BC in the 
cell_rhs when FENedelec systems are involved. 

I have attached my code here.
Thanks  










-- 
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/8c95e3cb-1643-4cbf-a79f-b7ebc7f02104n%40googlegroups.com.
#include 
#include 

#include 
#include 
#include 
#include 

#include 
#include 
#include 
#include 

#include 
#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; 

// Perfect conductor BC nxE=0
  class PerfConductor : public Function<3>
  {
  public:
PerfConductor(): Function<3>(3)  
{}

void vector_value_list(const std::vector> ,
  std::vector> &  values) const override;
  };

void PerfConductor::vector_value_list(const std::vector> ,
  std::vector> &  values) const 
 {
for (unsigned int i = 0; i < points.size(); ++i)
  {
values[i][0] = 0.0;
values[i][1] = 0.0;
values[i][2] = 0.0;
  }
  }


//Inlet BC nXE=1
  class Inlet: public Function<3>
  {
  public:
Inlet(): Function<3>(3)  
{}

void vector_value_list(const std::vector> ,
  std::vector> &  values) const override;
  };

  void Inlet::vector_value_list(const std::vector> ,
  std::vector> &  values) const 
 {
for (unsigned int i = 0; i < points.size(); ++i)
  {
values[i][0] = 1.0;
values[i][1] = 0.0;
values[i][2] = 0.0;
  }
  }



class Nedlec
{
  public:
  Nedlec();
  void run();

  private:
  void make_grid();
  void setup_system();
  void assemble_system();
  void solve();
  void output_results() const;


  Triangulation<3> triangulation;
  DoFHandler<3>   dof_handler;
  FESystem<3>  fe; 
  AffineConstraints constraints;

  SparsityPattern   sparsity_pattern;
  SparseMatrix  system_matrix;
  Vectorsolution;
  Vectorsystem_rhs;


};

//Constructor
Nedlec::Nedlec()   
: dof_handler(triangulation)
, fe(FE_Nedelec<3>(0), 2)  
{}


//makegrid
void Nedlec::make_grid()
{
  Point<3> p1(0,0,0) , p2 (0.2,0.1,1) ;
  
  GridGenerator::hyper_rectangle(triangulation,p1,p2); 
  triangulation.begin_active()->face(4)->set_boundary_id(1);  //set the 4th boundary id to 1 to apply dirichlet BC
  triangulation.refine_global(2);

  for (const auto  : dof_handler.active_cell_iterators())
  {
cell->set_refine_flag(RefinementCase<3>::cut_z);
  }

  triangulation.execute_coarsening_and_refinement();
  
for (const auto  : dof_handler.active_cell_iterators())
  {
cell->set_refine_flag(RefinementCase<3>::cut_z);
  }


  triangulation.execute_coarsening_and_refinement();
  
  std::ofstream out("grid-1.vtk");
  GridOut   grid_out;
  grid_out.write_vtk(triangulation, out);
  std::cout << "Grid written to grid-1.vtk" << std::endl;
  std::cout << "Number of active cells = " << triangulation.n_active_cells() << std::endl ;
 
}

void