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

[Daniel Arndt]

Hamed, 

[...]
> I am going to do the same as Bastian, namely periodic boundary condition 
> for displacement  such that u( 0, y ) = u( 1, y ) + *lambda ,*
> * so I used the recommended above code by Bastian for applying 
> inhomogenity to predefined periodicity, in order to have reletive 
> dispacement between right and left faces,*
> *but my results shows that both left and right faces are moving equally to 
> the same direction, not relative to each other.*
>
What you are doing is
dof_left = dof_right
dof_right = c
resulting in
dof_left = dof_right = c.
The solution is to set inhomogeneities for the dofs on the left boundary 
instead of the right boundary and thus
dof_left = dof_right +c.
In particular, you don't want to constraint dofs that are not already in 
the ConstraintMatrix after you called make_periodicity_constraints.
Hence, the code should then work without calling constraints.add_line().

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 Condition

[Daniel Arndt]

Hamed,

As you can see my assembly in the first post, my system_matrix and 
> system_rhs are not made directly from cell_system_matrix or 
> cell_system_rhs, but they are the result of some manipulations on 
> mass_matrix, laplase_matrix and nl_matrix (the same approach as step_25) 
> that have already been assembled (in the red part of the code). 
>
> To simplify the question, lets say we want to apply any sort of 
> constraints in the assembly part of step-25. I was wondering how to do use 
> ConstraintMatrix::distribute_local_to_global, while there is no 
> cell_system_matrix and cell_system_rhs.
>
This is only (directly) possible if you are not using inhomogeneous 
constraints. In that case, you can call 
ConstraintMatrix::distribute_local_to_global separately for matrix and 
right-hand side. 
In case you have inhomogeneous constraints, you can first create a vector 
that has the appropriate constraints:
inhom_constraints.distribute(u_inhom)
Afterwards you are just solving for the difference (u-u_inhom) that takes 
homogeneous constraints:
A*(u-u_inhom) = f-A*u_inhom
u=(u-u_inhom)+u_inhom

This approach should also work if A and f are linear combination of 
matrices resp. right-hand sides for which you used 
ConstraintMatrix::distribute_local_to_global individually.

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 condition in Meshworker.

[Daniel Arndt]

Sudarshan,

J-P is right. You can use constraints within the MeshWorker framework and 
hence also periodic boundary conditions.
In more detail, you have to provide the assembler with the 
ConstraintMatrix. In addition to step-45 for periodic boundary conditions,
you should also have a look to step-16 
[1] (LaplaceProblem::assemble_system ()) for using constraints with 
MeshWorker.

Best,
Daniel

[1] http://www.dealii.org/8.4.1/doxygen/deal.II/step_16.html

Am Mittwoch, 16. November 2016 22:43:55 UTC+1 schrieb Jean-Paul Pelteret:
>
> Sudarshan,
>
> Does MeshWorker take a constraints matrix? I suspect that it does 
> .
>  
> If so then yes, you should be able to implement periodic boundary 
> conditions.  I refer you to step-45 
>  to show 
> how this is done.
>
> Regards,
> Jean-Paul
>
> On Wednesday, November 16, 2016 at 9:43:23 PM UTC+1, Sudarshan Kumar wrote:
>>
>>
>> Is there any way to implement periodic  boundary condition with  
>> MeshWorker,?
>>
>> Any help is appreciated.
>>
>> Thanks a lot, 
>> Best
>>
>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"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: hp HDG in deal.II

[Daniel Arndt]

Samuel,

I am attempting to implement a hp hybridized dg-method for solving 
> ellipctic problems in deal.II. I am relying on steps 27 and 51 of the 
> deal.II tutorial programs. However, when I try to group some FE_FaceQ 
> elements in an hp::FeCollection and distribute the degrees of freedom, I 
> get an ExcNotImplemented() exception thrown by some method 
> hp_vertex_dof_identities. Am I getting something wrong? Or is this feature 
> just not implemented for
> FaceQ yet?
>
This feature is just not implemented yet. As FE_FaceQ is in fact 
discontinuous, you can try to just copy
hp_vertex_dof_identities, hp_line_dof_identities and hp_quad_dof_identities 
from source/fe/fe_dgq.cc to source/fe/fe_faceq.cc.
If you are successful with this, we would gladly accept a patch.  :-)

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 Condition

[Daniel Arndt]

Hamed,

