Re: [deal.II] Help with interpretation of PETSc and SLEPc error message.

[Daniel Arndt]

Dear David,

for debugging your code you should start with simplifying your code as much 
as possible to see what is going wrong.
In particular, you should try with not refining your mesh and confirm that 
the assembled matrices look as they should.
Can you modify the parameters in such a way that you would expect to 
reproduce the results of step-36?

Best,
Daniel


Am Montag, 16. Mai 2016 11:28:39 UTC+2 schrieb David:
>
> Hi, Tobi:
>
> Thank you so much for your suggestions. I have followed your direction and 
> checked my code again and again for all day long, but still
> could not resolve this problem. If you (or someone else) wouldn't mind, 
> could you please take a look at my code attached below?
>
> It is a really short code, of which most parts are exactly the same with 
> step-36, except some details I changed for learning purposes (I changed 
> this problem to a simple dynamical elasticity problem on a rectangular 
> domain). I know this might cause some unconveniences, but I really 
> need somesone's help. Since nobody around me uses dealii, this is the only 
> place I can turn to.
>
> Thanks in advance for any help!
>
> Best regards,
> David.
>
>
>
> On Monday, May 16, 2016 at 1:03:30 AM UTC+8, Tobi Young wrote:
>>
>> > Could anyone please take a look, it says: 
>> > 
>> > [0]PETSC ERROR: Object is in wrong state 
>> > [0]PETSC ERROR: Matrix is missing diagonal entry 8 
>> > 
>>
>> >The above tells us that the matrix is singular (missing a diagonal 
>> entry).  
>>
> > 
>
> >I can recommend you check first, if your matrix should be singular - 
>> >ie. if you are assembling him correctly; and second, check that the 
>> >SLEPc solver you are using is capable of solving a singular matrix 
>> >(not all are). You can get some information from the solver pages of 
>> >the SLEPc manual. 
>> >Hope that helps a little. 
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Help with interpretation of PETSc and SLEPc error message.

[Daniel Arndt]

Hey David,

you can always use the 'print' method of your matrix[1] to write its 
content to a file.
Doing so both for the reference code and your code, you can simply compare 
the two files.

Best,
Daniel

[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/classSparseMatrix.html#a57ee42a70d603759e17b8a204af3c4d5

Am Donnerstag, 19. Mai 2016 15:29:46 UTC+2 schrieb David:
>
> Hi, Daniel:
>
> Thank you for your suggestions and sorry about late reply. I decided to 
> simplify the starting mesh and display the resulting matrix using 
> *MatrixOut*, so that I can check whether I went wrong in assembling the 
> matrix. But after using the following code to get the output matrix 
> data file, I don't know how to display the the matrix in gnuplot with 
> their values showing in the corresponding positions.
>
> FullMatrix 
> <http://dealii.org/8.4.1/doxygen/deal.II/classFullMatrix.html> M;
>
> ...// fill matrix M with some values
>   // now write out M:
>
> MatrixOut <http://dealii.org/8.4.1/doxygen/deal.II/classMatrixOut.html> 
> matrix_out;
>
> std::ofstream out ("M.gnuplot");
> matrix_out.build_patches 
> <http://dealii.org/8.4.1/doxygen/deal.II/classMatrixOut.html#a0d01113fd37dfa352dc237f3d05c2182>
>  
> (M, "M");
> matrix_out.write_gnuplot 
> <http://dealii.org/8.4.1/doxygen/deal.II/classDataOutInterface.html#a471d836d307c6e888ab3b5ad81fe0416>
>  
> (out);
>
>
> I tried to use 
>
> plot 'M.gnuplot' matrix
>
> but all I get is some point without values showing up. Opening this file 
> with Emacs, the values can be found, but their orderings are not 
> convenient. 
> Do you know how to display this data file so that they can show the values 
> in correct order just like in MATLAB.
>
> Best regards,
> David.
>
>
>
> On Tuesday, May 17, 2016 at 3:30:22 AM UTC+8, Daniel Arndt wrote:
>>
>> >Dear David,
>> >
>> >for debugging your code you should start with simplifying your code as 
>> much as possible to see what is going wrong.
>> >In particular, you should try with not refining your mesh and confirm 
>> that the assembled matrices look as they should.
>> >Can you modify the parameters in such a way that you would expect to 
>> reproduce the results of step-36?
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Post-processing of solution (computing stress distribution)

[Daniel Arndt]

David,

assuming that you want to access the gradients at the support points of the 
element you are using, you can simply create a
appropriate Quadrature object to initialize FEValues[1]. In this case the 
quadrature points correspond to the local DoF indices.

Apart from that, you can use DataPostprocessor[2] for postprocessing. Its 
use is illustrated in step-29, step-32, and step-33.

Best,
Daniel


[1] 
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element
[2] https://www.dealii.org/8.4.0/doxygen/deal.II/classDataPostprocessor.html

Am Sonntag, 22. Mai 2016 05:00:41 UTC+2 schrieb David:
>
> Hi, everyone:
>
> I am doing a solid mechanics problem (plane stress), and I have already 
> got the displacement solution. Now I need to use them to compute stress
> distribution with the constitutive relation, namely:
>
>
> 
> So, I also need to compute shape gradient values at the nodal points, and 
> then attach these stress solutions to the corresponding nodes when output 
> results.
>
> But the problem here is that I need to know shape gradient values at 
> triangulation nodal points. It seems that I can only use FEValues to 
> achieve this, but FEValues can only give me
> shape gradient values at quadrature points. So I tried to create a second 
> order QGaussLobatto object, so that the quadrature points are the same with 
> nodal points, and the quadrature points of index 0, 1, 2 and 3 corresponds 
> respectively to each node of one cell. However it's still not that easy to 
> control the which node this quadrature object actually refers to when I 
> want to compute the stress at a given triangullation node.
>
> Dose anyone have some suggestions or better ideas about this? Any help 
> would be greatly appreciated!
>
> Best regards,
> David.
>
>
>
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Post-processing of solution (computing stress distribution)

[Daniel Arndt]

David,

the code snippet you posted is not going to work. You can't use a FEValues 
without initializing it with a cell_iterator.
If all the information you need for your postprocessing is contained in a 
data vector based on a single DoFHandler, then DataPostprocessor is what 
you want to use.
If you take a look at the interface at the signature of 
DataPostprocessor::compute_derived_quantities_scalar[1], you see that 
this provides the access you want.
Namely, you access at the vertices of the patches used for output.

Best,
Daniel

[1] 
https://www.dealii.org/8.4.0/doxygen/deal.II/classDataPostprocessor.html#af8f60ada7b7ae5417923613d87e16e0d


Am Sonntag, 22. Mai 2016 14:45:55 UTC+2 schrieb David:
>
> I decided to use DoFTool::map_dofs_to_support_points(), so my code above 
> will be like this:
>
> DoFTools::map_dofs_to_support_points(MappingQ1, dof_handler,
>   std::vector> all_nodes);
> std::vector principle_stress;
> principle_stress.resize(dof_handler.n_dofs()); //suppose dofs = 
> no_of_nodes here.
> for (unsigned int i = 0; i != dof_handler.n_dofs(); ++i)
> {
>Point p = all_nodes[i];
>Quadrature quad(p);
>FEValues fe_values (fe, quad, |update_gradient);
>
>... ... //compute stress with fe_values.shape_grad(i, q)
>principle_stress[i] = ... ...
> }
>
> //output stress results
> DataOut data_out;
> data_out.attach_dof_handler(dof_handler);
> data_out.add_data_vector(principle_stress,"principle stress")
>
> But here I'm not sure whether the ordering of points in the vector 
> "all_nodes" getting from using DoFTool::map_dofs_to_support_points() is the 
> same with dof_indices of the dof_handler. If so, there is no problem with 
> the above code. While if not so, the ordering of nodal stress solution 
> would also be different from dof_handler, as a result, some stress is 
> attached to wrong dofs by DataOut.
>
> So which case is correct, could anyone give some explanations? 
>
> Great thanks!
>
> Best,
> David.
>
> On Sunday, May 22, 2016 at 7:43:38 PM UTC+8, David wrote:
>>
>> Hi, Daniel:
>>
>> Thank you for your help. But this seems not to be convenient. I am 
>> thinking to run my code like the following (psuedocode)
>>
>> for (i=0; i != number_of_nodes; ++i){
>>Point pi = get_point_from_global_index(i);
>>Quadrature quad(pi);
>>FEValues fe_values (fe, quad, |update_gradients);
>>
>>... ...//do some computations with fe_values.shape_grad(i, q)
>> }
>>
>> In this way, I can easily run over all the triangulation nodes without 
>> worrying about worrying about recomputation at some nodes if I instead run 
>> over each cell.
>> But the difficullty here that I can't find a function which can do the 
>> same job as "get_point_from_global_index(i)", namely give me the point with 
>> global node index i of the triangulation.
>>
>> Is there any function in dealii with this functionality? 
>>
>> Great thanks!
>>
>> Best, 
>> David
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Post-processing of solution (computing stress distribution)

[Daniel Arndt]

Hey David,

Yes, this won't work since FEValues need a cell to reinit to get mapping 
> data. I looded at the class DataPostprocessor as you mentioned, but don't 
> quite understand the parameters passed to compute_derived_quantities_vector 
> method. I don't understand what exactly is the parameters "uh", "duh" ..., 
> I know they are the computed finite element solutions, but I don't how they 
> are composed, so I don't know how to use them. 
>
Copying from the documentation for 
DataPostprocessor::compute_derived_quantities_scalar:

"[...]uh is a reference to a vector of data values at all points, duh the 
same for gradients, dduh for second derivatives and normals is a reference 
to the normal vectors. Note, that the last four references will only 
contain valid data, if the respective flags are returned by 
get_needed_update_flags, otherwise those vectors will be in an unspecified 
state.[...]"


 

>
> As I have stated in the first post, I want to compute the stress using 
> dierivatives of the solution, namely
>
>
> 
>
> so how can I use it to compute u_k,k, u_i,j and u_j,i., does these 
> quantities the same with the parameter "duh" ? Could you (or anyone else) 
> please give some explanation?
>
Have a look at the postprocessing in compute_derived_quantities_vector in 
step-32 [1]. There you see how individual components of the solution vector 
are accessed to output the respective quantities as well as
how to compute derived quantities. In particular, the strain rate is 
computed. Modifying this to include the divergence should then be straight 
forward. 
The generated output is a std::vector. Hence, you need to decide in which 
component of this vector you want to store which component of your stress.

Best,
Daniel

[1] 
https://www.dealii.org/8.4.0/doxygen/deal.II/step_32.html#BoussinesqFlowProblemoutput_result
 

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Doubt in Mesh Movement Step 42

[Daniel Arndt]

Dear Rajat,

After the mesh partitioning is done, in a 
parallel::distributed::Triangulation the local meshes are independent of 
each other. In particular, moving a vertex on one process does not change 
anything for all the other processes.
Therefore, it is necessary to move the vertices of a given face on all 
processes that have a cell that has this face as it is done in step-42.

Of course, there is as well 
parallel::distributed::communicate_locally_moved_vertices[1] to synchronize 
changes in vertex positions if you don't keep track of this yourself.

Best,
Daniel

[1] 
https://www.dealii.org/8.4.0/doxygen/deal.II/classparallel_1_1distributed_1_1Triangulation.html#a247598f1323a9f847832e60d6c840469

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Doubt in Mesh Movement Step 42

[Daniel Arndt]

Rajat,

The IndexSet given by DoF_Handler::locally_owned_dofs gives a 
non-overlapping partitioning of the dofs.
Therefore, you are safe if you just modify values at locally owned dofs of 
a distributed vector.
Note, that you normally are not allowed to modify a ghosted vector at all, 
but need a non-ghosted one for that and
assign afterwards. If you write to the non-ghosted vector, you are only 
allowed to do so at the locally owned dofs and IndexSet::is_element
is a appropriate way to do so.

The mesh movement code does something different as the vertex positions are 
purely local values. In particular, they do not belong
to a distributed Vector as explained above.

Best,
Daniel

Am Samstag, 28. Mai 2016 22:34:41 UTC+2 schrieb RAJAT ARORA:
>
> Hi Daniel,
> Thanks for your prompt reply.
>
> While in case of mesh movement, the vertices are independent. However, if 
> want to add something to the solution vector at dofs of a vertex shared 
> between processors, 
> This mesh movement code wont be able to do that, It will add that value 
> twice for the dofs at the shared vertex.
>
> What could be an ideal way to do that?
>
> Will checking that the dof is a locally owned dof by using the is_element 
> of the indexset of locally owned dofs help ?
> Will this be slow ?
>
>
> On Wednesday, May 25, 2016 at 7:19:27 AM UTC-4, Daniel Arndt wrote:
>>
>> Dear Rajat,
>>
>> After the mesh partitioning is done, in a 
>> parallel::distributed::Triangulation the local meshes are independent of 
>> each other. In particular, moving a vertex on one process does not change 
>> anything for all the other processes.
>> Therefore, it is necessary to move the vertices of a given face on all 
>> processes that have a cell that has this face as it is done in step-42.
>>
>> Of course, there is as well 
>> parallel::distributed::communicate_locally_moved_vertices[1] to synchronize 
>> changes in vertex positions if you don't keep track of this yourself.
>>
>> Best,
>> Daniel
>>
>> [1] 
>> https://www.dealii.org/8.4.0/doxygen/deal.II/classparallel_1_1distributed_1_1Triangulation.html#a247598f1323a9f847832e60d6c840469
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Periodic boundary conditions seem to be not applied

[Daniel Arndt]

Bastian,

can you provide us with a minimal example that shows your problem? What you 
are trying to do should be possible even if it might not be the right
thing to do in your situation.

Best,
Daniel

Am Mittwoch, 1. Juni 2016 02:10:36 UTC+2 schrieb Bastien Lauras:
>
> Hi,
>
> I am working on a 2D square, a incompressible neo-hookean material. I've 
> used step 44 and modified it to work on a two field formulation (*u* and 
> p).
> I am trying to implement periodic boundary conditions for the vertical 
> displacement on lateral faces (right and left).
> My left face has boundary indicator 1, and my right face has boundary 
> indicator 3.
>
> In the make_constraints function, I've added, after constraints.clear() :
>
>DoFTools::make_hanging_node_constraints (dof_handler_ref, constraints); 
> // I have local refinement in my code
>
>   std::vector DoFHandler::cell_iterator> >
>   periodicity_vector;
>
>   const unsigned int direction = 0; // I've tried 0 and 1, I didn't 
> understan what it really was
>
>   GridTools::collect_periodic_faces(dof_handler_ref, 1, 3, direction,
> periodicity_vector);
>
> const FEValuesExtractors::Scalar x_displacement(0);
> const FEValuesExtractors::Scalar y_displacement(1);
>
>  DoFTools::make_periodicity_constraints >
>   (periodicity_vector, constraints, fe.component_mask(y_displacement));
>
>
> And then I apply Dirichlet boundary conditions on face 1 and 3 on the 
> horizontal displacement (x_displacement).
> I've also added the command :
>
> DoFTools::make_hanging_node_constraints (dof_handler_ref, constraints);
>
> before making my sparsity pattern.
> My triangulation is not a parallel::distributed one. I'm working on 
> multithreads with the workstream::run of step 44.
>
>
> The problem is that my periodic boundary conditions are not enforced. The 
> make_periodicity_constraints command does not change the number of 
> constrained degrees of freedom (that I compute using a loop over the DOFs 
> and calling constraints.is_constrained(i)).
> And moreover, I can see on my outputs that there is definitely no 
> periodicity on my vertical displacement.
>
> Where can the problem come from?
>
> Thanks a lot for your help beforehand!
>
> Bastien Lauras 
>  
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Periodic boundary conditions seem to be not applied

[Daniel Arndt]

It is quite difficult to say what is going wrong without seeing the code. 
You probably want to check first if you really want to impose that kind of 
boundary conditions and follow Jean-Paul's advice.

If the problem persists, it would be really helpful if you can give us a 
minimal compiling code that reproduces the problem.

Best,
Daniel

Am Mittwoch, 1. Juni 2016 18:05:32 UTC+2 schrieb Bastien Lauras:
>
> Hi Daniel,
>
> So, I'm modelling a unit square, with an incompressible neo-hookean 
> material. The formulation is (for the density of energy) : c_1 * (I_1 - 3) 
> - p (I_3 - 1) where c_1 = mu / 2, p is a lagrangian parameter, and I_1 and 
> I_3 are the first and third invariant of C = F^T * F.
> For me to see that my periodic boundary conditions are not applied, I put 
> a non-uniform mu aver the material (graded material). In the example 
> attached, I blocked the horizontal displacement (x_displacement) over left 
> and right faces, and blocked vertical displacement (y_displacement) over 
> the bottom face. And I (try to) apply a periodic boundary condition for the 
> vertical displacement over left and right faces.
> Attached is the output given to me after the first timestep (all timesteps 
> are the same since I'm not applying any displacement nor force constraint).
>
> We can obviously see that vertical displacement is not the same on the 
> left and right faces. Any idea of where it could come from?
>
> Many thanks!
> Best,
>
> Bastien
>
>
> On Wednesday, June 1, 2016 at 5:43:48 AM UTC-5, Daniel Arndt wrote:
>>
>> Bastian,
>>
>> can you provide us with a minimal example that shows your problem? What 
>> you are trying to do should be possible even if it might not be the right
>> thing to do in your situation.
>>
>> Best,
>> Daniel
>>
>> Am Mittwoch, 1. Juni 2016 02:10:36 UTC+2 schrieb Bastien Lauras:
>>>
>>> Hi,
>>>
>>> I am working on a 2D square, a incompressible neo-hookean material. I've 
>>> used step 44 and modified it to work on a two field formulation (*u* and 
>>> p).
>>> I am trying to implement periodic boundary conditions for the vertical 
>>> displacement on lateral faces (right and left).
>>> My left face has boundary indicator 1, and my right face has boundary 
>>> indicator 3.
>>>
>>> In the make_constraints function, I've added, after constraints.clear() :
>>>
>>>DoFTools::make_hanging_node_constraints (dof_handler_ref, 
>>> constraints); // I have local refinement in my code
>>>
>>>   std::vector>> DoFHandler::cell_iterator> >
>>>   periodicity_vector;
>>>
>>>   const unsigned int direction = 0; // I've tried 0 and 1, I didn't 
>>> understan what it really was
>>>
>>>   GridTools::collect_periodic_faces(dof_handler_ref, 1, 3, direction,
>>> periodicity_vector);
>>>
>>> const FEValuesExtractors::Scalar x_displacement(0);
>>> const FEValuesExtractors::Scalar y_displacement(1);
>>>
>>>  DoFTools::make_periodicity_constraints >
>>>   (periodicity_vector, constraints, 
>>> fe.component_mask(y_displacement));
>>>
>>>
>>> And then I apply Dirichlet boundary conditions on face 1 and 3 on the 
>>> horizontal displacement (x_displacement).
>>> I've also added the command :
>>>
>>> DoFTools::make_hanging_node_constraints (dof_handler_ref, constraints);
>>>
>>> before making my sparsity pattern.
>>> My triangulation is not a parallel::distributed one. I'm working on 
>>> multithreads with the workstream::run of step 44.
>>>
>>>
>>> The problem is that my periodic boundary conditions are not enforced. 
>>> The make_periodicity_constraints command does not change the number of 
>>> constrained degrees of freedom (that I compute using a loop over the DOFs 
>>> and calling constraints.is_constrained(i)).
>>> And moreover, I can see on my outputs that there is definitely no 
>>> periodicity on my vertical displacement.
>>>
>>> Where can the problem come from?
>>>
>>> Thanks a lot for your help beforehand!
>>>
>>> Bastien Lauras 
>>>  
>>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Thermoelastic Problem

[Daniel Arndt]

Hamed,

There shouldn't be a problem to do what you following step-22. In 
particular, you can modify step-18::get_strain as you proposed.
Did you face any difficulties so far?

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Thermoelastic Problem

[Daniel Arndt]

You're right, Hamed. There is FEValuesViews::Vector< dim, spacedim 
>::symmetric_gradient [1] 
and FEValuesViews::Vector< dim, spacedim >::gradient [2] instead. The first 
should be exactly what you are looking for.

Code that I have written normally uses multiple DoFHandlers if I solve for 
multiple equations separately.

Best,
Daniel

[1] 
https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html#a4e5dfbb49d284a368acac6ef5dd4b2ef
 
<https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html>
[2] 
https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html#ad7b4df64147989229f566eff541ddebc
 
<https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html>






Am Freitag, 3. Juni 2016 01:33:34 UTC+2 schrieb Hamed Babaei:
Daniel,

Since in the FEValuesViews::Vector< dim, spacedim > Class there isn't any 
shape_grad_component function, I am not sure that I can update the 
get_strain as I mentioned.
I was wondering if you have already written a code for a thermoelastic 
problem in dealii, solving heat equation and elastic equation separately.

Best,
Hamed

On Thursday, June 2, 2016 at 3:31:48 PM UTC-5, Daniel Arndt wrote:
Hamed,


There shouldn't be a problem to do what you following step-22. In 
particular, you can modify step-18::get_strain as you proposed.
Did you face any difficulties so far?


Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Periodic boundary conditions seem to be not applied

[Daniel Arndt]

You are right, Bastien. The function is not implemented in 1d.
However, you can just ask for the boundary dofs using 
DoFTools::extract_boundary_dofs [1].
This gives you two DoFs and you can just set the constraint using
constraint_matrix.add_line(dof_1)
constraint_matrix.add_entry(dof_1,dof_2,-1.)

Best,
Daniel

[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/namespaceDoFTools.html#ad3067f335a97de429176178689222f3a

Am Freitag, 3. Juni 2016 01:42:58 UTC+2 schrieb Bastien Lauras:
>
> Hi,
>
> Ok, after reviewing it, I am going to impose a constraint of the form : 
> u(0,y) = u(1,y) + lambda and v(0,y) = v(1,y) (lambda is the total 
> compression of my square). If I'm not mistaken, that is what Jean-Paul was 
> advicing.
> I have two example, the first one is the 2D incompressible neo-hookean 
> square, where I would like to implement the constraint just above. I don't 
> think it will solve the problem I had but I'll try and will get back to you.
>
> The second example is a 1D-beam on a linear elastic foundation. The only 
> degree of freedom is the vertical displacement. So everything is in one 
> dimension. I want to impose the constraint y(0) = y(1) (only vertical 
> displacement y allowed).
> Here is the declaration of my FESystem : 
> fe(FE_Q(parameters.poly_degree), dim)
>
> This needs to be changed to Hermitian cubic elements, but this is another 
> topic. My triangulation is :
> Triangulation   triangulation;
> with dim = 1. I then construct it, refine it, and mark the "faces". But 
> what are the faces in 1D?
> GridGenerator::hyper_cube(triangulation, 0.0, 1.0);
> triangulation.refine_global(std::max (1U, 
> parameters.global_refinement));
>
> const double length = GridTools::volume(triangulation);
> std::cout << "Grid:\n\t Length: " << length << std::endl;
>
> // We mark the sides of the beam in order to apply the boundary 
> conditions after
> typename Triangulation::active_cell_iterator cell =
>   triangulation.begin_active(), endc = triangulation.end();
> for (; cell != endc; ++cell)
>   for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f)
> if (cell->face(f)->at_boundary())
>   {
> const Point face_center = cell->face(f)->center();
> if (face_center[0] == 0) //Left
>   cell->face(f)->set_boundary_id (1);
> else if (face_center[0] == 1) //Right
>   cell->face(f)->set_boundary_id (2);
> else  //No way
>   cell->face(f)->set_boundary_id (3);
> // Is this working in 1D ?
>   }
>
> Then I make my sparsity pattern and all the like, and begin my timesteps, 
> and so on, my Newton-Raphson iterations. The function make_constraints is 
> as follows :
>  template 
>   void Solid::make_constraints(const int &it_nr)
>   {
> // Since the constraints are different at different Newton iterations, 
> we
> // need to clear the constraints matrix and completely rebuild
> // it. However, after the first iteration, the constraints remain the 
> same
> // and we can simply skip the rebuilding step if we do not clear it.
> if (it_nr > 1)
>   return;
>
> constraints.clear();
> std::vector DoFHandler::cell_iterator> >
> periodicity_vector;
>
> const int direction = 0;
> GridTools::collect_periodic_faces(dof_handler_ref, 1, 2, direction,
> periodicity_vector);
>
> const FEValuesExtractors::Scalar y_displacement(0);
>
> //DoFTools::make_periodicity_constraints >
> //  (periodicity_vector, constraints/*, 
> fe.component_mask(y_displacement)*/);
> DoFTools::make_periodicity_constraints (dof_handler_ref,1,2, 
> direction,
> constraints,
> fe.component_mask(y_displacement));
> // I tried both implementations of this function, but noone is working.
> constraints.close();
>   }
>
> When I try to compile my file, I get the following error :
> error: undefined reference to 'void 
> dealii::DoFTools::make_periodicity_constraints 
> >(dealii::DoFHandler<1, 1> const&, unsigned char, unsigned char, int, 
> dealii::ConstraintMatrix&, dealii::ComponentMask const&)'
>
> I think the problem comes from the fact I am working in 1D, doesn't it?
>
> Many thanks,
>
> Bastien
>
> On Wednesday, June 1, 2016 at 11:55:31 AM UTC-5, Daniel Arndt wrote:
>>
>> It is quite difficult to say what is going wrong without seeing the code. 
>> You probably want to check first if you really want to impose that kind 
>> of boundary condition

[deal.II] Re: Thermoelastic Problem

[Daniel Arndt]

Hamed,

assuming that both DoFHandlers are based on the same triangulation, you can 
do something like the following:

FEValues fe_stress_values(...);
stress_cell = dof_handler_stress.begin_active();
temperature_cell = dof_handler_temperature.begin_active();
for (; stress_cell != dof_handler_stress.end(); ++stress_cell, 
++temperature_cell)
{
  fe_stress_values.reinit(stress_cell);
  fe_stress_values.get_function_values(solution_stress, local_stress);
  //do all the temperature related stuff and use local_stress
}

Instead of having two iterators separated, you can also use 
SynchronousIterators [1]. They arec partly used in step-35 for that purpose.

Best,
Daniel

[1] https://dealii.org/8.4.1/doxygen/deal.II/structSynchronousIterators.html

Am Dienstag, 7. Juni 2016 20:30:27 UTC+2 schrieb Hamed Babaei:
>
> Hello Daniel,
>
> Since although we can solve two elastic and heat equations separately, we 
> still need to enter stress computed from elastic equation into the heat 
> equation in every time step and keep updating stress field.
> I was wondering when we have two different DoFHandlers for each of the 
> equations, how to use, for instance, stress field stored on a quadrature 
> point related to elastic problem DofHandler, in order to solve 
> stress-induced heat equation problem. In other words, I don't know how to 
> exchange information between two equation.
>
> Thanks,
> Hamed
>
> On Friday, June 3, 2016 at 3:40:28 AM UTC-5, Daniel Arndt wrote:
>>
>> You're right, Hamed. There is FEValuesViews::Vector< dim, spacedim 
>> >::symmetric_gradient [1] 
>> and FEValuesViews::Vector< dim, spacedim >::gradient [2] instead. The 
>> first should be exactly what you are looking for.
>>
>> Code that I have written normally uses multiple DoFHandlers if I solve 
>> for multiple equations separately.
>>
>> Best,
>> Daniel
>>
>> [1] 
>> https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html#a4e5dfbb49d284a368acac6ef5dd4b2ef
>>  
>> <https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html>
>> [2] 
>> https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html#ad7b4df64147989229f566eff541ddebc
>>  
>> <https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html>
>>
>>
>>
>>
>>
>>
>> Am Freitag, 3. Juni 2016 01:33:34 UTC+2 schrieb Hamed Babaei:
>> Daniel,
>>
>> Since in the FEValuesViews::Vector< dim, spacedim > Class there isn't any 
>> shape_grad_component function, I am not sure that I can update the 
>> get_strain as I mentioned.
>> I was wondering if you have already written a code for a thermoelastic 
>> problem in dealii, solving heat equation and elastic equation separately.
>>
>> Best,
>> Hamed
>>
>> On Thursday, June 2, 2016 at 3:31:48 PM UTC-5, Daniel Arndt wrote:
>> Hamed,
>>
>>
>> There shouldn't be a problem to do what you following step-22. In 
>> particular, you can modify step-18::get_strain as you proposed.
>> Did you face any difficulties so far?
>>
>>
>> Best,
>> Daniel
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] new shape function definition

