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.
>
> I don't know, and I don't know whether the people who wrote step-74 have 
> thought about the case. One would *hope* that the method also works for 
> the 
> case p=0. A good first step would be to look at how the solution looks if 
> you 
> visualize it -- that is often a good start towards debugging what is going 
> wrong.
>
> Best
> W.
>
> -- 
> ---

[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&utm_source=footer>
>> .
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/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 &fe_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/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&utm_source=footer>
>> .
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/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/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  
recently about something similar 
. (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&utm_source=footer>
>> .
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/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 &p) const override;
virtual Tensor<2, dim> gradient(const Point &p) const override;
};

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

template 
Tensor<2, dim> Exact::gradient(const Point &p) 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 &p) const override;
virtual Tensor<2, dim> gradient(const Point &p) const override;
};


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


Tensor<2, dim> Exact::gradient(const Point &p) 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&M 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 deal

[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 &fe_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> &grad)
  {
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> &grad_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> &grad_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

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 &cell : 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 &cell : dof_handler.active_cell_iterators())
for (const auto &face : cell->face_iterators())
{
fe_face_values.reinit(cell, face);
const auto &q_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.

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

. 

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> &points,
  std::vector> &  values) const override;
  };

  void BC_class::vector_value_list(const std::vector> &points,
  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> &points,
  std::vector> &  values) const override;
  };

void PerfConductor::vector_value_list(const std::vector> &points,
  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> &points,
  std::vector> &  values) const override;
  };

  void Inlet::vector_value_list(const std::vector> &points,
  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 &cell : dof_handler.active_cell_iterators())
  {
cell->set_refine_flag(RefinementCase<3>::cut_z);
  }

  triangulation.execute_coarsening_and_refinement();
  
for (const auto &cell : 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_ac