Am Dienstag, 15. November 2016 19:09:18 UTC+1 schrieb Hamed Babaei:
>
> Hi all,
>
> Following step-45, I am going to implement periodic boundary condition on 
> my code for a Thermo-elastic problem, in which the thermal and elastic 
> parts are solved separately (I have two different dof_handelers for 
> temprature and displacement fields without any block structure). I am 
> wondering how to apply periodic constrains, made in setup_system, in the 
> thermal field assemble_system for which I don't use 
> "constraints.distribute_local_to_global ()" function. Thermal 
> Assemble_system is the following: 
>
[...]
> 
> for (unsigned int i=0; i   {
> for (unsigned int j=0; j   {
>  mass_matrix.  add (local_dof_indices[i],
>  local_dof_indices[j],
>  cell_mass_matrix(i,j));
>
>laplace_matrix.add (local_dof_indices[i],
>local_dof_indices[j],
>  cell_laplace_matrix(i,j));
>
>nl_matrix. add (local_dof_indices[i],
>local_dof_indices[j],
>cell_nl_matrix(i,j));
>
>   }
> nl_term (local_dof_indices[i]) += cell_nl_term(i);
>   }
> I guessed perhaps I can replace the red part with the following, 
>
>  constraints.distribute_local_to_global (cell_mass_matrix,
>   
>  local_dof_indices,
>mass_matrix);
> constraints.distribute_local_to_global (cell_laplace_matrix,
>   
> local_dof_indices,
>laplace_matrix);
> constraints.distribute_local_to_global (cell_nl_matrix,
>   
> local_dof_indices,
>   nl_matrix);
> constraints.distribute_local_to_global (cell_nl_term,
> 
>  local_dof_indices,
>   nl_term);
>
Can you confirm that these approaches are the same if no constraints are 
specified? This should definitely hold.
 

>
> but the newton_rophson solver for thermal part diverged when I did so!!
>
It is difficult to help you with just this information. How do the 
constraints you are using actually look like?
Are they correct?

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: Very basic question on setting boundary id

[Daniel Arndt]

Jaekwang,

It seems that you are missing
constraints.distribute(solution);
after solving the linear system. The constrained dofs are set to zero 
during assembly and 
this call is necessary to get the correct values for these in the solution 
as well. You might want to have a look at
https://www.dealii.org/8.4.1/doxygen/deal.II/group__constraints.html or 
step-6.

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: storing sparse matrix in .txt and scan in Matlab

[Daniel Arndt]

Anup,

 
>
>> 1. if the way of writing the sparse matrix on a .txt file is appropriate 
>> for post-processing in matlab, if not, please suggest
>> a way;
>>
>
> There are three print functions in the deal.II SparseMatrix. As explained 
> in the documentation 
> ,
>  
> you've chosen one that dumps the entire matrix with all zero entries not in 
> the sparsity pattern printed as well. This function 
> 
>  
> is probably the one that you want.
>
In addition to what J-P is saying, you can use 
SparseMatrix::print_formatted providing the zero string "0." instead of "" 
for exporting to some file "matrix.txt" and then load "matrix.txt" in 
Matlab.
  

>  
>
>> 2.  how should I make matlab understand the size of the sparse matrix and 
>> arrangement of the elements.
>>
> Using SparseMatrix::print you get a list of entries starting at zero 
instead of one. You can use the "Import Data" feature with suitable 
delimiters "( " and ")"
to create a matrix in Matlab,
load "matrix.txt"
shift the first two columns by 1 via
matrix(:,1)=matrix(:,1)+1
matrix(:,2)=matrix(:,2)+1
and convert this to a sparse matrix via
spconvert(matrix)

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: Convergence problem arising for large number of DoFs

[Daniel Arndt]

Hamed,

Using the print function of Sparsematrix class (Thanks to Daniel for 
> letting me know that) I printed the elements of system_matrix for both 
> sequential and parallel codes. 
> I reduced the problem to only 54 DoFs. It seems that both system_matrices 
> are symmetric and Identical except for some components which despite 
> others, are very small numbers i.e. O(1e-14 or less).
>
Then I would expect that the solver should behave the same for both 
matrices. Are you still running into the same problems using just 4 cells 
with your parallel code?

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: Convergence problem arising for large number of DoFs

[Daniel Arndt]

Hamed,

Thanks for your help. I can reduce my problem to four 2D elements. I would 
> like to plot the system_matrix of sequential and parallel codes to compare 
> them but I am not sure how to do so.
> Can I loop over the elements of system_matrix and simply plot them?
>
All of the matrix classes also have a print method, see [1] for Trilinos.

Best,
Daniel

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

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"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: Visit error...?

[Daniel Arndt]

Kyusik,



Am Donnerstag, 3. November 2016 10:51:20 UTC+1 schrieb hank...@gmail.com:
>
> Hi all,
>
> I'm trying to calculate absolute values of (exact_solution - 
> numeric_solution) and visualize it using Visit.
>
> So, I added the following in my code.
>
> error = solution;
> error -= exact_solution;
>
> for(unsigned it i=0;i error(i)=std::fabs(error(i));
>
> ...
> data_out.add_data_vector(solution,"solution");
> data_out.add_data_vector(exact_solution,"exact_solution");
> data_out.add_data_vector(error,"error");
>
Yes, this gives you a representation of the discretization error.

And, I also calculated the {l2_norm of (exact_solution - numeric_solution) 
> / sqrt(#cells)} (I divided l2norm by sqrt(#cell) Because I want to know 
> kind of average value of difference between two solution.)
>
> i.e. I want to calculate  sqrt{sum_i(exact_sol_i - sol_i)^2 / (#cell)} (I 
> think it is kind of average value of difference between two solutions)
>
What do you mean by "average value" here? You probably just wnat to compute 
the error in the L2-norm, isn't it?


> So, I added the following in my code
>
> Vector difference_per_cell (triangulation.n_active_cells());
> VectorTools::integrate_difference(dof_handler,
>   solution,
>   ExactSolution(),
>   difference_per_cell,
>   QGauss(2),
>   
> VectorTools::L2_norm);
> const double L2_error = difference_per_cell.l2_norm();
>
> std::cout<<"l2norm(error)="<
> So, l2_norm(error) is 1.14461e-09 
>
> But the almost all values in the plot of error(Error_fe1.png) are about 
> 10^-4 that is 10^5 times higher than l2norm.
>
If the error is O(1.e-4) you expect the L^2 error to be approximately 
\sqrt{\int_\Omega 1.e-8} = |\Omega|^{1/2} 1.e-4.
Assuming |\Omega|~1, this would suggest that you use approximately 
1.e-4/1.e-9=1.e5 cells.

What I want to say: Consider just the L2-norm! Due to |\Omega|~1, this 
should give you an error in the magnitude of what you are observing in the 
plot.

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: Extra layer around mesh

[Daniel Arndt]

Joel,

Same error with your changes. I searched for subdivided in 
> http://www.dealii.org/developer/doxygen/deal.II/changes_between_8_3_and_8_4.html
>  
> and found under specific improvements, that subdivided_parallelepiped 
> produced the wrong boundary indicators. Could this be the reason, why it 
> doesnt work? Perhaps I need to update the version Im running.
>
Yes, this makes perfect sense. Either you update or you set the boundary 
indicators correctly yourself.

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: hp convergence for non-linear solver.

[Daniel Arndt]

Jaekwang, 

[...]
>
> VectorTools::integrate_difference 
> (dof_handler,solution,Solution(),
>
>difference_per_cell
> *,QGauss(degree+**2**)*,VectorTools::L2_norm);
>
somegthing that immediately comes to mind is that you don't use a 
non-default Mapping in here although your boundary is curved.
Make sure that you use an appropriate Mapping everywhere!

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] computing some other output quantity after solution is calculated.

[Daniel Arndt]

Kyusik,

[...]
> 3. I calculated l2 norm 
>  Difference = solution;
>  Difference -= exact_solution;
>  norm = Difference.l2norm();
>
> But, according to what you mentioned, I calculated l2norm of U_i-U_i_exact 
> not l2norm of u_h(x)-u_h_exact(x)
>
Yes, this gives you the l2-norm of U_i-U_i_exact and not the L2-norm of 
u_h(x)-u_h_exact(x).
VectorTools::integrate_difference[1] allows you to compute the L2-norm. 
This function returns a Vector containing the error on each cell.
The l2-norm of that Vector is the L2-norm you are looking for. Have also a 
look at step-7.

Best,
Daniel

[1] https://www.dealii.org/8.4.1/doxygen/deal.II/namespaceVectorTools.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: computing some other output quantity after solution is calculated.

[Daniel Arndt]

Kyusik,

 //by Han(To get psi0) 
>   Vector::iterator max;
>   max=std::max_element(solution.begin(),solution.end());
>
This can only work if you know that the maximum is attained at one of the 
support_points
of your FiniteElement. In general, you can use 
VectorTools::integrate_difference to find a good approximation
to the maximum value of your solution. 

cell_rhs(i) += (fe_values.shape_value(i, q_index)*
> right_hand_side.value(fe_values.quadrature_point(q_index),
> sol_tmp[q_index],*max) *
> fe_values.JxW (q_index));
>
I don't quite understand this line. Does right_hand_side.value really take 
three arguments? 
How does this resemble x*solution+solution_max/y?


> In other words, I just calculated  (pi_i,pi_j) x J_tor_j=(rhs_i,pi_i)   
>  (where pi is test function)
> It worked, the results are as I expected. But someone said to me it is 
> very inefficient way.
> I agree with that opinion
>
In principle, your approach looks good. You can have a look at the 
implementation of VectorTools::project for improving this.
Why would you think that what you are doing is extremely inefficient?

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: boundary conditions

[Daniel Arndt]

>
> so the uniform flux from left and right of the rectangle implies periodic 
> boundary condition. But the K ( hydraulic conductivity) is a function of 
> (x,y).
> If we want to enforce the periodic boundary condition, should we expect to 
> have a condition on K ? should K be periodic as well for example ?
>
If you want to use periodic boundary conditions, then you should prescribe 
them and not Neumann boundary conditions. By having something like
n \cdot \nabla H|_\Gamma_{in} = -n \cdot \nabla H|_\Gamma_{out},
you are not constraining the tangential direction. 
For periodic boundary conditions, K|_\Gamma{in}=K|_\Gamma{out} makes sense. 
Probably you don't want to have a jump in K across the periodic boundary 
faces.
You might want to have a look at step-45[1] for how to describe periodic 
boundary conditions in deal.II.


> for example, a specific distribution of K field, where  we imagine a block 
> of material near the right edge at the exit nodes with much higher 
> hydraulic conductivity, then the flux out of the rectangle can not be 
> uniform. the flow field will be distorted and will become non uniform. 
> enforcing a uniform flow works in an artificial sense in this case. 
> So it seems to me that the limitation on the specification of a fixed head 
> is more essential than just making the system unknown up to a fixed 
> constant.
>
I am quite not sure what the question is you want to be answered in the 
end. 

Best,
Daniel 
 
[1] https://www.dealii.org/8.4.1/doxygen/deal.II/step_45.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: Error with the content "dealii::numbers::is_finite(sum)"

[Daniel Arndt]

Benhour,

are you running in serial or with more than one process? In the first case, 
I would suspect that your matrix is singular.

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: boundary conditions

[Daniel Arndt]

OK. So do I understand you correctly, that you are imposing Neumann 
boundary conditions on all *four* boundaries?
In that case the solution is only unique up to a constant.

It seems to my that in the problem I posed, the boundary condition is not 
> correct.
>
Can you be a bit more specific? How does the weak formulation you consider 
look like?
What kind of problems do you observe?

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: boundary conditions

[Daniel Arndt]

Retired Replicant,

H can be hydraulic head and K is the hydraulic conductivity. 
> Or H can be temperature with k being diffusion coefficient.
>
Have a look at step-7, if you just consider a diffusion problem. There 
Neumann boundary conditions are used on part of the boundary which is also 
mathematically totally fine.

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: Trouble with VisIt 2.11

[Daniel Arndt]

Praveen,

It's hard to tell from your description what is causing the error. 
At least, Visit has some internal problems and if you can open with a 
previous version, chances are high that this is a Visit issue. Can you open 
the file with Paraview or any other suitable program? Can you provide us 
with the file or create a minimal example?

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.

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.

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,

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.

[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] 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 
> =
>
> _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] 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 =
>
> _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.

[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.

[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: 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,

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_index

[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.

[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: 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: 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.

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

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: 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.

[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: 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: 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: 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.

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.

[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; j

[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: 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: 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<dim, spacedim>::get_function_gradients(const 
>> InputVector&, std::vector<dealii::Tensor<1, spacedim, typename 
>> 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<std::vector<Tensor<1,dim>> 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.

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: [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.

[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: 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: 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: 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.

Re: [dealii-developers] step 47

[Daniel Arndt]

You might want to have a look at
http://journals.ub.uni-heidelberg.de/index.php/ans/article/view/22317/23996
This seems to be quite similar.

Best,
Daniel

Am Dienstag, 9. August 2016 19:28:51 UTC+2 schrieb bangerth:
>
> On 08/08/2016 11:55 PM, sella...@gmail.com  wrote: 
> > I am a new user of dealii and I came across the step-47 program which 
> implements level set method. But I could not find its complete 
> documentation anywhere. 
> > 
> > It would be very helpful if you could share the document explaining the 
> math behind this step. 
>
> There is none :-( The program was started in 2012, but we never finished 
> it, and nothing was ever written up. 
>
> I had always hoped that we would eventually finish it, but I don't see 
> that happening. It may be best to just delete it. 
>
> Best 
>   Wolfgang 
>

-- 
You received this message because you are subscribed to the Google Groups 
"deal.II developers" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii-developers+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Undefined reference to `tbb::interface5::internal::task_base::destroy(tbb::task&)

[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.

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.

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: 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: 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: 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.

[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 >([0], [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 

[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: 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: 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.

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: 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.

[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: 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.

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: 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.

[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 Vector(3));
>
> std::vector Vector(3));
>
> std::vector 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]

ot 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<parameters.local_refinement_cycles; 
>>> ++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<GridTools::PeriodicFacePair>> 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'

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: 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.

[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: 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: 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.

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.