[Daniel Arndt]

Bastien,

as there already is the Polynomial class HermiteInterpolation you should 
create a new FiniteElement based on a 
TensorProductPolynomial. A good way might be to copy 
FE_Q_Base and modify it accordingly.
Of course you need to think about how to interpolate values using this 
element, how to use constraints, etc.
You don't need to implement everything (like restriction and interpolation 
matrices) immediately.
If you are interested, just start and ask when you're stuck.

Best,
Daniel

Am Montag, 6. Juni 2016 20:04:04 UTC+2 schrieb Bastien Lauras:
>
> Dear Guido and dear all,
>
> I am trying to implement Hermitian Cubic elements on a 1D beam. I am 
> working on the vertical displacement y of the beam (only displacement 
> allowed, the beam is inextensible). I need a continuous first derivative of 
> the displacement, and this can easily be done with Hermitian Cubic elements.
> But I don't understand well the Polynomial module. I've seen that Hermite 
> Interpolation 
> 
>  exists 
> in deal.II. Though, I don't understand how I can define a polynomial space, 
> so that FE_Q would work with my Hermitian elements and not Lagrangian 
> elements.
>
> Thanks a lot for your help!
>
> Best,
>
> Bastien
>
> On Thursday, September 27, 2012 at 7:33:20 AM UTC-5, Guido Kanschat wrote:
>>
>> Dear Nicola, 
>>
>> typically, we define finite element shape function spaces in two steps: 
>>
>> 1. Define a polynomial space. Here, you can find examples in the 
>> module on Polynomials, in particular classes like 
>> TensorProductPolynomials or PolynomialsBDM. Your task would 
>> essentially be adding an additional constant function to 
>> TensorProductPolynomials. 
>>
>> 2. Create a finite element using your new polynomial space. FE_Q would 
>> be your model here.  Unfortunately, FE_Q has a lot of complicated 
>> hidden logic like the renumbering of polynomials from tensor product 
>> structure to internal deal.II enumeration. You would have to make 
>> sure, that your new polynomial at the end of the list stays at the 
>> end, since it is located in the interior of the cell. 
>>
>> I can help with the implementation if you get it started. 
>>
>> Best, 
>> Guido 
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Periodic boundary conditions seem to be not applied

[Daniel Arndt]

u may not be colouring any of your boundaries. This should be a 
>>tolerant check, something like "if (std::abs(face_center[1] -1.0) <= 
>>1e-6)" to take into account computational errors.
>>3. The offset vector should probably be all zero since the vertices 
>>align perfectly on either side of a hypercube. With an offset of 0.1 and 
>> a 
>>regular global refinement, if I'm not mistaken its impossible to get two 
>>vertices that are 0.1 units offset from one another.  From the 
>>documentation of the collect_periodic_faces function:
>>
>> The offset is a vector tangential to the faces that is added to the 
>>> location of vertices of the 'first' boundary when attempting to match them 
>>> to the corresponding vertices of the 'second' boundary
>>
>>
>> I hope that this helps,
>> J-P
>>
>> On Tuesday, June 7, 2016 at 6:47:51 PM UTC+2, Bastien Lauras wrote:
>>>
>>> Hi Daniel.
>>>
>>> Thank you for your answer, it worked well to enforce my periodic 
>>> boundary conditions on the 1D beam.
>>> Thus, I'm back to the 2D neo-hookean case. Here is the code of my 
>>> *make_grid()* function :
>>>
>>>   template 
>>>   void Solid::make_grid()
>>>   {
>>> GridGenerator::hyper_cube(triangulation, 0.0, 1.0);
>>>
>>> // We will refine the grid in five steps towards the inner circle of
>>> // the domain. See step 1 for details.
>>> for (unsigned int step=0; step>> ++step)
>>>   {
>>> typename Triangulation::active_cell_iterator
>>>  cell_ref = triangulation.begin_active(),
>>>  endc = triangulation.end();
>>> for (; cell_ref!=endc; ++cell_ref)
>>>  {
>>>for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f)
>>>  if (cell_ref->face(f)->at_boundary())
>>> {
>>>  const Point face_center = cell_ref->face(f)->center();
>>>  if (face_center[1] == 1) //Top
>>>cell_ref->set_refine_flag ();
>>> }
>>>  }
>>>
>>> // Now that we have marked all the cells that we want refined, we let
>>> // the triangulation actually do this refinement.
>>> triangulation.execute_coarsening_and_refinement ();
>>> triangulation.clear_user_flags();
>>>   }
>>>
>>> // We refine our mesh globally, at least once, to not too have a poor
>>> // refinement at the bottom of our square (like two cells)
>>> triangulation.refine_global(std::max (1U, 
>>> parameters.global_refinement));
>>>
>>> // We mark the surfaces in order to apply the boundary conditions 
>>> after
>>> typename Triangulation::active_cell_iterator cell =
>>>   triangulation.begin_active(), endc = triangulation.end();
>>> for (; cell != endc; ++cell)
>>>   for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; 
>>> ++f)
>>> if (cell->face(f)->at_boundary())
>>>  {
>>>const Point face_center = cell->face(f)->center();
>>>if (face_center[1] == 0)  //Bottom
>>>  cell->face(f)->set_boundary_id (0);
>>>else if (face_center[1] == 1) //Top
>>>  cell->face(f)->set_boundary_id (2);
>>>else if (face_center[0] == 0) //Left
>>>  cell->face(f)->set_boundary_id (1);
>>>else if (face_center[0] == 1) //Right
>>>  cell->face(f)->set_boundary_id (3);
>>>else  //No way
>>>  cell->face(f)->set_boundary_id (4);
>>>  }
>>>
>>>   std::vector>> Triangulation::cell_iterator> >
>>>   periodicity_vector;
>>>
>>>   const unsigned int direction = 0;
>>>
>>>   FullMatrix rotation_matrix(dim);
>>>   rotation_matrix[0][0]=1.;
>>>   rotation_matrix[1][1]=1.;
>>>
>>>   Tensor<1, dim> offset;
>>>   offset[0]=0.1;
>>>
>>>   GridTools::collect_periodic_faces(triangulation, 1, 3, direction,
>>> periodicity_vector, offset, 
>>> rotation_matrix);
>>>
>>>   std::cout << "\nSize of periodicity_vector : " << 
>>> periodicity_vector.size();
>>>   }
>>>
>>> The last output gives me "Size of periodicity_vector : 0". I don't 
>>> understant why. 

[deal.II] Re: Using the solution nodal values at previous time steps

[Daniel Arndt]

Mohammad,

first of all 

const double old_s = old_solution_values[q_point](1);

const double old_s = old_solution_values0[q_point](1);

const double old_s = old_solution_values1[q_point](1);

looks weird.  Are you redefining old_s all the time?


Apart from that it looks so far good to me. Can you confirm that 
old_solution_values1 is the same as old_solution_values in the previous 
time step?

Can you make any sense of what is stored in old_solution_values, 
old_solution_values1, 
old_solution_values2?


Best,

Daniel

Am Dienstag, 7. Juni 2016 11:29:45 UTC+2 schrieb Mohammad Sabawi:
>
> Dear all,
>
>  
>
> I am solving a 3x3 block system arising form DG time discretisation of the 
> semilinear parabolic equations in the form u_t - \Delta u = f(u) for the 
> block solution vector U which consists from three nodal values blocks U0, 
> U1 and U2. During the solution process I need just the solution values of 
> the second block U1 at the previous time steps (t_{n-1}, t_{n-2} and 
> t_{n-3})
>
> i.e. old_solution.block(1), old_old_solution.block(1) and 
> old_old_old_solution.block(1) for this task I used
>
>  
>
> std::vector old_solution_values (n_q_points, 
> Vector(3));
>
> std::vector old_solution_values0 (n_q_points, 
> Vector(3));
>
> std::vector old_solution_values1 (n_q_points, 
> Vector(3));
>
>  
>
> fe_values.get_function_values (old_solution, old_solution_values);
>
> fe_values.get_function_values (old_old_solution, old_solution_values0);
>
> fe_values.get_function_values (old_old_old_solution, old_solution_values1);
>
>  
>
>  
>
> In the solution process in the right hand side assembly I need just the 
> nodal values of second block U1 at the previous time steps (t_{n-1}, 
> t_{n-2} and t_{n-3}) i.e.  old_solution.block(1)
>
>  
>
>  For this reason, I used
>
>  
>
> const double old_s = old_solution_values[q_point](1);
>
> const double old_s = old_solution_values0[q_point](1);
>
> const double old_s = old_solution_values1[q_point](1);
>
>  
>
>  
>
> but I did not get required results and even when I changed 
> old_solution_values[q_point](1) to old_solution_values[q_point](0) nothing 
> has changed (just to check if there is any difference) and even when I 
> tried the wrong thing old_solution.block(1) nothing has changed and I got 
> the same results as with previous other cases.
>
>  
>
> I am confused about this case; did I do some thing wrong in defining the 
> nodal values of the second block of the solution vector at the previous 
> time steps? Any feedbacks and hints are appreciated.
>
>  
>
> Best regards
>
>  
>
> Mohammad Sabawi
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Periodic boundary conditions seem to be not applied

[Daniel Arndt]

Bastien,

Unfortunately, you won't be able to achieve what you want directly with 
make_periodicity_constraints.
The offset vector would help you with something like u(0,y)=u(1,y+lambda).
For what you want to do you would first create the constraints for 
u(0,y)=u(1,y) and modify the corresponding inhomogeneity afterwards.
To figure out the correct sign you would need to identify on which face the 
constrained DoF is located. For doing that 
DoFTools::map_dofs_to_support_points [1] might be useful.

Best,
Daniel

[1] 
https://www.dealii.org/8.4.0/doxygen/deal.II/namespaceDoFTools.html#a10af6ba7adc43214e7b9f6815588f2aa


Am Mittwoch, 8. Juni 2016 20:33:02 UTC+2 schrieb Bastien Lauras:
>
> Hi,
>
> Thanks a lot, it works well now! I don't know if I misread step-45, but I 
> don't think it is written in it that boundary_ids should be setup before 
> any refinement. And if you take something like step-44, boundry_ids are 
> setup after refining the mesh.
> So maybe it would be wise to detail it in step-45! Just my humble 
> thinking, as a beginner in deal.II.
>
> So, I can now enforce periodic boundary conditions like : u( 0, y ) = u( 
> 1, y ) and v ( 0, y ) = v( 1, y ).
> However, I would like something like : u( 0, y ) = u( 1, y ) + *lambda* 
> and v ( 0, y ) = v( 1, y ). I thought I could do that with the offset 
> vector, but this does not seem to work. Should I rather use this function 
> <https://dealii.org/developer/doxygen/deal.II/namespaceDoFTools.html#a4562db8677b300aeec1e3dd50ad01562>
>  and 
> use the two last arguments to do it?
>
> Thank you very much!
>
> Best,
>
> Bastien
>
> On Tuesday, June 7, 2016 at 11:41:53 PM UTC-5, Jean-Paul Pelteret wrote:
>>
>> I agree with Daniel. So in this case, you probably just need to change 
>> your order of operations and set the boundary_id's before doing any 
>> refinement.
>>
>> On Wednesday, June 8, 2016 at 2:07:09 AM UTC+2, Daniel Arndt wrote:
>>>
>>> Bastien,
>>>
>>> GridTools::collect_periodic_faces stores only the periodic faces on the 
>>> coarsest level.
>>> This information is used in DoFTools::make_periodicity_constraints to 
>>> create the corresponding constraints
>>> on the active set of DoFs.
>>> My best guess would be that you don't set the boundary_ids on the coarse 
>>> level correctly.
>>>
>>> Best,
>>> Daniel
>>>
>>> Am Mittwoch, 8. Juni 2016 00:47:56 UTC+2 schrieb Bastien Lauras:
>>>>
>>>> Hi Jean-Paul,
>>>>
>>>> Concerning the problem you've pointed out :
>>>> 1. Ok, I'll take it into account. From now on, let's try to solve this 
>>>> problem without local refinement. I do not do it anymore here.
>>>> 2. Thanks for the advice, I've changed my checks.
>>>> 3. I understand. However, I am almost sure that an offest vector non 
>>>> zero in the direction of the vector pointed by direction should not change 
>>>> our result. Well, I tried to let this offset be zero. It did not change 
>>>> anything.
>>>>
>>>> Bad news then! However, I've been working all day on this 
>>>> collect_periodic_faces function. I copied it from the git hub repository 
>>>> and played with it.
>>>>
>>>>   template
>>>>  void
>>>>  collect_periodic_faces
>>>>  (const MeshType&mesh,
>>>>   const types::boundary_id   b_id1,
>>>>   const types::boundary_id   b_id2,
>>>>   const int  direction,
>>>>   std::vector > &
>>>> matched_pairs,
>>>>   const Tensor<1,MeshType::space_dimension> &offset,
>>>>   const FullMatrix  &matrix)
>>>>  {
>>>>static const int dim = MeshType::dimension;
>>>>static const int space_dim = MeshType::space_dimension;
>>>>(void)dim;
>>>>(void)space_dim;
>>>>Assert (0<=direction && direction>>>ExcIndexRange (direction, 0, space_dim));
>>>>
>>>> // Loop over all cells on the highest level and collect all 
>>>> boundary
>>>>// faces belonging to b_id1 and b_id2:
>>>>
>>>> std::set 
>>>> > pairs1;
>>>>std::set > 
>>>> pairs2;
>>>>
>>>> for (typename MeshType::cell_iterator cell = mesh.begin(0);
>>>> cell != mesh.end(0); ++cel

[deal.II] Re: computing the buoyancy force

[Daniel Arndt]

Pankaj,

Hello,
> I am trying to write a code to calculate buoyancy force which is the 
> product of temperature and gravity. In my case temperature is available.
>
>
>
>   TrilinosWrappers::MPI::Vector temperature ;//given
>   const dealii::Point gravity = -( (dim == 2) ? (dealii::Point (
> 0,1)) :
> (dealii::Point (0,0,1)) );
>
>   TrilinosWrappers::MPI::Vector& Buoyancy_force = constant * temperature * 
> gravity ; //This is what i want to do
>
If you want to calculate your buyoancy as TrilinosWrappers::MPI::Vector, 
you need a corresponding DoFHandler.

 

>
>
>
> //below is the idea that i am not sure of
>
>typename DoFHandler::active_cell_iterator
>  cell   = fe_velocity_space.dh_begin_active(),
>  endc   = fe_velocity_space.dh_end(),
>  cell_t = fe_temp_space.dh_begin_active();
>for (; cell != endc; ++cell, ++cell_t) {
> if (cell->is_locally_owned() == false)
>   continue;
>
> b = 0;
>
> fe_velocity_values.reinit(cell);
> fe_temp_values.reinit(cell_t);
> const FEValuesViews::Vector& fe_vector_values = 
> fe_velocity_values[FEValuesExtractors::Vector(0)];
> fe_temp_values.get_function_values(temperature, temp_cell);
>
> double Beta = 0.04;
>
> for (unsigned int q=0; q const double& temp_q = temp_cell[q];
> const Tensor<1, dim> force = - Beta *(temp_q)*fe_values.
> JxW(q)*gravity ; 
> //Need help here
>
>
>
> cell->get_dof_indices(dof_indices);
>
> fe_velocity_space.get_constraints().distribute_local_to_global
> (b, dof_indices, buoyancy_force);
>
>}
>
I guess that in all places, in which you are using fe_velocity, you mean to 
use a FiniteElement for the buoyancy force? The variable b should be force, 
correct?
What you are doingso far is to compute the buoyancy force in each 
quadrature point. If your FiniteElement is interpolating in its 
support_points (as FE_Q), then
you want to use as quadrature points in the FEValues objects 
fe_velocity_values and fe_temp_values the unit_support_points of the 
FiniteElement [1].
The number of quadrature points you get after reinitializing the FEValues 
objects equals the number of DoFs in each component multiplied by the 
number of components and contains the mapped support points. Then you can 
simply assign the values you computed above. In particular, you can do 
something like:

Quadrature q(fe.get_unit_support_points()); 
FEValues fe_values (..., q, update_q_points);
...
for (unsigned int q=0; q force = - Beta *(temp_q)*fe_values.JxW(q)*gravity ;
const unsigned int component = 
FiniteElement::system_to_component_index(q).first;
b(q)=force(component)
}
...  

Note that in this snippet, the temperature in the same quadrature point is 
calculated multiply.

Best,
Daniel

[1] 
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Read serial vector into a parallel vector

[Daniel Arndt]

Praveen,

I solve a PDE in serial and save the solution to file.
>
Do you mean that your triangulation is not a 
parallel::distributed::Triangulation or that you are simply running with 
one process?
 

> Is there a way to now read this solution into a 
> TrilinosWrappers::MPI::Vector object ?
>
If you already save the solution from a TrilinosWrappers::MPI::Vector 
object, then parallel::distributed::SolutionTransfer should work for you as 
well.
Note that you are allowed to choose the number of processes when reading 
different from the number of processes you used for writing.
If your original vector is not a TrilinosWrappers::MPI::Vector object, how 
did you write the vector to a file?

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Read serial vector into a parallel vector

[Daniel Arndt]

Praveen,

Am Montag, 20. Juni 2016 18:05:12 UTC+2 schrieb Praveen C:
>
> Hello Daniel
>
> I have a normal Triangulation and I solve on this using a Vector. 
> I could save this to file calling Vector::print.
>
> Is it now possible to read this into a TrilinosWrappers::MPI::Vector ?
>
No, there is no such function. A Vector does only make sense in connection 
with a DoFHandler.
Therefore the question is what you want to achieve in the first place.
If you want to write a TrilinosWrappers::Vector to a file and initialize a 
TrilinosWrappers::MPI::Vector with it aftwerwards,
why don't you use a TrilinosWrappers::MPI::Vector from the beginning?

Something you can also try is to associate the DoFs in the serial vector 
with the ones in the parallel vector (e.g. by support_points and component) 
and assign
the values accordingly.
Note that there is also the possibility to initialize a (parallel) 
TrilinosWrappers::MPI::Vector from a (serial) TrilinosWrappers::Vector 
directly [1].

Best,
Daniel

[1] 
https://dealii.org/8.4.1/doxygen/deal.II/classTrilinosWrappers_1_1MPI_1_1Vector.html#a814279778da76eb9a17eedce154a61f7

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Extract boundary DoFs using DoFTools::extract_boundary_dofs() and Indexset for parallel distributed triangulation

[Daniel Arndt]

Dear ?,

dealii::IndexSet::ElementIterator::operator*() returns the stored 
global_dof_index [1]. 
Therefore, all you have to do is to dereference index:

dealii::types::global_dof_index gdi = *index

Best,
Daniel

[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/classIndexSet_1_1ElementIterator.html#a6218ab9e8e76699fc998834bf8e006ee

Am Donnerstag, 30. Juni 2016 15:12:00 UTC+2 schrieb dealii...@gmail.com:
>
> Hi
>
> I am trying to extract all degrees of freedom which are at the boundary 
> and belong to specified components using DoFTools::extract_boundary_dofs() 
> function and IndexSet for a parallel distributed triangulation. Now, how I 
> can iterate over the IndexSet and how I can access to the Dof index in it?
>
> Something like this:
>
> dealii::IndexSet::ElementIterator index = in_set.begin(), endind 
> = in_set.end();
>
>  for (; index!=endind; ++index)
>  {
>  how to access to dof index here?
>  }
>
> Thank you
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Extract boundary DoFs using DoFTools::extract_boundary_dofs() and Indexset for parallel distributed triangulation

[Daniel Arndt]

Pasha,

First of all: The sum of values in the right-hand side vector only 
corresponds to a sum of values of the represented right-hand side function, 
if you use interpolating Finite Elements (such as FE_Q).
Furthermore, you have to be sure that you are really interested in the 
values and not in the integral over the boundary or something similar.

Apart from the last issue you mentioned, your approach seems to be correct 
but is not really efficient.
For summing up all the values in the right-hand side vector corresponding 
to degrees of freedom on the boundary, you would just want to use your 
inner part with small modifications:

dealii::IndexSet::ElementIterator index = in_set.begin(),
dealii::IndexSet::ElementIterator endind = in_set.end();

dealii::IndexSet locally_owned_dofs = 
dof_handler.locally_owned_dofs();

for (; index!=endind; ++index)
{
dealii::types::global_dof_index gdi = *index;
// check if this DoF is locally owned
if (locally_owned_dofs.is_element(gdi))
  local_load -= saved_residual(gdi);
}

Note that checking whether the DoF is locally owned, you circumvent the 
problem of taking the same DoF in to account multiple times on different 
MPI processes.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Question about Dealii

[Daniel Arndt]

Hamed,

Have a look at how temporary non-ghosted vectors are used in step-40.
Similar to VectorTools::interpolate
the solution vector in SolverCG::solve must not be ghosted.

Best,
Daniel

Am 07/01/2016 um 01:15 AM schrieb Hamed Babaei:
> Hi Daniel,
>
> I hope you are doing well.
> I got an error in my parallelized code at the run part where I want to
> fill my solution vector with its initial values:
>
>  VectorTools::interpolate (dof_handler,
>InitialValues(),
>solution);
> The error says: "You are trying an operation on a vector that is only
> allowed if the vector has no ghost elements, but the vector you are
> operating on does have ghost elements. Specifically, vectors with
> ghost elements are read-only and cannot appear in operations that
> write into these vectors."
>
> In fact, the solution vector is a ghosted vector and I red one of your
> suggestions in forum about using temporary non-ghosted vectors but I
> don't know how to do so.
> It would be appreciated if you could give me an example or point me to
> the Dealii documentation where I can find the remedy.
>
> Thank you for your help and kindness.
>
> Best Regards, 

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range

[Daniel Arndt]

Ehsan,

All I can say: After switching the order of arguments in SparseMatrix::add, 
your code runs for me with a recent developer version and Trilinos at least.

Best,
Daniel

Am Montag, 4. Juli 2016 18:59:05 UTC+2 schrieb Ehsan Esfahani:
>
> Dear Professor Bangerth,
>
> Thanks for your response. Yes, I did. As I mentioned, I got a backtrace in 
> the debugger (eclipse) and I find out that the problem is in the line I 
> have mentioned but I couldn't find out what's the problem in that line of 
> the code which causes segmentation violation.
>
> Best,
> Ehsan
>
>
> On Sunday, July 3, 2016 at 4:32:16 PM UTC-5, bangerth wrote:
>>
>> On 07/03/2016 03:50 PM, Ehsan Esfahani wrote: 
>> > Dear All, 
>> > 
>> > Greetings. I changed step-25 (minor changes) in order to solve my 
>> problem. Now 
>> > I want to change this code for parallel computation based on the method 
>> > mentioned in step-40. I got several errors and solved them one by one 
>> but the 
>> > following error: 
>> > 
>> > /Number of active cells: 1024 
>> > //   Total number of cells: 1365 
>> > //{[0,4224]}// 
>> > //Time step #1; advancing to t = 0.1. 
>>  > [...] 
>> > //[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation 
>> Violation, 
>> > probably memory access out of range 
>>  > [...] 
>> > 
>> > / 
>> > / 
>> > Eclipse gives me a backtrace in the following line of the code: 
>> > /solver.solve (system_matrix, 
>> completely_distributed_solution_update, 
>> > system_rhs,/ 
>> > /  preconditioner);/ 
>> > I have no idea why I got this error. The code is running properly for 
>> /fe(1) 
>> > /and /n_global_refinements (4)/ but when I change them to /fe(2)/ and 
>> > n_global_refinments (4) I got that error related to /Segmentation 
>> Violation. 
>> > /Do you know what's going on? Also, I have attached the code here . 
>> Thanks in 
>> > advance for your help. 
>>
>> Ehsan, 
>>
>> segmentation violations (SEGV) are typically easy to debug because you 
>> can get 
>> a backtrave in the debugger of the exact place where it happens, and you 
>> can 
>> then look at the local variables to see why this may have happened. Have 
>> you 
>> tried to run the program in a debugger and see what is going on? 
>>
>> Best 
>>   W. 
>>
>> -- 
>>  
>> Wolfgang Bangerth   email:bang...@math.tamu.edu 
>>  www: http://www.math.tamu.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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: problems with installation ; Unknown CMake command "DEAL_II_ADD_LIBRARY"

[Daniel Arndt]

Pablo,

Did you see https://groups.google.com/forum/#!topic/dealii/rTx7Fea65dM ?
Can you try without Opencascade and a clean build directory? Do you get 
this error also with the latest release tarball?

Best,
Daniel

Am Montag, 4. Juli 2016 10:15:34 UTC+2 schrieb Pablo Perez del Castillo:
>
> Hello Timo,
>
> I got same Error. No idea.
>
> Pablo
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Geometry and boundary conditions

[Daniel Arndt]

Benhour,

Have a look at the implementation of GridGenerator::half_hyper_ball [1]. 
You probably just want to use the first 4 vertices and then 
p,
p+Point<2>(0,-1) *(radius/std::sqrt(2.0)*a),
p+Point<2>(0,-1) *radius,
additionally.

Can you specify in formulas what symmetry boundary conditions mean for you?

Best,
Daniel


[1] 
https://www.dealii.org/developer/doxygen/deal.II/grid__generator_8cc_source.html

Am Donnerstag, 7. Juli 2016 00:08:29 UTC+2 schrieb benhour@gmail.com:
>
> Dear All,
> I want to define a new geometry (a quarter of a circle) and apply symmetry 
> boundary condition on two sides of the quarter. It would be very kind of 
> you if you help me for that.
>
> Best,
> Benhour
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: geometry

[Daniel Arndt]

Benhour,

Did you try what I suggested in 
https://groups.google.com/d/msg/dealii/QPLvZIQx7yw/LQ7t0-xWBgAJ ?

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: geometry

[Daniel Arndt]

Benhour,

It would be helpful if you tell us what the problem with the newly created 
problem is.
Looking at the implementation of GridGenerator::half_hyper_ball [1],
you can observe that there are cells which have vertices both with positive 
and negative
y-components. This will likely lead to undesired results just deleting 
cells.

It might still be easiest to modify [1] yourself for your purpose.

Best,
Daniel

[1] 
https://www.dealii.org/8.4.0/doxygen/deal.II/grid__generator_8cc_source.html



Am Samstag, 16. Juli 2016 04:10:33 UTC+2 schrieb benhour@gmail.com:
>
> Dear J-P,
> Thanks for your response. I used the code for creating a quarter of a 
> circle that comes as follow:
>
> Triangulation<2> triangulation;
> const Point<2> center;
> const double radius = 1.;
> GridGenerator::half_hyper_ball(triangulation, center, radius);
>
> Triangulation<2>::active_cell_iterator
> cell = triangulation.begin_active(),
> endc = triangulation.end();
> std::set::active_cell_iterator> remove;
> // Triangulation<2>::active_cell_iterator cell;
> for (cell=triangulation.begin_active(); cell!=triangulation.end(); ++cell)
> {
> Point<2> v=cell->center();
> // If it is in the bottom part of the geometry, remove it
> if (v(1) <= 0)
> remove.insert(cell);
> }
>  // remove the marked cells
>  GridGenerator::create_triangulation_with_removed_cells(triangulation, 
> removal, triangulation);
>
> Unfortunately, It does not create the right geometry. I do really 
> appreciate your kindness if you let me know which parts of this code go 
> wrong. Looking forward to hearing from you.
>
> Best,
> Benhour
>
> On Wednesday, July 13, 2016 at 6:15:04 AM UTC-5, Jean-Paul Pelteret wrote:
>>
>> Dear Benhour,
>>
>> You could use one of the GridGenerator options to create a half circle, 
>> and then use GridGenerator::create_triangulation_with_removed_cells 
>> 
>>  to 
>> remove the excess cells that yo udon't want.
>>
>> Regards,
>> J-P
>>
>> On Wednesday, July 13, 2016 at 3:40:02 AM UTC+2, benhour@gmail.com 
>> wrote:
>>>
>>> Dear All,
>>> I want to create a quarter of a circle in deal ii. In  this case, I 
>>> could create half of a circle with half_hyper_ball. I have tried 
>>> GridTools::delete_unused_vertices or GridTools::transformation to remove or 
>>> shift the same vertices to eliminate half of the circle. It would be very 
>>> kind of you if you let me know how to define this geometry in Dealii. It 
>>> should be noted that I also used quarter_hyper_shell and put inner radius 
>>> to zero but it did not work for this case.
>>>
>>> Best,
>>> Benhour
>>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Dirichlet constraint to a specific node

[Daniel Arndt]

Anup,

if you know which node you want to constrain, you just have to find the 
corresponding degree of freedom and add a constraint to the 
ConstraintMatrix.
For finding this DoF, DoFTools::map_dofs_to_support_points[1] might be 
useful. The alternative is to loop through all the cells and ask each DoF 
for its support_point [2] and component[3].

Best,
Daniel

[1] 
https://www.dealii.org/8.4.0/doxygen/deal.II/namespaceDoFTools.html#a10af6ba7adc43214e7b9f6815588f2aa
 

[2] 
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element
[3] 
https://dealii.org/8.4.1/doxygen/deal.II/classFiniteElement.html#a27220a135402b96c7e6eecbb04acda56

Am Samstag, 16. Juli 2016 03:55:29 UTC+2 schrieb Anup Basak:
>
> Hello,
>
> I have a vector valued (say, displacement) problem and I want to apply 
> Dirichlet boundary condition
> corresponding to one component of the displacement vector to a single node 
> (whose position is know)
> on a boundary.
>
> I shall be thankful if someone can  tell me how can I implement it in 
> dealii.
>
> Thank you,
>
> Anup.
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: geometry

[Daniel Arndt]

Benhour,

you can find a modification of step-3 in which the domain is a quarter 
circle attached.
The four cells from GridGenerator::half_hyper_ball have been replaced by 
just three. 
The vertex (0,-radius) has been removed and all the remaining vertices
with a negative y-component have been moved to y=0 while keeping the 
x-component.
Finally, the inner vertices have been moved a bit to increase the quality 
of the mesh.

You see that the modifications really are minor.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 1999 - 2015 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 at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Authors: Wolfgang Bangerth, 1999,
 *  Guido Kanschat, 2011
 */



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



class Step3
{
public:
  Step3 ();

  void run ();


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

  Triangulation<2> tria;
  FE_Q<2>  fe;
  DoFHandler<2>dof_handler;

  SparsityPattern  sparsity_pattern;
  SparseMatrix system_matrix;

  Vector   solution;
  Vector   system_rhs;
};


Step3::Step3 ()
  :
  fe (1),
  dof_handler (tria)
{}



void Step3::make_grid ()
{
  const unsigned int dim = 2;
  const double radius = 1.;
  const Point p;

  // equilibrate cell sizes at
  // transition from the inner part
  // to the radial cells
  const Point<2> vertices[7] = { p+Point<2>(0,0) *radius,
 p+Point<2>(+1,0) *radius,
 p+Point<2>(+1,0) *(radius/2),
 p+Point<2>(0,+1) *(radius/2),
 p+Point<2>(+1,+1) *(radius/(2*sqrt(2.0))),
 p+Point<2>(0,+1) *radius,
 p+Point<2>(+1,+1) *(radius/std::sqrt(2.0))
   };

  const int cell_vertices[3][4] = {{0, 2, 3, 4},
{1, 6, 2, 4},
{5, 3, 6, 4}
  };

  std::vector > cells (3, CellData<2>());

  for (unsigned int i=0; i<3; ++i)
{
  for (unsigned int j=0; j<4; ++j)
cells[i].vertices[j] = cell_vertices[i][j];
  cells[i].material_id = 0;
};

  tria.create_triangulation (
std::vector >(&vertices[0], &vertices[7]),
cells,
SubCellData());   // no boundary information

  Triangulation<2>::cell_iterator cell = tria.begin();
  Triangulation<2>::cell_iterator end = tria.end();

  while (cell != end)
{
  for (unsigned int i=0; i::faces_per_cell; ++i)
{
  if (cell->face(i)->boundary_id() == numbers::internal_face_boundary_id)
continue;

  // If x is zero, then this is part of the plane
  if (cell->face(i)->center()(0) < p(0)+1.e-5 * radius
  || cell->face(i)->center()(1) < p(1)+1.e-5 * radius)
cell->face(i)->set_boundary_id(1);
}
  ++cell;
}

  static const HyperBallBoundary boundary;
  tria.set_boundary (0, boundary);

  tria.refine_global (6);


  std::cout << "Number of active cells: "
<< tria.n_active_cells()
<< std::endl;
}




void Step3::setup_system ()
{
  dof_handler.distribute_dofs (fe);
  std::cout << "Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;

  DynamicSparsityPattern dsp(dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern (dof_handler, dsp);
  sparsity_pattern.copy_from(dsp);

  system_matrix.reinit (sparsity_pattern);

  solution.reinit (dof_handler.n_dofs());
  system_rhs.reinit (dof_handler.n_dofs());
}



void Step3::assemble_system ()
{
  QGauss<2>  quadrature_formula(2);
  FEValues<2> fe_values (fe, quadrature_formula,
 update_values | update_gradients | update_JxW_values);

  const unsigned int   dofs_per_cell

[deal.II] Re: Transferring solutions in distributed computing

[Daniel Arndt]

Junchao,

You want to use parallel::distributed::SolutionTransfer instead if you are 
on a parallel::distributed::Triangulation
Executing
$ grep -r "parallel::distributed::SolutionTransfer" .
in the examples folder, tells me that this object is used in step-32, 
step-42 and step-48.
Have for example a look at how this is done in step-42::refine_grid[1].

Best,
Daniel

[1] 
https://www.dealii.org/8.4.0/doxygen/deal.II/step_42.html#PlasticityContactProblemrefine_grid

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Transferring solutions in distributed computing

[Daniel Arndt]

Junchao,

It seems that the documentation is outdated for this piece of information.
In fact, neither PETScWrapper::MPI::Vector nor TrilinosWrappers::MPI::Vector
does have update_ghost_values.
What you should do is exactly what is done in the few lines of step-42 you 
referenced.
"solution = distributed_solution" imports ghost values while assigning the 
local values.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to apply boundary constraints in time-dependent and adaptive code?

[Daniel Arndt]

Junchao,

For what matters here, you are just interesting in how to deal with 
boundary constraints for a time-dependent problem.

I am not quite sure what your problem is. You have time-dependent boundary 
conditions so you want to use boundary
conditions corresponding to the correct time step. Why would you want to 
change the boundary conditions for an old solution?

The interpolated solution should already have the correct (old) boundary 
conditions.

Best,
Daniel


Am Freitag, 22. Juli 2016 00:34:10 UTC+2 schrieb Junchao Zhang:
>
> Hello,
>   I want to know how apply boundary constraints in a time-dependent, 
> adaptively refined and distributed memory code. In a time-step, I want 
> multiple refinements.  I could not find such an example in deal.II 
> tutorial. 
>   I have the following code, with questionable code in red.  Basically, 
> when transferring old solution, I need boundary constraints for one 
> time-step back. When assembling the system, I need boundary constraints for 
> the current time-step. But I don't know how to satisfy both.
>   
>   for (int time_step = 0; time_step < n_time_steps; time_step++) {
>time += dt;
>old_locally_relevant_solution = locally_relevant_solution;
>
>   for (int refine_step = 0; refine_step < n_refine_steps; 
> refine_step++) {
>  assemble_system(); // this function accesses 
> old_locally_relevant_solution and calls 
> constraints.distribute_local_to_global (...) to add local matrix to global 
> matrix.
>   // The boundary constraints 
> were computed for time - dt, i.e., one time-step back, however, here what 
> we need is boundary constraints for the current time step.
>  solve(); // update locally_relevant_solution
>
>  if (refine_step+1 < n_refine_steps) {
> // Estimate Kelly errors;
>...
> // Do triangulation, and transfer old solution from old mesh 
> to new mesh
>triangulation.prepare_coarsening_and_refinement();
>parallel::distributed::SolutionTransfer PETScWrappers::MPI::Vector> soltrans(dof_handler);
>   
>  
> soltrans.prepare_for_coarsening_and_refinement(old_locally_relevant_solution);
>triangulation.execute_coarsening_and_refinement ();
>
>setup_system(); // reinit vector, matrix and constraints for 
> the new mesh. The boundary function is computed for time - dt, since that 
> is the time when old_locally_relevant_solution was computed
>
>PETScWrappers::MPI::Vector 
> interpolated_solution(locally_owned_dofs, mpi_communicator);
>soltrans.interpolate(interpolated_solution);
>
>constraints.distribute(interpolated_solution); 
>old_locally_relevant_solution = interpolated_solution;
>  }
> }
>
>
>Thank you!
>
> --Junchao Zhang
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to apply boundary constraints in time-dependent and adaptive code?

[Daniel Arndt]

Junchao,


>> The interpolated solution should already have the correct (old) boundary 
>> conditions.
>>
> Is it true? As you see, I call 
> constraints.distribute(interpolated_solution) to impose boundary 
> constraints, which are computed in setup_system() as follows
>
What I meant to say is that the interpolated vector has the old boundary 
values if you don't apply additional constraints.
On cells that you refine, you can improve the approximation of the boundary 
values by applying a ConstraintMatrix (at time-dt) and hanging node 
constraints.

For the current time step, you then initialize a different ConstraintMatrix 
containing both Dirichlet boundary values (at time instead of time-dt) and 
hanging node constraints
and use this when assembling your system.

Best,
Daniel

>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to solve each block separately when working with TrilinosWrappers::BlockSparseMatrix and AMG precondition.

[Daniel Arndt]

Pasha, 

there are tutorial programs, such as step-22[1], that use BlockMatrices and 
solve for individual blocks.
What exactly is the problem you are facing? Does it have anything to do 
with Trilinos objects or is it a conceptual problem?

Can you tell us how exactly you want to split solving your system?

Best,
Daniel

[1] https://www.dealii.org/8.4.0/doxygen/deal.II/step_22.html

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to add values of solutions from different FE spaces (separate DoFHandlers)

[Daniel Arndt]

Jiaqi,

what you really want is that
(u_n_plus_1, v_h) = (u_stat-dt*dp)
holds for all v_h in your velocity ansatz space. This means inverting a 
mass matrix for the velocity
and may come down to just assigning function values depending on what your 
ansatz spaces are.

In that case, you can loop over all cells and get on each of them via 
fe_values.get_function_values
and fe_values.get_function_gradients all the information you need at all 
the velocity support points.
See [1] for how to setup an appropriate quadrature rule.

Otherwise you have to assemble the mass matrix and the right-hand side.

In any case, you want to use two FEValues objects and two cell iterators 
for velocity and pressure.
In particular, use fe_values_pressure.get_function gradients for the 
FEValues object initialized with the pressure cell_iterator
and fe_values_velocity.get_function_values for the FEValues obect 
initialized with the velocity cell_iterator.
The error (probably) tells you that you are trying to reinitialize an a 
FEValues object with a cell_iterator that does notr correspnd to the 
FiniteElement
it has originally been set up with.

Best,
Daniel

[1] 
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element

Am Montag, 25. Juli 2016 19:33:48 UTC+2 schrieb Jiaqi ZHANG:
>
> Dear all,
>
> I am trying to solve the NS equation using project methods, and follow the 
> idea of Step-35 to treat velocity and pressure separately,
>
> but when I am not sure how to implement the correction of the velocity 
> "u_n_plus_1 = u_star - dt * dp",
>
> where "u_star" and  "dp" are from the diffusion step and the projection 
> step respectively, 
>
> since velocity and pressure belong to different finite element spaces. 
>
> I tried "get_function_gradients", but end up with 
> The violated condition was: 
> static_cast&>(*this->fe) == 
> static_cast&>(cell->get_fe())
>
> Thanks in advance.
>
>
> Best,
> Jiaqi
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to solve each block separately when working with TrilinosWrappers::BlockSparseMatrix and AMG precondition.

[Daniel Arndt]

Pasha,


> *if* (direct_solver)
>
> {
>
> SolverControl cn;
>
> TrilinosWrappers::SolverDirect solver(cn);
>
> solver.*solve*(system_pde_matrix.block(0,0),
>
> newton_update.block(0), system_pde_residual.block(0));
>
>
> constraints_update.distribute(newton_update);
>
>
> *return* 1;
>
> }
>
>
>
> block(0,0) of system matrix and block(0) of residual and Newton update 
> vectors are selected but all block is solved simultaneously!! So I cant 
> solve specific block using direct solver and it is the first problem.
>
This is really strange. None of the objects you are using should see the 
other blocks. Can you confirm that the first block really does not contain 
the whole system? Can you just check the sizes of these objects?
 

>
> The second problem is when working with iterative solver. I cant set an 
> AMG Precondition for specific block of the system. For example I define a 
> class like this:
>
>
> *template* <*class* *PreconditionerA*>
>
> *class* BlockDiagonalPreconditioner_one
>
> {
>
> *public*:
>
> *BlockDiagonalPreconditioner_one*(*const* LA::MPI::BlockSparseMatrix &M,
>
> *const* *PreconditionerA* &pre_A)
>
> : matrix(M),
>
> prec_A (pre_A)
>
> {
>
> }
>
>
> *void* *vmult* (LA::MPI::BlockVector &dst,
>
> *const* LA::MPI::BlockVector &src) *const*
>
> {
>
> prec_A.vmult(dst.block(0), src.block(0));
>
> }
>
>
>
> *const* LA::MPI::BlockSparseMatrix &matrix;
>
> *const* *PreconditionerA* &prec_A;
>
> };
>
>
> for solving block(0,0) as,
>
>
> BlockDiagonalPreconditioner_one
>
> preconditioner(system_pde_matrix.block(0,0), preconditioner_solid);
>
>
> solver.*solve*(system_pde_matrix.block(0,0), newton_update.block(0),
>
> system_pde_residual.block(0), preconditioner);
>
>
> But there are two errors:
>
>
> 1) no matching function for call to 
> ‘BlockDiagonalPreconditioner_one::BlockDiagonalPreconditioner_one(dealii::BlockMatrixBase::BlockType&,
>  
> dealii::LinearAlgebraTrilinos: 
>
:MPI::PreconditionAMG&)’
>
The type of an individual block is Trilinos::Wrappers::SparseMatrix and not 
LA::MPI::BlockSparseMatrix. You need to change the definition in your 
Preconditioner class.
 

>
> 2) Invalid arguments '
>
> Candidates are:
>
> void solve(const #1 &, dealii::TrilinosWrappers::MPI::BlockVector &, 
> const dealii::TrilinosWrappers::MPI::BlockVector &, const #10001 &)
>
This hopefully goes away after you fixed the first error.
 

>
> the first error may be due to sending block(0,0) of TrilinosWrappers 
> BlockMatrix as first argument of *BlockDiagonalPreconditioner_one* 
> constructor. 
>
> Finally how I can set precondition for specific block and how I can solve 
> that block when working with TrilinosWrappers?
>
Doing a `$ grep -r "TrilinosWrappers::BlockSparseMatrix" .` tells me that 
step-31, step-32, step-43 and step-45 use 
TrilinosWrappers::BlockSparseMatrix
and all of these use preconditioners for specific blocks. That said, I 
would assume that what you are doing should work as well.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Step-42 -- Problem with Trilinos and P4EST

[Daniel Arndt]

Ehsan,

can you show us what your "detailed.log" in the build directory looks like?
This file stores the configuration that is used for building deal.II.
` -DCMAKE_INSTALL_PREFIX= /home/General_for_All_Users/deal_II_install_dir` 
should tell deal.II to copy the build library into that folder
when you are invoking `$make install`.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Using create_point_source_vector() function with hp::DoFHandler

[Daniel Arndt]

Deepak,

could you provide us with a minimal example that shows your problem?
Are you running in debug mode? What exactly does the program tell you when 
aborting?

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Using create_point_source_vector() function with hp::DoFHandler

[Daniel Arndt]

Deepak,

without a working example I can run, I don't see an obvious error.
Note that there are two tests, namely "numerics/create_point_source_hp" and 
"hp/vectors_point_source_hp_01"
that test `VectorTools::create_point_source` for hp. Maybe this helpful for 
debugging.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: A question about velocity correction in Step-35

[Daniel Arndt]

Jiaqi,

the algorithm in step-35 is the one that is proposed in 
"L. Timmermans, P. Minev, and F. Van De Vosse, “An approximate projec-
tion scheme for incompressible flow using spectral elements”, International
journal for numerical methods in fluids, vol. 22, no. 7, pp. 673–688, 1996."

The velocity that is solved for is the auxiliary one. You can recover the 
"true"
weakly divergence-free velocity by the correction you mentioned.
Notice that the divergence-free velocity is eliminated from the equation 
for the
auxiliary velocity.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Using the solution from one problem as a boundary condition in another problem with matching mesh on the boundary

[Daniel Arndt]

krei,

If you want to solve different PDEs on different domains that can be 
discretized by a common mesh, the preferred approach is to use a hp-vector 
finite element.
This means that on each of your subdomains all blocks of your finite 
element but one are of type FENothing.
You might want to have a look into step-46 for how to do this.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Floquet periodic conditions and complex valued algebra

[Daniel Arndt]

Daniel,

1.) You can use complex valued algebra if you split your problem into two 
equations. Apart from that most of the algebra objects can be used with 
std::complex.
At the moment, this is not true for constraints given by a ConstraintMatrix 
object, but this is WIP (https://github.com/dealii/dealii/issues/1760).
2.) Using SLEPc it should be able to compute eigenvalues for a complex 
valued system
3.) As already mentioned at the moment you can't use std::complex in a 
ConstraintMatrix. If however, you are able to formulate your problem as a 
system of real valued problems 
you can use whatever (linear) boundary condition you want. You can use 
make_periodicity_constraints to setup the DoFs you want to constraint and 
modify the inhomogeneity afterwards suitably.
4.) If you can formulate your problem as PDE (in not more than 3D) you can 
likely use deal.II for solving it numerically.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: A question about velocity correction in Step-35

[Daniel Arndt]

Wolfgang,

Interesting that this isn't even mentioned in the introduction to that 
> program, or anywhere else, for that matter. How did you find out? 
>
A major part of my PhD thesis is based on these projection methods.
Reading the literature this seems to be the first paper in which the 
rotational incremental
form has been considered.
 

> It's probably worthwhile adding this reference somewhere. I can do that if 
> you 
> want me to. 
>
I will create a PR next week including some more references to the analysis 
for that scheme.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: time-dependent boundary conditions for distributed computations

[Daniel Arndt]

Marek,

You would just set up your ConstraintMatrix anew. This is also the way it 
is done in step-44.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Step-42 -- Problem with Trilinos and P4EST

[Daniel Arndt]

Ehsan,

Good that you found the error. Another possibility is to use ccmake.
This works for me quite well.

Best,
Daniel


Am Mittwoch, 3. August 2016 11:37:14 UTC+2 schrieb Ehsan:
>
> Finally I found the problem after installing and uninstalling for more 
> than ten times.
> The problem was the space after "=" !!
> There should be NO space before and after "=" in setting cmake variable.
>
>
> On Monday, August 1, 2016 at 3:49:56 PM UTC+2, Jean-Paul Pelteret wrote:
>>
>> Dear Eshan,
>>
>> Sometimes it is necessary to manually remove CMakeCache.txt before 
>> reconfiguring a project. Can you try to do this and see if you have any 
>> further success?
>>
>> Regards,
>> J-P
>>
>> On Monday, August 1, 2016 at 2:42:25 PM UTC+2, Ehsan wrote:
>>>
>>> Dear Daniel,
>>> I checked the "detailed.log" file and noticed that DEAL_II_WITH_MPI, 
>>> DEAL_II_WITH_P4EST and DEAL_II_WITH_TRILINOS all three are off.
>>> my cmake is:
>>> cmake 
>>> -DCMAKE_INSTALL_PREFIX=/home/General_for_All_Users/deal_II_install_dir 
>>> -DDEAL_II_WITH_64BIT_INDICES= ON 
>>> -DDEAL_II_WITH_MPI= ON 
>>> -DDEAL_II_WITH_P4EST= ON 
>>> -DDEAL_II_WITH_TRILINOS= ON 
>>> -DP4EST_DIR=/home/General_for_All_Users/P4EST_install_dir/ 
>>> -DTRILINOS_DIR=/home/General_for_All_Users/Trilinos_install_dir/trilinos_path/
>>>  
>>>
>>> ..
>>>
>>> I also attached the detailed.log file.
>>> When I put -DCMAKE_INSTALL_PREFIX before options related to MPI, P4EST 
>>> and TRILINOS, it is implemented.
>>> So something is wrong with setting MPI, P4EST or TRILINOS but I do not 
>>> know what !!!
>>>
>>> best regards.
>>> Ehsan
>>>
>>>
>>> On Monday, August 1, 2016 at 2:32:22 PM UTC+2, Ehsan wrote:
>>>>
>>>> Dear Daniel,
>>>>
>>>> When I run cmake in terminal and the end I receive below warning:
>>>>
>>>> CMake Warning:
>>>>   Manually-specified variables were not used by the project:
>>>>
>>>> P4EST_DIR
>>>> TRILINOS_DIR
>>>>
>>>> How should I solve this problem.
>>>> Best regards.
>>>> Ehsan
>>>>
>>>> On Monday, August 1, 2016 at 2:21:31 PM UTC+2, Ehsan wrote:
>>>>>
>>>>> Dear Daniel,
>>>>>
>>>>> In the "detailed.log" the "CMAKE_INSTALL_PREFIX:" is empty but I am 
>>>>> sure that I defined it in cmake:
>>>>> I have a generat question.
>>>>> is it possible to compile deal ii with enabled Trilinos, P4EST, and 
>>>>> MPI options?
>>>>>
>>>>> Best regards.
>>>>> Ehsan
>>>>>
>>>>>
>>>>>
>>>>> On Wednesday, July 27, 2016 at 8:52:36 PM UTC+2, Daniel Arndt wrote:
>>>>>>
>>>>>> Ehsan,
>>>>>>
>>>>>> can you show us what your "detailed.log" in the build directory looks 
>>>>>> like?
>>>>>> This file stores the configuration that is used for building deal.II.
>>>>>> ` -DCMAKE_INSTALL_PREFIX= /home/General_for_All_Users/
>>>>>> deal_II_install_dir` should tell deal.II to copy the build library 
>>>>>> into that folder
>>>>>> when you are invoking `$make install`.
>>>>>>
>>>>>> Best,
>>>>>> Daniel
>>>>>>
>>>>>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

dealii@googlegroups.com

[Daniel Arndt]

It seems that there is something wrong with your TBB installation. Which 
version of TBB are you using?
Could you send your "detailed.log" file located in the deal.II build 
directory?

TBB is also bundled with deal.II and you can force to use that via
$ cmake -DDEAL_II_FORCE_BUNDLED_THREADS=ON .
in your build directory. Does this work for you?

Best,
Daniel

Am Samstag, 6. August 2016 06:17:11 UTC+2 schrieb deal.II newbie:
>
> Hi!
>
> After installing Deal.II into my machine, it has not passed the quick 
> tests. I have checked the quicktests.log, the error message for all four 
> tests are as below,
>
> Start 1: step.debug
> 1/4 Test #1: step.debug ...***Failed   12.15 sec
> Test step.debug: BUILD
> ===   OUTPUT BEGIN 
>  ===
> step.debug: BUILD failed. Output:
> [  0%] Built target kill-step.debug-OK
> [  0%] Built target expand_instantiations_exe
> [  0%] Built target obj_opencascade.inst
> [  0%] Built target obj_opencascade.debug
> [  5%] Built target obj_boost_serialization.debug
> [  5%] Built target obj_boost_system.debug
> [  7%] Built target obj_boost_iostreams.debug
> [  9%] Built target obj_amd_int.debug
> [ 18%] Built target obj_umfpack_I_UMF.debug
> [ 25%] Built target obj_umfpack_L_UMF.debug
> [ 31%] Built target obj_umfpack_I_UMFPACK.debug
> [ 35%] Built target obj_umfpack_L_UMFPACK.debug
> [ 35%] Built target obj_umfpack_DI_TSOLVE.debug
> [ 35%] Built target obj_umfpack_DI_TRIPLET_MAP_NOX.debug
> [ 35%] Built target obj_umfpack_DI_TRIPLET_MAP_X.debug
> [ 37%] Built target obj_umfpack_DI_TRIPLET_NOMAP_X.debug
> [ 37%] Built target obj_umfpack_DI_STORE.debug
> [ 37%] Built target obj_umfpack_DI_ASSEMBLE.debug
> [ 37%] Built target obj_umfpack_DI_SOLVE.debug
> [ 37%] Built target obj_umfpack_DL_TSOLVE.debug
> [ 37%] Built target obj_umfpack_DL_TRIPLET_MAP_NOX.debug
> [ 37%] Built target obj_umfpack_DL_TRIPLET_MAP_X.debug
> [ 37%] Built target obj_umfpack_DL_TRIPLET_NOMAP_X.debug
> [ 38%] Built target obj_umfpack_DL_STORE.debug
> [ 38%] Built target obj_umfpack_DL_ASSEMBLE.debug
> [ 38%] Built target obj_umfpack_DL_SOLVE.debug
> [ 38%] Built target obj_umfpack_GENERIC.debug
> [ 40%] Built target obj_amd_long.debug
> [ 40%] Built target obj_amd_global.debug
> [ 42%] Built target obj_muparser.debug
> [ 48%] Built target obj_fe.inst
> [ 55%] Built target obj_fe.debug
> [ 57%] Built target obj_dofs.inst
> [ 59%] Built target obj_dofs.debug
> [ 62%] Built target obj_numerics.inst
> [ 70%] Built target obj_numerics.debug
> [ 72%] Built target obj_lac.inst
> [ 77%] Built target obj_lac.debug
> [ 79%] Built target obj_base.inst
> [ 88%] Built target obj_base.debug
> [ 90%] Built target obj_grid.inst
> [ 94%] Built target obj_grid.debug
> [ 94%] Built target obj_hp.inst
> [ 96%] Built target obj_hp.debug
> [ 98%] Built target obj_multigrid.inst
> [ 98%] Built target obj_multigrid.debug
> [100%] Built target obj_distributed.inst
> [100%] Built target obj_distributed.debug
> [100%] Built target obj_algorithms.inst
> [100%] Built target obj_algorithms.debug
> [100%] Built target obj_integrators.debug
> [100%] Built target obj_matrix_free.inst
> [100%] Built target obj_matrix_free.debug
> [100%] Built target obj_meshworker.inst
> [100%] Built target obj_meshworker.debug
> [100%] Built target deal_II.g
> Linking CXX executable step.debug
> ../../lib/libdeal_II.g.so.8.4.1: undefined reference to 
> `tbb::interface5::internal::task_base::destroy(tbb::task&)'
> collect2: error: ld returned 1 exit status
> gmake[7]: *** [tests/quick_tests/step.debug] Error 1
> gmake[6]: *** [tests/quick_tests/CMakeFiles/step.debug.dir/all] Error 2
> gmake[5]: *** [tests/quick_tests/CMakeFiles/step.debug.run.dir/rule] Error 
> 2
> gmake[4]: *** [step.debug.run] Error 2
>
>
> step.debug: **BUILD failed***
>
> ===OUTPUT END   
> ===
> Expected stage PASSED - aborting
> CMake Error at /root/dealii-8.4.1/cmake/scripts/run_test.cmake:140 
> (MESSAGE):
>   *** abort
>
> As indicated in the log file, I guess there may be some problem with the 
> TBB. I would like to ask if anyone could please teach me how to fix this 
> error? I very much appreciate any help and hints, and thank you very much 
> in advance!
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Post-processing with different cell material properties

[Daniel Arndt]

Pete,

I remember someone having something like this in mind. This is the link to 
that discussion:
https://groups.google.com/forum/#!topic/dealii/QKuDndKojug

A different way would be to derive from DataOut and overwrite 
DataOut::first_cell and DataOut::next_cell by skipping all the cells that 
don't have a specific material_id.
If you do this for each material_id, you could create an output file for 
each material_id and join them afterwards. This would be similar to how you 
combine output from different MPI processes using a pvtu file (step-40).

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Embedding a Python script

[Daniel Arndt]

Oded,

If you include 

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

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

into your CMakeLists.txt, you can switch between debug and release mode 
using
$make debug
and
$make release

Best,
Daniel

Am Dienstag, 9. August 2016 10:36:39 UTC+2 schrieb Oded Yaakobi:
>
>
> Hi Denis,
>
> I tried your suggestion to drop DEAL_II_INVOKE_AUTOPILOT(), and it worked.
> However, if I want to keep the functionality of switching the build types 
> between "Debug" and "Release" versions, how can I do it now?
>
> Thanks,
> Oded
>
>
> On Monday, August 8, 2016 at 11:08:14 PM UTC+3, Denis Davydov wrote:
>>
>> Hi Ode,
>>
>> First, Spack is a tool to compile and install packages, which I don't 
>> think you are using as you said that
>> you play around with Luca's build. So the link to this wiki should not be 
>> relevant here.
>>
>> Second, you CMake input is wrong. you need to drop
>>
>> DEAL_II_INVOKE_AUTOPILOT()
>>
>> see the very first example here 
>> https://www.dealii.org/developer/users/cmakelists.html#cmakesimple
>>
>> My understanding is that DEAL_II_INVOKE_AUTOPILOT is for very simple 
>> cases like examples/ which don't
>> need to link against other than deal.II libraries and where one needs a 
>> single executable.
>> If you are working on something like you PhD code, then it's most likely 
>> better to keps things are
>> you library + unit tests + executable which calls you main class with 
>> parameter file.
>>
>> Regards,
>> Denis.
>>
>> On Monday, August 8, 2016 at 12:54:02 PM UTC+2, Oded Yaakobi wrote:
>>>
>>> Dear group,
>>>
>>>  
>>>
>>> I apologize if my question below has already been addressed in the one 
>>> of deal.ii’s documentation resources, but I didn't understand how to answer 
>>> it using the pieces of documentation that I have read such as
>>>
>>> “How to use CMake to configure your projects with deal.II” 
>>> https://www.dealii.org/developer/users/cmakelists.html#cmakeauto
>>>
>>>  
>>>
>>> and 
>>>
>>> “deal.II in Spack” 
>>> https://github.com/dealii/dealii/wiki/deal.II-in-Spack
>>>
>>>  
>>>
>>> I am trying to embed a short Python script in my code, but get an error 
>>> message after running “make release”.  Attached are an example code with a 
>>> minor addition to the “main” of Step-1 and the “CMakeLists.txt” file that I 
>>> tried to adapt to support Python.
>>>
>>>  
>>>
>>> I am running my code on a MAC OS Yosemite 10.10.5 with the special 
>>> version of deal.ii that Luca prepared for me 
>>> https://groups.google.com/forum/#!searchin/dealii/%22Make$20test%22$20failure$20following$20Trilinos$20installation$20with$20Mesquite%7Csort:relevance/dealii/Y7MH_BE0yko/1zvnHGR7AQAJ
>>>  
>>> 
>>>
>>>  
>>>
>>>
>>> Here are the error messages that I get:
>>>
>>>  
>>>
>>> [ yaakobioy L02029080  
>>> ~/dealii/dealii-8-5-0pre-v2/dealii-8-5-0pre-v2-Install/examples/step-1-mod 
>>> ]$ make release
>>>
>>> CMake Error at 
>>> /Users/yaakobioy/dealii/dealii-8-5-0pre-v2/dealii-8-5-0pre-v2-Install/share/deal.II/macros/macro_deal_ii_invoke_autopilot.cmake:57
>>>  
>>> (ADD_EXECUTABLE):
>>>
>>>   add_executable cannot create target "step-1" because another target 
>>> with
>>>
>>>   the same name already exists.  The existing target is an executable 
>>> created
>>>
>>>   in source directory
>>>
>>>   
>>> "/Users/yaakobioy/dealii/dealii-8-5-0pre-v2/dealii-8-5-0pre-v2-Install/examples/step-1-mod".
>>>
>>>   See documentation for policy CMP0002 for more details.
>>>
>>> Call Stack (most recent call first):
>>>
>>>   CMakeLists.txt:49 (DEAL_II_INVOKE_AUTOPILOT)
>>>
>>>  
>>>
>>>  
>>>
>>> -- Configuring incomplete, errors occurred!
>>>
>>> See also 
>>> "/Users/yaakobioy/dealii/dealii-8-5-0pre-v2/dealii-8-5-0pre-v2-Install/examples/step-1-mod/CMakeFiles/CMakeOutput.log".
>>>
>>> make: *** [cmake_check_build_system] Error 1
>>>
>>>  
>>>
>>>  
>>>
>>>  
>>>
>>> I would appreciate if someone could explain what should I do to overcome 
>>> this problem.
>>>
>>>  
>>>
>>> Thanks in advance,
>>> Oded
>>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Embedding a Python script

[Daniel Arndt]

Oded,

Your CMakeLists.txt file works for me. Maybe you just didn't run cmake 
again?
Can you try on a clean directory just containing CMakeLists.txt and 
step-1.cc?

Best,
Daniel

Am Dienstag, 9. August 2016 11:38:54 UTC+2 schrieb Oded Yaakobi:
>
> HI Daniel,
>
> I tried your suggestion but it didn't work for me. See the attached 
> "CMakeLists.txt" file. I get the following error messages:
>
> [ yaakobioy L02029080  
> ~/dealii/dealii-8-5-0pre-v2/dealii-8-5-0pre-v2-Install/examples/step-1-mod 
> ]$ make debug
>
> make: *** No rule to make target `debug'.  Stop.
>
> [ yaakobioy L02029080  
> ~/dealii/dealii-8-5-0pre-v2/dealii-8-5-0pre-v2-Install/examples/step-1-mod 
> ]$ make release
>
> make: *** No rule to make target `release'.  Stop.
>
> Thanks,
> Oded
>
>
> On Tuesday, August 9, 2016 at 12:00:22 PM UTC+3, Daniel Arndt wrote:
>>
>> Oded,
>>
>> If you include 
>>
>> ADD_CUSTOM_TARGET(debug
>>   COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR}
>>   COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
>>   COMMENT "Switch CMAKE_BUILD_TYPE to Debug"
>>   )
>>
>> ADD_CUSTOM_TARGET(release
>>   COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR}
>>   COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all
>>   COMMENT "Switch CMAKE_BUILD_TYPE to Release"
>>   )
>>
>> into your CMakeLists.txt, you can switch between debug and release mode 
>> using
>> $make debug
>> and
>> $make release
>>
>> Best,
>> Daniel
>>
>> Am Dienstag, 9. August 2016 10:36:39 UTC+2 schrieb Oded Yaakobi:
>>>
>>>
>>> Hi Denis,
>>>
>>> I tried your suggestion to drop DEAL_II_INVOKE_AUTOPILOT(), and it 
>>> worked.
>>> However, if I want to keep the functionality of switching the build 
>>> types between "Debug" and "Release" versions, how can I do it now?
>>>
>>> Thanks,
>>> Oded
>>>
>>>
>>> On Monday, August 8, 2016 at 11:08:14 PM UTC+3, Denis Davydov wrote:
>>>>
>>>> Hi Ode,
>>>>
>>>> First, Spack is a tool to compile and install packages, which I don't 
>>>> think you are using as you said that
>>>> you play around with Luca's build. So the link to this wiki should not 
>>>> be relevant here.
>>>>
>>>> Second, you CMake input is wrong. you need to drop
>>>>
>>>> DEAL_II_INVOKE_AUTOPILOT()
>>>>
>>>> see the very first example here 
>>>> https://www.dealii.org/developer/users/cmakelists.html#cmakesimple
>>>>
>>>> My understanding is that DEAL_II_INVOKE_AUTOPILOT is for very simple 
>>>> cases like examples/ which don't
>>>> need to link against other than deal.II libraries and where one needs a 
>>>> single executable.
>>>> If you are working on something like you PhD code, then it's most 
>>>> likely better to keps things are
>>>> you library + unit tests + executable which calls you main class with 
>>>> parameter file.
>>>>
>>>> Regards,
>>>> Denis.
>>>>
>>>> On Monday, August 8, 2016 at 12:54:02 PM UTC+2, Oded Yaakobi wrote:
>>>>>
>>>>> Dear group,
>>>>>
>>>>>  
>>>>>
>>>>> I apologize if my question below has already been addressed in the one 
>>>>> of deal.ii’s documentation resources, but I didn't understand how to 
>>>>> answer 
>>>>> it using the pieces of documentation that I have read such as
>>>>>
>>>>> “How to use CMake to configure your projects with deal.II” 
>>>>> https://www.dealii.org/developer/users/cmakelists.html#cmakeauto
>>>>>
>>>>>  
>>>>>
>>>>> and 
>>>>>
>>>>> “deal.II in Spack” 
>>>>> https://github.com/dealii/dealii/wiki/deal.II-in-Spack
>>>>>
>>>>>  
>>>>>
>>>>> I am trying to embed a short Python script in my code, but get an 
>>>>> error message after running “make release”.  Attached are an example code 
>>>>> with a minor addition to the “main” of Step-1 and the “CMakeLists.txt” 
>>>>> file 
>>>>> that I tried to adapt to support Python.
>>>>>
>&

[deal.II] Re: Using the solution from one problem as a boundary condition in another problem with matching mesh on the boundary

[Daniel Arndt]

krei,

If your emission current boundary conditions do not depend linearly on the 
electric field, the whole problem becomes non-linear and hence you can't 
solve the whole problem directly.
What you can do is to first solve for the electric field and afterwards for 
the metal part. In particular, this means to neglect in the first assembly 
all the cells that a related to the metal part.
After solving you would then loop over the interface to find out which dofs 
on the metal part you want to constrain. For this you would initialize a 
Quadrature object with the unit support_points [1].
At the same time you can use FEValues::get_function_* to evaluate the 
electric field on the locations of these dofs. This should be much faster 
than VectorTools::PointGradient.
Finally, assemble and solve for the metal part taking the computed 
constraints into account.

Does this make sense to you?

Best,
Daniel

[1] 
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element

Am Dienstag, 9. August 2016 16:06:04 UTC+2 schrieb krei:
>
> Hello,
>
> I mostly implemented the hp-vector finite element approach (according to 
> step-46), but alas, I think it might not be applicable. (I simplified the 
> boundary condition in original post a bit.) In my case I need to apply an 
> emission current boundary condition to the electric currents in copper 
> based on the electric field in vacuum. The emission currents are not 
> expressed through the electric field by a simple manner, more specifically, 
> J ~ E^2, J ~ 1/E and J ~ exp(a*E)  and I might want to use interpolation 
> from a grid to evaluate J(E). (J - emission current density, E - electric 
> field at the surface).
>
> If I want to solve everything in one big system (as is done in step-46), 
> then I don't have access to the electric field values during the system 
> assembly, I can access shape functions but I can't evaluate the emission 
> current through them.
>
> Perhaps I have misunderstood and I can somehow evaluate the boundary 
> conditions? Or what would a better approach be? 
>
> On Saturday, July 30, 2016 at 1:38:06 PM UTC+3, Daniel Arndt wrote:
>>
>> krei,
>>
>> If you want to solve different PDEs on different domains that can be 
>> discretized by a common mesh, the preferred approach is to use a hp-vector 
>> finite element.
>> This means that on each of your subdomains all blocks of your finite 
>> element but one are of type FENothing.
>> You might want to have a look into step-46 for how to do this.
>>
>> Best,
>> Daniel
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: tests fail

[Daniel Arndt]

Giorgis,

This looks like a problem with your PETSc installation. Can you confirm 
that it is installed correctly?
Do any problems appear if you compile deal.II without PETSc?

Best,
Daniel

Am Donnerstag, 11. August 2016 22:16:49 UTC+2 schrieb Giorgos Kourakos:
>
> Hi,
>
> I'm trying to compile a recent version of deal and I'm getting the 
> following error during the tests as well as when I compile the step-1 
>
> ../../lib/libdeal_II.g.so.8.4.1: error: undefined reference to 
> 'PCHYPRESetType'
>
>
> Any ideas what is wrong?
>
> thank you
>
> Giorgos
>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

Does the equation you want to solve for have multiple components? Otherwise 
it would not make sense to multiply in the right hand side something with a 
gradient.
If you want to project the gradient into a vector valued finite element 
ansatz space, then the matrix you are assembling looks weird.
In what you are doing you are also considering the coupling between all 
components.

You are calculating the gradient of your scalar finite element function 
correctly, but you need to find the correct component to use.
For doing this you can use something like:

const unsigned int component = fe.system_to_component_index(i).first;
cell_rhs(i) += ((fe_values.shape_value (i, q_index) *
fe_function_grad[q_index][component])*
fe_values.JxW (q_index));

If you want to assemble a mass matrix for a vector-valued finite element, 
you have also to restrict assembling for the matrix to the case
where the component for the ith and jth ansatz function are the same.

You might want to have a look at step-8 [1] and the module "Handling vector 
valued problems"[2].

Best,
Daniel

[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/step_8.html#ElasticProblemassemble_system
[2] https://www.dealii.org/8.4.1/doxygen/deal.II/group__vector__valued.html

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

Yes the matrix you are assembling is a vector-valued mass matrix now.
For me your code is failing with

void dealii::FEValuesBase::get_function_gradients(const 
InputVector&, std::vector >&) const [with InputVector = 
dealii::Vector; int dim = 3; int spacedim = 3; typename 
InputVector::value_type = double]
The violated condition was: 
(fe->n_components()) == (1)
The name and call sequence of the exception was:
dealii::ExcDimensionMismatch((fe->n_components()),(1))
Additional Information: 
Dimension 3 not equal to 1

and this is not surprising. What you are doing is to extract the gradients 
of a vector-valued finite element solution. The object that should store 
this should therefore be a std::vector> as you 
want to store for each quadrature point and each component a Tensor<1,dim>.

What you want to do is really that in "Do not work". As you have two 
DoFHandlers, you should also have two FEValues objects and two 
cell_iterator corresponding to the correct DoFHandler. Then you can extract 
the gradient of your scalar field and project it onto the ansatz space for 
the vector-valued field. This should look like:

typename DoFHandler::active_cell_iterator cell_scalar = 
dof_handler_scalar.begin_active();
typename DoFHandler::active_cell_iterator cell = 
dof_handler.begin_active();
...

 for (; cell!=endc; ++cell, ++cell_vector)
{
  fe_values.reinit (cell);
  fe_values_scalar.reinit (cell_scalar);
  fe_values_scalar.get_function_gradients(test,fe_function_grad);
  ...
}

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Data exchanges in IO

[Daniel Arndt]

Bruce,

The p_local_solution vector is filled with values as the result of solving 
> a linear equation, p_local_rho is filled with values by looping over 
> individual elements and assigning values calculated from other quantities. 
> Part of the reason the p_local_rho vector doesn't have ghost nodes is that 
> you can't access individual elements for a ghosted vector.
>
You can't write to a ghosted vector, but you can read from it. Is this what 
you mean?

The DataOut object appears to only handle ghosted vectors. It also appears 
> that you need to call compress on non-ghosted vectors (this was determined 
> empirically based on error messages). 
>
Yes, for proper output you should use ghosted vectors that store the 
solution on all the dofs of locally owned cells. If you modify a parallel 
vector you have to call compress as mentioned in step-40 [1] for example. 
Would you want to have this information at some different place as well. If 
yes, where?
 

> The vecScale function is just designed to multiply all values of an 
> un-ghosted vector by a constant. 
>
You could also use operator*= [2] for this.

If you don't change the values of the vectors using vecScale then the data 
> for "Phi" and "Rho" appears smooth (but the scale is unwieldy). If you 
> change the values of the vectors using vecScale, then "Phi" looks okay but 
> "Rho" appears to have unscaled values at locations that look like they are 
> at processor boundaries. I'm puzzled about how this is happening, since my 
> understanding is that a data exchange happens when you do the assignment 
> rho = p_local_rho. The vecScale function only operates on un-ghosted 
> vectors, but is it possible that it is missing some values somewhere?
>
This sounds weird as you do the exact same thing to (rho, p_local_rho) and 
(p_local_solution, phi). So you observe that p_local_solution looks wrong 
as well? Do you have a minimal working example showing this problem?

Best,
Daniel

[1] 
https://dealii.org/8.4.0/doxygen/deal.II/step_40.html#LaplaceProblemassemble_system
 

[2] 
https://www.dealii.org/8.4.0/doxygen/deal.II/classPETScWrappers_1_1VectorBase.html#a37900779c6049418c39bacc1d44f4260

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Derivative of a scalar function in the code

[Daniel Arndt]

Benhour,

What exactly is your problem? Lots of the examples programs, i.e. step-14, 
step-15, step-18, step-31, step-32, step-35, step-43 and step-44, use 
FEValues::get_function_gradients.

Best,
Daniel

Am Montag, 15. August 2016 03:13:14 UTC+2 schrieb benhour@gmail.com:
>
> Dear All,
> I want to extract the first derivative and the gradient of a scalar 
> function in dealii. Maybe I should use the void FEValuesViews::Scalar 
> <
>  
> dim, spacedim >::get_function_gradients. However, I really do not kn ow how 
> I should use it in a real code. I do really appreciate it if you help me to 
> write it.
>
> Best,
> Benhour
>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Data exchanges in IO

[Daniel Arndt]

Bruce,

[...]  3. Like in 2., but for adding values to individual elements. Use 
> VectorOperation::add.
>
> I guess 3) is obscure. Does this mean that you are adding a new element 
> where none existed before or that you are incrementing an existing element? 
> If the former, aren't the total number of elements in a vector specified 
> when you declare the vector size so how do you add more elements?
>
This does not mean that you want to modify the size of a vector or the 
underlying distribution pattern. It just means that after you modify some 
element of a vector using its original value, e.g. vec(5)+=2, you need to 
call compress(VectorOperation::add).
 

>
> The results coming from p_local_solution look good (smooth) and appear to 
> have been scaled, only the results for p_local_rho have a problem. I've 
> tried to create a simple reproducer but so far haven't had any success. 
> I'll keep trying.
>
As you already figured, DataOut needs a ghosted vector containing all the 
elements on locally owned cells for proper output. If I read your code 
correctly, only p_local_solution and rho have ghost elements while 
p_local_rho and phi are completely distributed, i.e. the don't have ghost 
elements. Hence, I am not surprised that the visualization of p_local_rho 
looks incorrect.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Data exchanges in IO

[Daniel Arndt]

Bruce,

You were just missing
p_rho.compress(VectorOperation::insert);
after assemble_system().

Best,
Daniel

Am Mittwoch, 17. August 2016 21:01:22 UTC+2 schrieb bjpal...@gmail.com:
>
>
>
> On Tuesday, August 16, 2016 at 2:47:18 AM UTC-7, Daniel Arndt wrote:
>>
>> Bruce,
>>
>> [...]  3. Like in 2., but for adding values to individual elements. Use 
>>> VectorOperation::add.
>>>
>>> I guess 3) is obscure. Does this mean that you are adding a new element 
>>> where none existed before or that you are incrementing an existing element? 
>>> If the former, aren't the total number of elements in a vector specified 
>>> when you declare the vector size so how do you add more elements?
>>>
>> This does not mean that you want to modify the size of a vector or the 
>> underlying distribution pattern. It just means that after you modify some 
>> element of a vector using its original value, e.g. vec(5)+=2, you need to 
>> call compress(VectorOperation::add).
>>  
>>
>>>
>>> The results coming from p_local_solution look good (smooth) and appear 
>>> to have been scaled, only the results for p_local_rho have a problem. I've 
>>> tried to create a simple reproducer but so far haven't had any success. 
>>> I'll keep trying.
>>>
>> As you already figured, DataOut needs a ghosted vector containing all the 
>> elements on locally owned cells for proper output. If I read your code 
>> correctly, only p_local_solution and rho have ghost elements while 
>> p_local_rho and phi are completely distributed, i.e. the don't have ghost 
>> elements. Hence, I am not surprised that the visualization of p_local_rho 
>> looks incorrect.
>>
> Actually, I'm adding rho to the DataOut object, not p_local_rho, so all 
> vectors added to the DataOut object are ghosted. The problem seems to be 
> when I scale values. If I use *= to modify the values, then it works okay, 
> but if I use my own hand-rolled code to do it, it fails to update the 
> processor boundaries. I've attached a code that shows what is happening. It 
> is a little on the long side but the important stuff is in the vecScale and 
> write_grid routines. At this point, I have a working solution, but I'd be 
> interested in understanding why my approach fails.
>
> Bruce
>
>
>
>> Best,
>> Daniel
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: [dealii-developers] dealii.org website down?

[Daniel Arndt]

dealii.org is online again.

Am Mittwoch, 17. August 2016 19:33:48 UTC+2 schrieb Wolfgang Bangerth:
>
> On 08/17/2016 11:30 AM, thomas stephens wrote: 
> > Currently unable to load https://www.dealii.org/ along with tutorials 
> > (problem began at approximately 12:30pm ETD) 
>
> Yes, there is a power outage on campus in Heidelberg, Germany, where the 
> website is hosted. Guido says that people are working on it -- expect 
> things to come back up online in the (hopefully) not too far future. 
>
> 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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Getting number of hanging support points

[Daniel Arndt]

Deepak,

I am aware of the fact that ConstraintMatrix.n_constraints() gives the 
> number of hanging nodes for h-refinement, but if I only use p-refinement, 
> can it give the number of hanging support points occurring due to the 
> difference in the order of bases between two adjacent elements? I am 
> expecting this based on the description of the hp-refinement document.
>
> I have been trying this but always get 0 hanging support points, although 
> I know that the some elements of the mesh use Q1, others Q2.
>
That sounds definitely strange. step-27 uses a hp::DoFHandler and 
make_hanging_node_constraints. Do you see that constraints are created 
there?
If you can't find a solution looking at that example program, can you 
reduce your code to a minimal example showing that no constraints are 
created for p-refinement?

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

You should not rely on any particular sorting of the dofs in the solution 
vector. Instead you can ask the FiniteElement on each cell for the 
component of a local dof by
const unsigned int component = fe.system_to_component_index(i).first; 

Apart from that the DataOut object knows about components and you can do 
some postprocessing on the solution vector using DataPostprocessor [1]

Best,
Daniel

 [1] https://dealii.org/8.4.1/doxygen/deal.II/classDataPostprocessor.html

Am Montag, 22. August 2016 17:05:08 UTC+2 schrieb Joel Davidsson:
>
> Dear Daniel,
>
> Thank you for a very good answer, adding the fe_values_scalar and 
> cell_scalar fixed the problem.
>
> I have a follow-up question about the solution I get out, how is the data 
> organized in the solution vector? Say for example I want to loop over all 
> the x components, how would I do that?
>
> Best,
> Joel
>
> On Friday, August 12, 2016 at 5:35:07 PM UTC+2, Daniel Arndt wrote:
>>
>> Joel,
>>
>> Yes the matrix you are assembling is a vector-valued mass matrix now.
>> For me your code is failing with
>>
>> void dealii::FEValuesBase::get_function_gradients(const 
>> InputVector&, std::vector> InputVector::value_type> >&) const [with InputVector = 
>> dealii::Vector; int dim = 3; int spacedim = 3; typename 
>> InputVector::value_type = double]
>> The violated condition was: 
>> (fe->n_components()) == (1)
>> The name and call sequence of the exception was:
>> dealii::ExcDimensionMismatch((fe->n_components()),(1))
>> Additional Information: 
>> Dimension 3 not equal to 1
>>
>> and this is not surprising. What you are doing is to extract the 
>> gradients of a vector-valued finite element solution. The object that 
>> should store this should therefore be a 
>> std::vector> as you want to store for each 
>> quadrature point and each component a Tensor<1,dim>.
>>
>> What you want to do is really that in "Do not work". As you have two 
>> DoFHandlers, you should also have two FEValues objects and two 
>> cell_iterator corresponding to the correct DoFHandler. Then you can extract 
>> the gradient of your scalar field and project it onto the ansatz space for 
>> the vector-valued field. This should look like:
>>
>> typename DoFHandler::active_cell_iterator cell_scalar = 
>> dof_handler_scalar.begin_active();
>> typename DoFHandler::active_cell_iterator cell = 
>> dof_handler.begin_active();
>> ...
>>
>>  for (; cell!=endc; ++cell, ++cell_vector)
>> {
>>   fe_values.reinit (cell);
>>   fe_values_scalar.reinit (cell_scalar);
>>   fe_values_scalar.get_function_gradients(test,fe_function_grad);
>>   ...
>> }
>>
>> Best,
>> Daniel
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: where and how should I use the MeshWorker class? how should I do to implement DG without MeshWorker?

[Daniel Arndt]

An old version of step-12[1] has an implementation that doesn't use 
MeshWorker.
In my opinion, it is covenient to use MeshWorker if you can when you want 
to consider discontinuous ansatz spaces due to the
different cases that can occur for face contributions. If you want or need 
to use more than one DoFHandler, this more or less the only case when you 
can't use MeshWorker directly. You might however, be able to combine 
different equations in a single DoFHandler.
What exactly is it that you want to do?

Best,
Daniel

[1] https://www.dealii.org/7.1.0/doxygen/deal.II/step_12.html 

Am Freitag, 26. August 2016 02:52:44 UTC+2 schrieb yocz...@gmail.com:
>
> Hi, 
> I'm reading the step-12, and i get confused that where and how should I 
> use the MeshWorker class? 
> The MeshWorker class seems just using in the DG examples. 
>
> But I do not think it is a good idea to use MeshWorker in my problems, the 
> problems, I will treat,  maybe a little  complicated, it combines 
> DG and couples different kinds of equations in different parts of the 
> domain. If I don't want to use the MeshWorker, how should I do to implement 
> DG?
>
> Are there any examples using DG without MeshWorker?
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: assembling individual block step 43

[Daniel Arndt]

Dear Franck,


> I am looking at step-43 using different method. I do a semi discretization 
> in space and end up with an ode for the saturation and use different 
> technique for time integration for the saturation.
>
> For this I will need to assemble the matrix H after solving for the 
> velocity in the darcy system.  Individual entries of H are given by:
>
> \mathbf{H}_{ij} & = - \Delta t^{(n)}_c \Big( F\left(S^{(n-1)}\right) 
> \mathbf{v}_i,\nabla\phi_j\Big)_{\Omega} 
>
>
> I am struggling with sparsity pattern that I have to define for this 
> matrix. I try this:
>
> assemble_matrix_H.clear();
> DynamicSparsityPattern dsp (n_s, n_u);
> DoFTools::make_sparsity_pattern (saturation_dof_handler, 
> darcy_dof_handler, dsp,
> saturation_constraints, false);
> ssemble_matrix_H.reinit(dsp);
>
>
>
> dealii is complaining and is saying
>
> /home/franckm/apps/candi/deal.II-toolchain/deal.II-v8.4.0/include/deal.II/dofs/dof_tools.h:563:3:
>  
> note: template void 
> dealii::DoFTools::make_sparsity_pattern(const DoFHandlerType&, const 
> DoFHandlerType&, SparsityPatternType&)
>make_sparsity_pattern (const DoFHandlerType &dof_row,
>^
> /home/franckm/apps/candi/deal.II-toolchain/deal.II-v8.4.0/include/deal.II/dofs/dof_tools.h:563:3:
>  
> note:   template argument deduction/substitution failed:
> /home/franckm/candi-examples/step-43b/step-43b.cc:918:68: note:   
> candidate expects 3 arguments, 5 provided
>saturation_constraints, false);
> ^
> make[3]: *** [CMakeFiles/step-43b.dir/step-43b.cc.o] Error 1
>
The error is not surprising as there isn't a variant of 
DoFTools::make_sparsity_pattern that takes two DoFHandlers.
 

> How one define a sparsity pattern for assembling H in the above scenario?
>
What you can do, is to consider all the blocks of your system and use 
\widetilde{H}=[0 H;0 0] instead. In this case your DoFHandler would have to 
have two blocks and 
you can use the variant [1] instead with an appropriate coupling.
A different approach is to build the SparsityPattern yourself by simulating 
assembling the matrix. This means that instead of using 
ConstraintMatrix.distribute_local_to_global
you would call SparsityPattern.add [2].

The question is if H really is the matrix you want to assemble. How does 
the equation you are considering really look like? 

Best,
Daniel
 
[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/group__constraints.html#ga33ec1ee2ac7c86b16700022b704b69c8
 
[2] 
https://www.dealii.org/8.4.0/doxygen/deal.II/classSparsityPattern.html#aa8dd1f7f393272cc5ae7789ac6be3413

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: where and how should I use the MeshWorker class? how should I do to implement DG without MeshWorker?

[Daniel Arndt]

Yongchao,

If you are still interested in using MeshWorker, there are some more 
examples (and in particular Stokes) at
https://bitbucket.org/guidokanschat/amandus

Best,
Daniel


Am Montag, 29. August 2016 11:15:22 UTC+2 schrieb yocz...@gmail.com:
>
> Dear Wolfgang ,
> Thanks for your advice, I have read the 6.2.1 version of step-12, and I 
> know how to deal DG without using MeshWorker.
> What's more there is the  6.2.1 version of step-12 without using 
> MeshWorker, if someone needs.
> https://www.dealii.org/6.2.1/doxygen/deal.II/step_12.html
>
>
> 在 2016年8月29日星期一 UTC+8上午11:38:34，Wolfgang Bangerth写道：
>>
>> On 08/27/2016 08:45 PM, yocz...@gmail.com wrote: 
>> > 
>> > Actually, I also want to use the MeshWorker,  what worried me is that I 
>> can't 
>> > control MeshWorker by myself, it something likes I know  the MeshWorker 
>> can do 
>> > lots of thing, but I don't know how to use MeshWorker. I am just 
>> beginning to 
>> > use deal.ii, and deal.ii seems a little complex to me.  So I want to 
>> find some 
>> > examples of MeshWorker.  And I am wondering if it is there any more 
>> advices? 
>>
>> Yongchao -- the issue with MeshWorker is that when it was implemented, we 
>> didn't write enough documentation/tutorial programs/examples that show 
>> how to 
>> use it in more complicated situations. Development then also stopped, and 
>> so 
>> there are valid cases where it would be nice to use MeshWorker for which 
>> nobody today knows any more whether they can or can not be done with 
>> MeshWorker. The only person who still knows how this all works 
>> unfortunately 
>> also doesn't regularly participate in the help forum any more :-( 
>>
>> In other words, as in other threads about MeshWorker on this forum, you 
>> can of 
>> course ask questions about how this or that can be done, but you're 
>> likely not 
>> going to get a lot of helpful answers because nobody knows any more how 
>> this 
>> works. It's really an unfortunate situation we got ourselves into with 
>> this 
>> small corner of the library. 
>>
>> 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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: assembling individual block step 43

[Daniel Arndt]

>
>
> There actually *is* a function that takes two DoFHandler objects, but it 
> doesn't take a constraint matrix or the last boolean argument.
>
[...] 

> You can of course see what variants of the make_sparsity_pattern function 
> there are by looking here: 
>https://www.dealii.org/8.4.1/doxygen/deal.II/namespaceDoFTools.html

Yes, I must have overlooked that... Well hidden from the other 
make_sparsity_pattern_functions.

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Neumann condtitions in the mixed space setting

[Daniel Arndt]

Eldar,

As you already found out, interpolate_boundary_values can only be used with 
a ComponentMask if the element you are using is primitive.
Therefore, the workaround you are using seems to be suitable for the moment.

Since BDM1 is div-conforming, project_boundary_values and 
project_boundary_values should have the same effect.

Best,
Daniel

Am Mittwoch, 31. August 2016 00:55:33 UTC+2 schrieb Eldar Khattatov:
>
> I managed to imply the Neumann conditions on stress only by first 
> projecting the boundary values to the entire mixed space (affecting 
> rotations), and then manually removing the constraints associated with 
> rotations. The code is below if someone is interested.
> VectorTools::project_boundary_values (dof_handler,
>   stress_boundary_functions,
>   QGauss(2),
>   boundary_values_stress);
>
> FEValuesExtractors::Scalar rotation(dim*dim + dim);
> VectorTools::interpolate_boundary_values (dof_handler,
>   0,
>   ZeroFunction(dim*dim + dim 
> + 0.5*dim*(dim-1)),
>   boundary_values_rotation,
>   fe.component_mask(rotation));
> 
> for (std::map::iterator it_r=
> boundary_values_rotation.begin();
>  it_r!=boundary_values_rotation.end();
>  ++it_r)
>boundary_values_stress.erase(it_r->first);
>
> MatrixTools::apply_boundary_values (boundary_values_stress,
> system_matrix,
> solution,
> system_rhs);
>
>
> I still wonder, though, if there is a nicer way to do it.
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: nearest neighbor

[Daniel Arndt]

Joel,

In GridTools::find_cells_adjacent_to_vertex(dof_handler,center_id) you need 
to give the number of the vertex_number as center_id, not the number of the 
degree of freedom.
Then, you want to use 

Quadrature  
quadrature_formula(dof_handler.get_fe().get_unit_support_points());

to create a Quadrature object that gives you the support points on each 
cell.
Finally, you can output the global dof number and their support point on 
each cell by: 

for (unsigned int j=0; jhttp://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] singularity error due to 1/0

[Daniel Arndt]

Kyusik,

Can you summarize what you have tried and which errors/problems you got 
with these approaches?
What exactly do you mean with "csc(theta) or sec(theta) can't be used in 
cylindrical coord system"?

Best,
Daniel

Am Freitag, 9. September 2016 06:17:44 UTC+2 schrieb hank...@gmail.com:
>
> Thanks for your answer. But, I still don't know how I can deal with this 
> error.
>
> Anyway, Thank you very much.
>
> Kyusik. 
>
> 2016년 9월 5일 월요일 오후 9시 47분 21초 UTC+9, Wolfgang Bangerth 님의 말:
>>
>> On 09/05/2016 06:04 AM, hank...@gmail.com wrote: 
>> > 
>> > K_Inv[0][0]=2*r/(cos(theta)*cos(theta)), 
>> > K_Inv[1][1]=2*r/(sin(theta)*sin(theta)), K_Inv[2][2]=r 
>> > 
>> > r is calculated by "r=sqrt(x^2+y^2) ", and theta is calculated by 
>> "theta=x/y" 
>> > 
>> > But, as you can expect that ,on the points where |x|<0.001 or 
>> > |y|<0.01,  cos(theta) or sin(theta) is almost zero. So, It seems 
>> that It 
>> > causes the singularity. 
>> > 
>> > So, At first I tried to change the above 2 element in K_Inv like 
>> this... 
>> > 
>> > K_Inv[0][0]=2*r/(cos(theta)*cos(theta)+del) 
>> >  K_Inv[1][1]=2*r/(sin(theta)*sin(theta)+del), 
>>
>> Instead of this approach, you should use the following: if r>delta, use 
>> the 
>> formula above. If r<=delta, use l'Hopital's rule to get an expression 
>> that 
>> avoids the division by zero. 
>>
>> That said, you still have a problem for those values of x,y where theta 
>> is 
>> close to zero or one, but r is not. For example, for x=0, y=1, you have 
>>r=1 
>>theta=pi/2 
>>cos(theta)=0 
>>K_inv[0,0] = 2/0 
>> You need to think about what you want to do with this situation. 
>>
>> 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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Applying Multi Point Constraints

[Daniel Arndt]

Rajat,

Since you are using p4est and PETSc, I suspect that you want to run with 
more than one MPI process.
If you are using a parallel::distributed::Triangulation you will probably 
run into trouble when the pair of degrees of freedom
you want to identify are not known by the same process.
You might want to have a look into 
https://groups.google.com/d/topic/dealii/aX41KqrWE1g/discussion
for how you can deal wit this situation by sending the required dofs to the 
respective MPI processes.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Using constraints for already assembled RHS

[Daniel Arndt]

Rajat,

The first equation is the usual equilibrium equation. It is solved to give 
> displacements. The stiffness matrix is then multiplied by this global 
> displacement vector to give the global traction vector which forms the rhs 
> for the second equation.
>
The easiest solution to deal with constraints to deal with them while 
assembling the matrix. Therefore, I would not compute the right-hand side 
for the second equation beforehand.
Instead just do it while assembling the second equation and use 
constraints.distribute_local_to_global().

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to find all the partial derivatives of the solution polynomial in deal.II

[Daniel Arndt]

Jiaqi,

Currently you can access first, second and third partial derivatives in 
FEValues. The FiniteElement class also implements fourth derivatives.
If this is sufficient for your purposes, then using 
FEValues*::third_derivative should be fine for you. I don't see a easier 
way at the moment.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Can we say that the higher order method, the more accurate?

[Daniel Arndt]

Jaekwank,

You would indeed expect significantly better results with higher polynomial 
degree on the same mesh, if the solution is sufficiently regular.
What kind of mesh are you using? Are you resolving all the boundary layers 
suitably? Are you using a Mapping of the same order as your velocity ansatz 
space is?

Best,
Daniel

Am Sonntag, 18. September 2016 17:57:16 UTC+2 schrieb JAEKWANG KIM:
>
>
> Hello, I am a starter of dealii and am learning a lot these days with the 
> help of video lectures and tutorial examples. 
>
> I modified step-22 code (stokes flow code) into my own problem, the flow 
> around sphere.
>
> and I intend to evaluate the drag force (which is analytically given by 
> stokes equation) 
>
> My code reached quite close to the value since the absolute error  : 
> abs(drag_calculated-drag_exact)/drag_exact is around 10^(-3)
>
> However, I expected that if I input higher 'degree' I will receive more 
> accurate result, but it didn't
>
> Obviously Q2 is better than Q1. and Q3 is better than Q2. But Q4 or Q4 is 
> not better than Q2 or Q3? 
>
> Is there any reason on this? 
>
> (To be specific, if i say degree 2 , that mean I use (2+1) for velocity, 
> (2) for pressure, and (2+2) for Gauss integral
>
>
> Thank you 
>
> Jaekwang Kim  
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: A*diag(V) with mmult?

[Daniel Arndt]

Franck,

One possibility is to use A.mmult(C,B,V) where B is an IdentityMatrix.
Another way is to create an IdentityMatrix B, modify its entries suitably 
and use A.mmult(C,B).
A third approach would be to just write the the scaling yourself.

Best,
Daniel

Am Sonntag, 18. September 2016 18:01:17 UTC+2 schrieb franck75:
>
> I am having a SparseMatrix A which has a certain sparsity pattern and a 
> vector V.
>
> How to perform the matrix multiplication 
>
> A*diag(V)
>
> where diag(V) is a diagonal matrix with the vector V in the main diagonal.
>
> is there any such function in dealii?
>
> how to create the SparseMatrix diag(V) for a given vector V?
>
> to my knowledge dealii on provided 
>
> A.mmult(C,B,V)  for   C= A*diag(V)*B or 
>
> A.mmult(C,B)  for   C= A*B  
>
> Cheers
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

> [...]
>  for (unsigned int q_index=0; q_indexfor (unsigned int i=0; i  {
> velocity[local_dof_indices[i]] = local_velocity_values[q_index];
>  }
>
> Here you assign the same Tensor<1,dim> at all DoFs multiple times and you 
end with velocity[*]=local_velocity_values[n_q_points]-1.
I don't think that this is what you want to do, is it? If I understand you 
correctly, you try to assign the magnitude of the velocity in each degree 
of freedom
to a different scalar Vector that represents a scalar FiniteElement Vector.
In this case, you need to use a Quadrature object that uses the 
unit_support_points of the scalar FE [1] and your assignment would read

for (unsigned int i=0; i >.

Best,
Daniel

[1] 
https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-to-get-the-mapped-position-of-support-points-of-my-element

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] VectorTools::integrate_difference function

[Daniel Arndt]

Jaekwang,

This seems to be unrelated. Did you setup the Triangulation before calling 
fe.distribute_dofs()?

Best,
Daniel

Am Dienstag, 20. September 2016 17:04:02 UTC+2 schrieb JAEKWANG KIM:
>
> yes...I tried that before, but I get an error message as follow.. 
>
> As you mentioned, I first learn this function in step 7 tutorial. 
> But the problem is that I want to integrate my vector-valued-solution to 
> exact-vector-valued solution, while step 7 compares scalar-solution. 
>
> Or is there any ways to compare each solution component separately? as it 
> did in step-7. 
>
> I am really thank you for your advice!
>
>
> An error occurred in line <235> of file 
>  in 
> function
>
> static types::global_dof_index 
> dealii::internal::DoFHandler::Policy::Implementation::distribute_dofs(const 
> types::global_dof_index, const types::subdomain_id, DoFHandler spacedim> &) [dim = 2, spacedim = 2]
>
> The violated condition was: 
>
> tria.n_levels() > 0
>
> The name and call sequence of the exception was:
>
> ExcMessage("Empty triangulation")
>
> Additional Information: 
>
> Empty triangulation
>
>
> Stacktrace:
>
> ---
>
> #0  2   libdeal_II.g.8.4.1.dylib0x000104ef307c 
> _ZN6dealii8internal10DoFHandler6Policy14Implementation15distribute_dofsILi2ELi2EEEjjjRNS_10DoFHandlerIXT_EXT0_EEE
>  
> + 300: 2   libdeal_II.g.8.4.1.dylib0x000104ef307c 
> _ZN6dealii8internal10DoFHandler6Policy14Implementation15distribute_dofsILi2ELi2EEEjjjRNS_10DoFHandlerIXT_EXT0_EEE
>  
>
> #1  3   libdeal_II.g.8.4.1.dylib0x000104ef2d97 
> _ZNK6dealii8internal10DoFHandler6Policy10SequentialILi2ELi2EE15distribute_dofsERNS_10DoFHandlerILi2ELi2EEERNS1_11NumberCacheE
>  
> + 39: 3   libdeal_II.g.8.4.1.dylib0x000104ef2d97 
> _ZNK6dealii8internal10DoFHandler6Policy10SequentialILi2ELi2EE15distribute_dofsERNS_10DoFHandlerILi2ELi2EEERNS1_11NumberCacheE
>  
>
> #2  4   libdeal_II.g.8.4.1.dylib0x000104eb6797 
> _ZN6dealii10DoFHandlerILi2ELi2EE15distribute_dofsERKNS_13FiniteElementILi2ELi2EEE
>  
> + 135: 4   libdeal_II.g.8.4.1.dylib0x000104eb6797 
> _ZN6dealii10DoFHandlerILi2ELi2EE15distribute_dofsERKNS_13FiniteElementILi2ELi2EEE
>  
>
> #3  5   mystokes0x000101350e7d 
> _ZN8MyStokes13StokesProblemILi2EE10setup_dofsEv + 397: 5   mystokes
> 0x000101350e7d 
> _ZN8MyStokes13StokesProblemILi2EE10setup_dofsEv 
>
> #4  6   mystokes0x000101341a52 
> _ZN8MyStokes13StokesProblemILi2EE3runEv + 386: 6   mystokes
> 0x000101341a52 _ZN8MyStokes13StokesProblemILi2EE3runEv 
>
> #5  7   mystokes0x000101340aa5 main + 69: 
> 7   mystokes0x000101340aa5 main 
>
> #6  8   libdyld.dylib   0x7fff95b575ad start + 1: 
> 8   libdyld.dylib   0x7fff95b575ad start 
>
> #7  9   ??? 0x0001 0x0 + 1: 9 
>   ??? 0x0001 0x0 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] VectorTools::integrate_difference function

[Daniel Arndt]

...of course I meant dof_handler.distribute_dofs(). Does the DoFhandler 
know at this point about the Trinagulation?

Best,
Daniel

Am Dienstag, 20. September 2016 17:58:02 UTC+2 schrieb Daniel Arndt:
>
> Jaekwang,
>
> This seems to be unrelated. Did you setup the Triangulation before calling 
> fe.distribute_dofs()?
>
> Best,
> Daniel
>
> Am Dienstag, 20. September 2016 17:04:02 UTC+2 schrieb JAEKWANG KIM:
>>
>> yes...I tried that before, but I get an error message as follow.. 
>>
>> As you mentioned, I first learn this function in step 7 tutorial. 
>> But the problem is that I want to integrate my vector-valued-solution to 
>> exact-vector-valued solution, while step 7 compares scalar-solution. 
>>
>> Or is there any ways to compare each solution component separately? as it 
>> did in step-7. 
>>
>> I am really thank you for your advice!
>>
>>
>> An error occurred in line <235> of file 
>>  in 
>> function
>>
>> static types::global_dof_index 
>> dealii::internal::DoFHandler::Policy::Implementation::distribute_dofs(const 
>> types::global_dof_index, const types::subdomain_id, DoFHandler> spacedim> &) [dim = 2, spacedim = 2]
>>
>> The violated condition was: 
>>
>> tria.n_levels() > 0
>>
>> The name and call sequence of the exception was:
>>
>> ExcMessage("Empty triangulation")
>>
>> Additional Information: 
>>
>> Empty triangulation
>>
>>
>> Stacktrace:
>>
>> ---
>>
>> #0  2   libdeal_II.g.8.4.1.dylib0x000104ef307c 
>> _ZN6dealii8internal10DoFHandler6Policy14Implementation15distribute_dofsILi2ELi2EEEjjjRNS_10DoFHandlerIXT_EXT0_EEE
>>  
>> + 300: 2   libdeal_II.g.8.4.1.dylib0x000104ef307c 
>> _ZN6dealii8internal10DoFHandler6Policy14Implementation15distribute_dofsILi2ELi2EEEjjjRNS_10DoFHandlerIXT_EXT0_EEE
>>  
>>
>> #1  3   libdeal_II.g.8.4.1.dylib0x000104ef2d97 
>> _ZNK6dealii8internal10DoFHandler6Policy10SequentialILi2ELi2EE15distribute_dofsERNS_10DoFHandlerILi2ELi2EEERNS1_11NumberCacheE
>>  
>> + 39: 3   libdeal_II.g.8.4.1.dylib0x000104ef2d97 
>> _ZNK6dealii8internal10DoFHandler6Policy10SequentialILi2ELi2EE15distribute_dofsERNS_10DoFHandlerILi2ELi2EEERNS1_11NumberCacheE
>>  
>>
>> #2  4   libdeal_II.g.8.4.1.dylib0x000104eb6797 
>> _ZN6dealii10DoFHandlerILi2ELi2EE15distribute_dofsERKNS_13FiniteElementILi2ELi2EEE
>>  
>> + 135: 4   libdeal_II.g.8.4.1.dylib0x000104eb6797 
>> _ZN6dealii10DoFHandlerILi2ELi2EE15distribute_dofsERKNS_13FiniteElementILi2ELi2EEE
>>  
>>
>> #3  5   mystokes0x000101350e7d 
>> _ZN8MyStokes13StokesProblemILi2EE10setup_dofsEv + 397: 5   mystokes
>> 0x000101350e7d 
>> _ZN8MyStokes13StokesProblemILi2EE10setup_dofsEv 
>>
>> #4  6   mystokes0x000101341a52 
>> _ZN8MyStokes13StokesProblemILi2EE3runEv + 386: 6   mystokes
>> 0x000101341a52 _ZN8MyStokes13StokesProblemILi2EE3runEv 
>>
>> #5  7   mystokes0x000101340aa5 main + 69: 
>> 7   mystokes0x000101340aa5 main 
>>
>> #6  8   libdyld.dylib   0x7fff95b575ad start + 1: 
>> 8   libdyld.dylib   0x7fff95b575ad start 
>>
>> #7  9   ??? 0x0001 0x0 + 1: 9 
>>   ??? 0x0001 0x0 
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Can we say that the higher order method, the more accurate?

[Daniel Arndt]

Jaekwang,

Have a look at step-10. There you see how the order of the mapping 
influences the error.
Basically, you have to set the mapping for your FEValues objects, when you 
compute error via integrate_difference
and when you project or interpolate some Function to an FE Vector.

Best,
Daniel

Am Dienstag, 20. September 2016 18:10:11 UTC+2 schrieb JAEKWANG KIM:
>
>
>
> 2016년 9월 18일 일요일 오전 10시 57분 16초 UTC-5, JAEKWANG KIM 님의 말:
>>
>>
>> Hello, I am a starter of dealii and am learning a lot these days with the 
>> help of video lectures and tutorial examples. 
>>
>> I modified step-22 code (stokes flow code) into my own problem, the flow 
>> around sphere.
>>
>> and I intend to evaluate the drag force (which is analytically given by 
>> stokes equation) 
>>
>> My code reached quite close to the value since the absolute error  : 
>> abs(drag_calculated-drag_exact)/drag_exact is around 10^(-3)
>>
>> However, I expected that if I input higher 'degree' I will receive more 
>> accurate result, but it didn't
>>
>> Obviously Q2 is better than Q1. and Q3 is better than Q2. But Q4 or Q4 is 
>> not better than Q2 or Q3? 
>>
>> Is there any reason on this? 
>>
>> (To be specific, if i say degree 2 , that mean I use (2+1) for velocity, 
>> (2) for pressure, and (2+2) for Gauss integral
>>
>>
>> Thank you 
>>
>> Jaekwang Kim  
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Input Q2 element generated by Gmsh using Gridin

[Daniel Arndt]

Chenchen,

We only support to read in meshes that are not refined. Translated to your 
situation this means that you are only able to use Q1 meshes as initial 
grid.
Of course, you can then use triangulation.set_manifold to define the 
correct behavior when you refine your mesh within dealii along with a 
Mapping of appropriate order.

Best,
Daniel

Am Mittwoch, 21. September 2016 17:15:05 UTC+2 schrieb Chenchen Liu:
>
> Hi all,
>
>
> I have to use quadratic shape function (Q2 elements) to do my project 
> (2D). Previously the code works well with Q1 linear elements, when I want 
> to input the gird with Q2 elements generated by Gmsh, it gives the 
> following error:
>
> An error occurred in line <1418> of file <../source/grid/grid_in.cc> in 
> function
>
> void dealii::GridIn<2, 2>::read_msh(std::istream &) [dim = 2, spacedim 
> = 2]
>
> The violated condition was: 
>
> false
>
> The name and call sequence of the exception was:
>
> ExcGmshUnsupportedGeometry(cell_type)
>
> Additional Information: 
>
> The Element Identifier <8> is not supported in the deal.II library when 
> reading meshes in 2 dimensions.
>
>
> And, it says the supported quad element should have 4 nodes and 4 edges. 
> Just want to check whether deal.ii support Q2 elements by Gmsh.
>
> My mesh looks like the following
>
>
> 
>
>
>
> If deal.ii does not support higher order element by Gmsh, what other 
> format is recommended? Thanks!
>
>
> Best,
>
> Chenchen
>
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: periodic boundary conditions

[Daniel Arndt]

Juan,

[...]
> 1) My first question is about how the "direction" works, is it 1 
> correspond to "x" and 2 corresponds to "y"?
>
No, direction 0/1/2 means that the support_points of the dofs to be 
identified just differ in the x-/y-/z-component.
What you are doing should be perfectly fine.
 

> 2) When I do this I obtain fairly reasonable solutions, however I am 
> getting a weird Dirichlet-type of boundary conditions on boundaries 1 and 3 
> (see the attached figures), that distort the fields.
> I seem unable to find what is wrong, so I wonder if anyone has experienced 
> this kind of behavior.
>
No, I haven't observed this problem. Can you check in your ConstraintMatrix 
that you don't use any other constraints? For doing so use 
ConstraintMatrix.print(std::cout)
and check if there is a DoF that is not constrained to another one.
 
Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: how to get 3D structured mesh from rotated 2D mesh in dealii

[Daniel Arndt]

Marek,


[...] I would like to process it in dealii by rotating it around y-axis, to 
> get 3D grid.
> Is it possible in dealii?
>
No, there is no such function yet.
 

> [...]
> My objective is the mesh in the file the attached file perf_chamber.png, 
> which was obtained by merging
> the inner cylinder and outer cylinder shell. The mesh seems however to 
> have some "double" vertices at
> the interface between the merged entities, which possibly cause, that the 
> computation fails, see image
> velocity.png
>
This is the approach I would have taken for creating the desired geometry. 
The problem here seems to be 
that the faces at the interface don't match as the inner cylinder has more 
cells. This is the reason why
GridGenerator::merge_triangulations doesn't work as expected.
How do you create the cylinder and the cylinder shell and how do they 
actually look like at the interface?

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: how to get 3D structured mesh from rotated 2D mesh in dealii

[Daniel Arndt]

Marek,

the attached file should be close to what you want to do.
Basically, you need to create two geometries: One for the cylinder shell 
and one for the cylinder.
Here, we take advantage of the fact that GridGenerator::cylinder always 
creates four subdivisions w.r.t height.
Furthermore, the cylinder needs to be refined once to match the vertices of 
the cylinder shell. Since we can only merge
coarse triangulations, we need to flatten the resulting triangulation.
If you want to create a gemoetry with a different height, you can delete 
cells in the cylinder triangulation using 
GridGenerator::create_triangulation_with_removed_cells
or add more cells by merging another cylinder.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.
Point<3> transformation2(const Point<3> &in)
{
  Point<3> out = in;
  double tmp=out(0);
  out(0)=out(2);
  out(2)=1-tmp;
  return out;
}

void grid_4()
{
  const CylinderBoundary<3> boundary (.5,2);
  Triangulation<3> tria1, tria;
  Triangulation<3> tria3, tria3_final;

  GridGenerator::cylinder_shell (tria1, 0.5, 0.5, 1, 8);
  GridGenerator::cylinder (tria3, 0.5, 1);
  GridTools::transform(transformation2, tria3);
  tria3.set_boundary (0, boundary);
  tria3.refine_global(1);
  GridGenerator::flatten_triangulation(tria3, tria3_final);

  std::ofstream out1 ("grid-4a.vtk");
  std::ofstream out3 ("grid-4b.vtk");

  GridOut grid_out;
  grid_out.write_vtk (tria1, out1);
  grid_out.write_vtk (tria3, out3);

  GridGenerator::merge_triangulations (tria1, tria3_final, tria);
  std::ofstream out ("grid-4.vtk");
  grid_out.write_vtk (tria, out);
}

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

Its been a week since my latest post, I was just wondering if you could 
> help me.
>
Did you try what I suggested? What were the problems?

Looking at the latest files you sent, local_velocity_values has still not 
the type std::vector >
and you are still using

for (unsigned int q_index=0; q_indexhttp://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: More power outages

[Daniel Arndt]

Everything should be available again.

Am Mittwoch, 28. September 2016 15:19:15 UTC+2 schrieb sotelo...@gmail.com:
>
> it is already 28th, the server still  down. is it going to take longer 
> than announced?
> Edith
>
> On Monday, September 12, 2016 at 11:21:13 AM UTC-5, Guido Kanschat wrote:
>>
>> Dear all,
>>
>> We are expecting more power outages on September 26 to 28. They are 
>> trying to fix the emergency supply such that the recent failures don't 
>> repeat. Let's keep our fingers crossed.
>>
>> The dealii web site will be down during these outages, which may last a 
>> few hours. Please be patient. 
>>
>> I apologize for the inconvenience,
>> Guido 
>>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

I did try it, but I cant just change the local_velocity_values from a 
> std::vector> to std::vector. Then the 
> fe_values[velocities].get_function_values function doesnt compile. 
>
That's true, but you should use std::vector > and not 
std::vector, see also [1].
 

> (correct me if wrong) This method of using norm_sqr() would only produce 
> the size, not the individual x,y and z components. I need the x,y and 
> z components separately so I can get the size but also the difference 
> between two solutions and that error of these two solutions. Hence my 
> previous message. I want to be able to manipulate the vector components. 
> How would I do this?
>
Yes, but since you get all the values at the local dofs, you can play 
around with them as you like and store the result in velocity. If this is 
not sufficient, you can also make velocity a std::vector and use something 
like
velocity[component][local_dof_indices[i]] = 
local_velocity_values[i][component]

Best,
Daniel

>
[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/classFEValuesBase.html#a396f1430dd3f5716a9fe6ef2762edb5d

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Neumann bc in step-20

[Daniel Arndt]

Erik,

Am Donnerstag, 29. September 2016 07:32:04 UTC+2 schrieb Erik Svensson:
>
> Hi, 
>
> I’m trying to implement homogeneous Neumann bc on part of the boundary in 
> tutorial example step-20. 
>
> The Neumann bc should be imposed strongly (as for Dirichlet bc in normal, 
> non-mixed, formulation). Is this correct?
>
> In any case, I’m trying to impose the Neumann bc strongly using the 
> technique from step-22.
>
What exactly are you doing? As far as I see step-22 just uses weak Neumann 
bc and if theses are homogeneous you normally don't have to implement 
anything additional.


> Running program I get the error:
>
> ---
> An error occurred in line <1878> of file 
> 
>  
> in function
> void 
> dealii::VectorTools::{anonymous}::do_interpolate_boundary_values(const 
> M_or_MC&, 
> const DoFHandlerType&, const typename dealii::FunctionMap space_dimension>::type&, std::map&, const 
> dealii::ComponentMask&, dealii::internal::int2type) [with 
> DoFHandlerType = dealii::DoFHandler<2>; M_or_MC = dealii::Mapping; int dim_ 
> = 2; typename dealii::FunctionMap::type = 
> std::map*, 
> std::less, std::allocator const dealii::Function<2, double>*> > >]
> The violated condition was: 
> cell->get_fe().is_primitive (i)
> The name and call sequence of the exception was:
> ExcMessage ("This function can only deal with requested boundary " 
> "values that correspond to primitive (scalar) base " "elements")
> Additional Information: 
> This function can only deal with requested boundary values that correspond 
> to primitive (scalar) base elements
> ---
>
For non-primitive FiniteElements (e.g. FE_RaviartThomas) you can't use 
VectorTools::interpolate_boundary_values. Use 
VectorTools::project_boundary_values [1] instead.

Best,
Daniel

[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/group__constraints.html#ga8981786b2369d1c6db8f5c363e290edc
 

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Vector-valued gradient of solution vector

[Daniel Arndt]

Joel,

Yes, this looks about right. I would really replace "std::vector double 
x,y,z" by "Vector x,y,z".
In the end, entries in these vectors only make sense with the corresponding 
DoFHandler.
This would also allow you to use these Vectors for postprocessing in 
deal.II.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: how to get 3D structured mesh from rotated 2D mesh in dealii

[Daniel Arndt]

Marek,

[...] 
>
I want to ask a question, which could be interesting for other users of 
> dealii-
> do we have some bulletproof procedure, how to produce (not so simple) 
> geometries,
> meshes in gmsh and then merge them "safely" in the dealii? Safely meaning- 
> without
> the creation of the doubled vertices at the boundary between merged 
> geoetries, meshes?
> Is there some procedure in dealii to cope with this problem?
> The function merge_triangulation requieres "matching" triangulation, could 
> it be possible
> to define a function, which would somehow merge vertices at the boundary 
> of two not-so-nice
> matching trianguations (possibly within some toleration, i.e. if two 
> vertices are less then \eps
> away from each other, than merge them)
> Does my idea make sense?
>
This is already done in GridTools::delete_duplicated_vertices [1]. The 
problem you experienced results from the fact that a coarse mesh cannot 
contain any hanging nodes in deal.II. Therefore, the faces at the interface 
have to match, i.e. there is 1-1 correspondence. This is something you have 
ensure when creating the meshes you want to merge
and I don't see an easy way to work around this automatically.

Best,
Daniel

[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/namespaceGridTools.html#a79820d66b08b06e94541817ba78d6213
 

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Derivatives (gradients) of Symmetric Tensors

[Daniel Arndt]

Metin,

[...]
> For the vectors, the implementation (for the else case) is as following;
>
>   *for* (*unsigned* *int* shape_function=0;
>
>shape_function
> {
>
> *for* (*unsigned* *int* d=0; d<*spacedim*; ++d)
>
>   *if* 
> (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
>
> {
>
>   *const* dealii::Tensor<1,*spacedim*> 
> *shape_gradient_ptr =
>
> &shape_gradients[shape_function_data[shape_function].
>
>  row_index[d]][0];
>
>   *for* (*unsigned* *int* q_point=0; 
> q_point
> divergences[q_point] += value * 
> (*shape_gradient_ptr++)[d];
>
> }
>
>   }
>
> To me, it's not clear how the "value and shape_gradient_ptr" is used or in 
> general how the (e.g.) u_x+v_y+w_z is done..
>
Let me just comment on this: 
The divergence at a quadrature point is decomposed into a sum over all 
shape functions and all of their components.
Since the value at a certain quadrature point q_p of a FE function u can be 
written as u(q_p)=\sum_i u_i \phi_i(q_p),
it holds \nabla \cdot u(q_p) = \sum_i u_i \nabla \cdot \phi_i(q_p) = \sum_i 
u_i \sum_c \partial_c \phi_i^c(q_p) where \phi_i^c is the cth component of 
the ith shape function. 
This is exactly what is done in the last line.
shape_gradient_ptr is initialized with the gradient of the dth component of 
the considered shape_function at the first quadrature point.
By "*shape_gradient_ptr++", the gradient at this quadrature point is 
returned and the pointer incremented such that it points to the
gradient at the next quadrature point.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: how to get 3D structured mesh from rotated 2D mesh in dealii

[Daniel Arndt]

Marek,


Only the vertices with indices in considered_vertices are tested for 
>> equality. This speeds up the algorithm, which is quadratic and thus quite 
>> slow to begin with. However, if you wish to consider all vertices, simply 
>> pass an empty vector.
>>
>
> i.e. I have assumed it is called by default for all vertices.
>
Even if it is called for all vertices, this would not help you with your 
problem. You would still have the problem with creating hanging nodes in 
your coarse mesh.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Derivatives (gradients) of Symmetric Tensors

[Daniel Arndt]

Metin,

I understand that the components are in space (x/y/z), 
> and (*shape_gradient_ptr++)[d] would be derivative of phi in each 
> direction; is that correct?
>
Yes, that is the derivative of the dth component in direction d.

>
> (compared to vectors) the implementation of divergence function for the 
> tensors is as following;
>
>   *if* (snc != -1)
>
> {
>
>   *const* *unsigned* *int* comp =
>
> 
> shape_function_data[shape_function].single_nonzero_component_index;
>
>
>   *const* dealii::Tensor < 1, *spacedim*> *shape_gradient_ptr 
> =
>
> &shape_gradients[snc][0];
>
>
>   *const* TableIndices<2> indices = dealii::Tensor<2,
> *spacedim*>::unrolled_to_component_indices(comp);
>
>   *const* *unsigned* *int* ii = indices[0];
>
>   *const* *unsigned* *int* jj = indices[1];
>
>
>   *for* (*unsigned* *int* q_point = 0; q_point < 
> n_quadrature_points;
>
>++q_point, ++shape_gradient_ptr)
>
> {
>
>   divergences[q_point][jj] += value * 
> (*shape_gradient_ptr)[ii];
>
> }
>
> }
>
This is similar to the previous one. Note that the ith component of the 
divergence of a rank-2-tensor is defined as sum_i \partial_i T_{ij}.
This is what is implemented for the general (non-symmetric) Tensor class. 
You have to compare this to the Vector implementation for primitive shape 
functions.
Basically, ii corresponds to the first index of the non-zero component and 
jj corresponds to the second index of the non-zero component.
We figure the correct component of the divergence add add the appropriate 
derivative to it.
For a SymmetricTensor, only T_{ij} for i\leq j is stored.Therefore, we need 
to have the line

if (ii != jj)
  divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii];

additionally to account for the contribution of T_{ji} add the same time.


> There is only the "snc != 1" case, and it seems that the component now 
> means the tensor components and not the directions in space. Is this 
> function (for tensors) not yet implemented? or is it sth related to 
> non-primitive shape functions (I would appreciate if you could comment on 
> this as well)?
>
This is only implemented for primitive shape functions.

Best,
Daniel 

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] how to use solution.

[Daniel Arndt]

Kyusik,

As Wolfgang already mentioned, you want to use 
FEValues::get_function_values() 
to extract from a global FE Vector the local values at the quadrature 
points in the cell you are on.
These quadrature points normally don't coincide with with the mapped unit 
support points of the
FiniteElement you are using.
Nevertheless, the global DoF indices of the local DoFs can be obtrained as 
usual by
cell->get_dof_indices()

Have a look at some of the time-dependen steps, e.g. 
step-35::assemble_one_cell_of_advection,
or just search for get_function_values in the examples folder.

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to solve for the null space of system matrix?

[Daniel Arndt]

Edith,

If I understand you correctly, you basically consider step-7 and replace 
the Neumann boundary conditions by inhomogeneous Dirichlet boundary 
conditions.
The system matrix is symmetric and positive definite and hence CG should 
work.

Can you verify that the solution is correct if you use a direct solver 
(SparseDirectUMFPACK [1])?

Best,
Daniel

[1] 
https://www.dealii.org/8.4.1/doxygen/deal.II/classSparseDirectUMFPACK.html

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Output parameters other than DoFs

[Daniel Arndt]

Lisa,


> [...]
> I do not know how to output below values:
>
> 1- u_r
> 2- sigma_xx, sigma_yy
> 3- sigma_rr , sigma_tt
>
For computing these values you should have a look at the DataPostProcessor 
class [1]. Its use is explained in step-29.
What you basically want to do for u_r in there is something like

double r = 
std::sqrt(evaluation_points[0]*evaluation_points[0]+evaluation_points[0]*evaluation_points[0])
computed_quantities[i](0) = uh[i](0)*evaluation_points[0]/r + 
uh[i](1)*evaluation_points[1]/r

and similar things for sigma_*.

Best,
Daniel


[1] https://www.dealii.org/8.4.1/doxygen/deal.II/classDataPostprocessor.html

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Daniel Arndt]

Hamed,


>>> It still works with Petsc and Trilinos is not implemented at all. I 
>>> appreciate it if you could let me know what changes should be made in 
>>> step-40 to make it work with Trilinos.
>>>
>> What exactly, do you mean by that? So you have both PETSc and Trilinos 
installed and step-40 works. Are you saying that PETSc objects are used 
even if you use 
"using namespace ::LinearAlgebraTrilinos;" ? How did you find out?

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Daniel Arndt]

In case there is an error with Trilinos, what is it?

Am Sonntag, 16. Oktober 2016 17:27:26 UTC+2 schrieb Daniel Arndt:
>
> Hamed,
>
>
>>>> It still works with Petsc and Trilinos is not implemented at all. I 
>>>> appreciate it if you could let me know what changes should be made in 
>>>> step-40 to make it work with Trilinos.
>>>>
>>> What exactly, do you mean by that? So you have both PETSc and Trilinos 
> installed and step-40 works. Are you saying that PETSc objects are used 
> even if you use 
> "using namespace ::LinearAlgebraTrilinos;" ? How did you find out?
>
> Best,
> Daniel
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Easy way to calculate second-deviatoric tensor

[Daniel Arndt]

Jaekwang,

>From this symmetric gradient tensor, do we have one easy way to calculate 
> the following "the second invariant of rate-of-strain tensor"? 
>
There is SymmetricTensor::second_invariant [1] which is defined as I2 = 
1/2[ (trace sigma)^2 - trace (sigma^2) ]
and SymmetricTensor::first_invariant [1] which just the trace I_1=trace 
sigma. 
According to Wikipedia [2], your \Pi_\gamma can then be expressed as 
\Pi_\gamma^2 = 1/2(I1^2 -2 I2).

Best,
Daniel

[1] 
https://www.dealii.org/8.4.0/doxygen/deal.II/classSymmetricTensor.html#a7bbdfc57da6931de6a1757a0fa7ee982
 

[2] https://en.wikipedia.org/wiki/Invariants_of_tensors

>
> 
>
>
>
> Regards, 
>
> Jaekwang Kim 
>  
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Daniel Arndt]

Hamed,

running your code with a Trilinos installation alone fails for me in line 
1617 with

An error occurred in line <279> of file 
 in function
void dealii::TrilinosWrappers::SolverBase::do_solve(const 
dealii::TrilinosWrappers::PreconditionBase&)
The violated condition was:
false
Additional information:
AztecOO::Iterate error code -2: numerical breakdown

and the same for no preconditioner. The same holds for PETSc. It seems that 
you should review if this matrix is assembled correctly and that CG is 
appropriate, i.e. is this matrix symmetric and positive definite?

I can't reproduce your error

"PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, 
probably memory access out of range"

running with a developer version.

Best,
Daniel


Am Montag, 17. Oktober 2016 04:47:32 UTC+2 schrieb Wolfgang Bangerth:
>
> On 10/16/2016 03:46 PM, Hamed Babaei wrote: 
> > Another point I need to mention is that I had the same problem when as 
> my 
> > first try of parallelization, I paralleled a much simpler code, the 
> step-25, 
> > based on the step-40 . 
> > I received the same Segmentation Violation error from Petsc. At that 
> time, I 
> > replaced PreconditionerAMG with PreconditionerJacobi and the problem 
> resolved. 
>
> This does not help right now, but as a general rule, it is far simpler to 
> debug things when you have a small, simple program that when you have a 
> large, 
> complicated one. In your case, you had a problem you didn't understand, 
> and 
> you chose to address it in a way that papered over it, rather than 
> properly 
> fixed it based on an understanding of what is going on. It is a truism 
> that 
> this sort of issue will come back at some time. 
>
> In other words, if you have a problem, debug it until you understand what 
> the 
> problem is, and then fix it the right way. Don't paper over it. 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

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

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Daniel Arndt]

Hamed,

which solvers do you think I can examin and see if they can be used instead 
> of SolverCG?
>

Running with the deal.II CG solver gives me


An error occurred in line <424> of file 
 in function
void 
PhaseField::Material_Compressible_Neo_Hook_Three_Field::update_material_data(const
 
dealii::Tensor<2, dim, double>&, const dealii::Tensor<2, dim, double>&, 
double, double, double, double, double) [with int dim = 3]
The violated condition was: 
det_F > 0

so there seem to be some more problems... 
Is everything working as expected with a direct solver (like 
SparseDirectUMFPACK)?

Best,
Daniel 

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Daniel Arndt]

Hamed,

At least I can't open your file in a way to clearly see what your 
bilinearform is, but could just guess from the implementation. Can you 
either write it here (using LaTeX or similar) or upload the file as pdf 
again?

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to use sacado package

[Daniel Arndt]

Jaekwang,

It seems that no fortran compiler is found. You should be able to configure 
Trilinos without fortran using
-DTrilinos_ENABLE_Fortran:BOOL=OFF

After that it's the usual make install [1].

Best,
Daniel

[1] https://www.dealii.org/developer/external-libs/trilinos.html

Am Samstag, 22. Oktober 2016 01:44:26 UTC+2 schrieb Jaekwang Kim:
>
>
> I'd like to use sacado package for automatic differentiation. 
>
>
> However, I have installed deal.ii on my computer with trilinos 
> configuration set off. 
>
>
> I looked back deal.ii installation manual and got to know that I have to 
> install trilnos package first. 
>
>
> However, when I try to install trilnos following instruction, 
>
>
> cd trilinos-12.4.2
> mkdir build
> cd build
>
> cmake\
> -DTrilinos_ENABLE_Amesos=ON  \
> -DTrilinos_ENABLE_Epetra=ON  \
> -DTrilinos_ENABLE_Ifpack=ON  \
> -DTrilinos_ENABLE_AztecOO=ON \
> -DTrilinos_ENABLE_Sacado=ON  \
> -DTrilinos_ENABLE_Teuchos=ON \
> -DTrilinos_ENABLE_MueLu=ON   \
> -DTrilinos_ENABLE_ML=ON  \
> -DTrilinos_VERBOSE_CONFIGURE=OFF \
> -DTPL_ENABLE_MPI=ON  \
> -DBUILD_SHARED_LIBS=ON   \
> -DCMAKE_VERBOSE_MAKEFILE=OFF \
> -DCMAKE_BUILD_TYPE=RELEASE   \
> -DCMAKE_INSTALL_PREFIX:PATH=$HOME/share/trilinos \
> ../
>
> make install
>
>
> I got error message as follows...
>
>
>
>
> "
>
> The Fortran compiler identification is unknown
>
> CMake Error at 
> cmake/tribits/core/package_arch/TribitsGlobalMacros.cmake:1638 
> (ENABLE_LANGUAGE):
>
>   No CMAKE_Fortran_COMPILER could be found.
>
>
>   Tell CMake where to find the compiler by setting either the environment
>
>   variable "FC" or the CMake cache entry CMAKE_Fortran_COMPILER to the full
>
>   path to the compiler, or to the compiler name if it is in the PATH.
>
> Call Stack (most recent call first):
>
>   cmake/tribits/core/package_arch/TribitsProjectImpl.cmake:188 
> (TRIBITS_SETUP_ENV)
>
>   cmake/tribits/core/package_arch/TribitsProject.cmake:93 
> (TRIBITS_PROJECT_IMPL)
>
>   CMakeLists.txt:89 (TRIBITS_PROJECT)
>
> "
>
>
> What's the problem? 
>
> and I would very thank if someone tell me how to change deal.ii 
> configuration to set on Trilinos after that
>
>
> Thank you! 
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
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.
For more options, visit https://groups.google.com/d/optout.