Re: [deal.II] Re: deal.ii and Trilinos

[Wolfgang Bangerth]

On 06/08/2016 09:26 AM, Marco Delchini wrote:


Our objective is to solve for two-phase flows in multi-D using a fully
implicit temporal scheme. In the short term, we would need the Trilinos
solver/perconditioner capabilities and also are very interested in the Sacado
package that does automatic differentiation.


As already mentioned, all of these have deal.II wrappers (except Sacado, which 
doesn't need that).


You may want to read through step-33, which does all of the things you want to 
do for a different equation.


Best
 W.

--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Transfer of quadrature point history data after mesh refinement.

[Wolfgang Bangerth]

On 06/14/2016 10:44 AM, Denis Davydov wrote:

As for the solution transfers, i would speculate that you need two
different ones: one for the main field you solve for and another
for quadrature data.


Correct. You need to define a separate DoFHandler for the DG field to 
which you project your quadrature point data. Then you need a separate 
(second) SolutionTransfer object for the second DoFHandler.


Best
 W.
--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Fast FE evaluation over several points and several functions

[Wolfgang Bangerth]

I am solving an eigenvalue problem similar than step-36. After solving for the
eigenpairs, I evaluate the eigenfunctions in the standard way:

 VectorTools::point_value ( mapping, dof_handler, efun[m],
q_points[j], Uq );

where: efun[m] is the m-th eigenfunction from step-36,
q_points[j] are selected quadrature points,
Uq is where I store the FE evaluation.

As expected, this operation is awfully slow! it takes seconds for a single
point evaluation with a decent discretization and having m eigenfunctions
makes it worse!


Yes. But there is a much more efficient way to do this:
  FEValues fe_values (...);
  std::vector sol_at_q_points (...);
  for (cell=...)
{
   fe_values.reinit (cell);
   fe_values.get_function_values (efun[m], sol_at_q_points);

This gives you the values of efun[m] at all quadrature points on the current 
cell at once, and this approach is efficient and independent of the overall 
size of the problem.


Best
 W.

--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Adding dirichlet values to ConstraintMatrix

[Wolfgang Bangerth]

On 06/23/2016 02:46 AM, Praveen C wrote:


void add_dirichlet_constraints (const
std::map<types::global_dof_index,double> ,
ConstraintMatrix
 )
{
   for (const auto  : values)
   {
  constraints.add_line (pair.first);
  constraints.set_inhomogeneity (pair.first, pair.second);
   }
}



Yes, this would work. I would add an assertion that pair.first is not already 
constrained. You'll also have to avoid using the C++11-style for loop. Other 
than that, this looks good!


Best
 W.

--
----
Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] On the use of MappingFEField

[Wolfgang Bangerth]

Praveen,


I change to

MappingFEField<dim,dim,TrilinosWrappers::MPI::Vector> map(dh_euler, 
euler_vector);

But I get a linking error

Undefined symbols for architecture x86_64:
   "dealii::MappingFEField<2, 2, dealii::TrilinosWrappers::MPI::Vector,
dealii::DoFHandler<2, 2> >::MappingFEField(dealii::DoFHandler<2, 2> const&,
dealii::TrilinosWrappers::MPI::Vector const&, dealii::ComponentMask)",
referenced from:
   _main in main.cc.o


Indeed. You need this patch:
  https://github.com/dealii/dealii/pull/2717
Best
 W.

--
----
Wolfgang Bangerth   email:bange...@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.

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

[Wolfgang Bangerth]

Thank you for your suggestions and sorry about late reply. I decided to
simplify the starting mesh and display the resulting matrix using
*MatrixOut*, so that I can check whether I went wrong in assembling the
matrix. But after using the following code to get the output matrix
data file, I don't know how to display the the matrix in gnuplot with their
values showing in the corresponding positions.

|
FullMatrix
<http://dealii.org/8.4.1/doxygen/deal.II/classFullMatrix.html>M;

...// fill matrix M with some values
// now write out M:

MatrixOut 
<http://dealii.org/8.4.1/doxygen/deal.II/classMatrixOut.html>matrix_out;

std::ofstream out("M.gnuplot");
matrix_out.build_patches
<http://dealii.org/8.4.1/doxygen/deal.II/classMatrixOut.html#a0d01113fd37dfa352dc237f3d05c2182>(M,"M");
matrix_out.write_gnuplot
<http://dealii.org/8.4.1/doxygen/deal.II/classDataOutInterface.html#a471d836d307c6e888ab3b5ad81fe0416>(out);

|

I tried to use

|
plot 'M.gnuplot'matrix
|

but all I get is some point without values showing up. Opening this file with
Emacs, the values can be found, but their orderings are not convenient.
Do you know how to display this data file so that they can show the values in
correct order just like in MATLAB.


You can do this in gnuplot:
  set style data lines
  set hidden3d
  splot "M.gnuplot" using 1:2:3

Best
 W.

--
----
Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Segmentation fault in WorkStream::run() for colored iterators in parallel(MPI)

[Wolfgang Bangerth]

Julius,


the current MeshWorker::loop() uses the WorkStream::run() which copies in 
serial.
I tried to use the WorkStream::run() version which uses colored iterators and
therefore assembles in parallel, in the sense of a tbb::parallel_for().

My code runs without using MPI as well as using an arbitrary amount of MPI
procs with a single thread each. If I use 2 MPI procs with 2 threads each, it
crashes due to a segmentation fault.
For the construction of the colored iterators I am using the
make_graph_coloring() in such a way that two cells sharing a face have
different colors.

The segmentation fault always appears while calling assemble() of the
ResidualSimple in simple.h.

Does anybody have a clue why I get a segmentation fault after adding local
data to the global dst vector inside assemble(), although I avoid write
conflicts by colored iterators?


Nothing particularly obvious comes to mind, without having tried to debug this 
in more detail.


Have you tried running the program under a debugger? Segmentation faults are 
often relatively easy to debug because something is accessing an invalid 
address, and the debugger can tell you what object this address corresponds 
to, and how this address got into the current function. Try the approach 
outlined here, for example:


https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#how-do-i-debug-mpi-programs

using the `mpirun -np 4 xterm -e gdb ./my_executable` approach.

Best
 W.

--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] question about compute_projection_from_quadrature_points_matrix function

[Wolfgang Bangerth]

On 07/26/2016 06:21 AM, dealii.gr...@gmail.com wrote:


I use a discontinuous field to transfer quadrature points data during mesh
refinement. In addition I use same method to write quadrature data into output
file. Why compute_projection_from_quadrature_points_matrix function sometime
creates negative data for a quadrature data which is always positive (such as
effective stress or equivalent plastic strain)?


Because it does a polynomial interpolation. Imagine you are trying to put a 
(quadratic) parabola through the points

  x=0,   y=0
  x=1/2, y=0
  x=1,   y=1
Here, all data values (y) are non-negative. But if you put a parabola through 
all data points, the parabola takes on negative values between x=0 and x=1/2.


You can generate similar examples also with linear elements if you choose the 
x- and y-values appropriately.


Best
 W.

--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Problem using WorkStream with a hp::DoFHandler

[Wolfgang Bangerth]

On 07/26/2016 08:54 AM, Claire wrote:


I hope that I provide enough information.


Bruno suggestion was a good guess, but not the cause. The real cause is that 
your local assemble function takes the following arguments:

  const cell_iterator&,
  Assembly::ScratchData<3>&,
  Assembly::CopyData&
But, because you only work on active cells, the compiler is confused that you 
are saying in the first call to WorkStream::run() that you want to use active 
iterators, but that the local assembly function takes a (not necessarily 
active) iterator argument.


Change the first argument of the local assembly function to 
active_cell_iterator and it should all work.


Best
 W.

--
----
Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Issue on Windows x64 Debug linkage

[Wolfgang Bangerth]

On 07/27/2016 08:05 PM, Tulio Ligneul wrote:

In my system, deal_II.lib is about 400MB and deal_II.g.lib is almost 2GB.


How much memory/RAM do you have in your system? I would not be surprised if 
you needed 10GB to link a 2GB library. If you don't have that much, the 
operating system will start to swap out data and then it can take forever.




Maybe the issue is somehow related to these reports:


https://connect.microsoft.com/VisualStudio/feedback/details/995447/linker-is-hanging-while-linking-some-project-in-debug-x64-configuration

<https://urldefense.proofpoint.com/v2/url?u=https-3A__connect.microsoft.com_VisualStudio_feedback_details_995447_linker-2Dis-2Dhanging-2Dwhile-2Dlinking-2Dsome-2Dproject-2Din-2Ddebug-2Dx64-2Dconfiguration=CwMFaQ=ODFT-G5SujMiGrKuoJJjVg=bn1clvQAiDQrfAC4yKbN0PlSr7RWRs-U3WJ0zRBB2qM=c_xwp0Rh13-GKNF7fVcUZoqlcZj2VZXblnm-2ZUnrhw=reRA7jaZJcHsEqlwnTZVrPH_wspoTOncSAPUu52tDiM=>

https://www.gittprogram.com/question/1740272_x64-linker-hanging.html

<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.gittprogram.com_question_1740272-5Fx64-2Dlinker-2Dhanging.html=CwMFaQ=ODFT-G5SujMiGrKuoJJjVg=bn1clvQAiDQrfAC4yKbN0PlSr7RWRs-U3WJ0zRBB2qM=c_xwp0Rh13-GKNF7fVcUZoqlcZj2VZXblnm-2ZUnrhw=qm7zUzMXCYYVlnMjutGub8L_miInfoGNbU39xKQQIXc=>


http://stackoverflow.com/questions/23494142/ms-c-2012-linker-hangs-but-only-in-debug-mode-x64

<https://urldefense.proofpoint.com/v2/url?u=http-3A__stackoverflow.com_questions_23494142_ms-2Dc-2D2012-2Dlinker-2Dhangs-2Dbut-2Donly-2Din-2Ddebug-2Dmode-2Dx64=CwMFaQ=ODFT-G5SujMiGrKuoJJjVg=bn1clvQAiDQrfAC4yKbN0PlSr7RWRs-U3WJ0zRBB2qM=c_xwp0Rh13-GKNF7fVcUZoqlcZj2VZXblnm-2ZUnrhw=ht7u_cKJxrXj5gSPwhkpJq9p8ZVyUHqVV7k71iMZULM=>


Hard to tell whether these are related.



I think i'll wait for the patch allowing to build deal.II as a dynamic
library.  I usually develop in Ubuntu but i also need my program to run on
Windows. For now, compiling only in Release mode is good enough.


OK, good to know that you have an alternative way to work.

Best
 W.


--
----
Wolfgang Bangerth   email:bange...@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: Development activity over the past year

[Wolfgang Bangerth]

I'm writing my annual report to the NSF and wanted to share a few statistics I
thought were interesting:

* 40,000 new lines of code in deal.II (excluding testsuite) over the past year
* 77 new testcases
* 4 new tutorial programs
* 1540 emails to dealii@googlegroups.com between January to July alone
* 725 people are subscribed to the mailing list (+190 compared to last year)


I could also add this:
* On average, more than 1,000 downloads per month (excluding those who work
  directly from the github repository)
* We merged 3.4 pull request per day on average, weekday or weekend.
* 54 distinct people contributed code and documentation over the last
  year.

I am particularly proud of the last number!

Cheers
 W.

--

Wolfgang Bangerth   email:bange...@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] Development activity over the past year

[Wolfgang Bangerth]

All,
I'm writing my annual report to the NSF and wanted to share a few statistics I 
thought were interesting:


* 40,000 new lines of code in deal.II (excluding testsuite) over the past year
* 77 new testcases
* 4 new tutorial programs
* 1540 emails to dealii@googlegroups.com between January to July alone
* 725 people are subscribed to the mailing list (+190 compared to last year)

What this shows to me is that this is a vibrant community that is doing 
awesome things and helping each other! Thanks to you all for participating in 
the many ways in which you contribute!


Best
 Wolfgang

PS: If you haven't actively participated so far, we would love to have you be 
a more active part of this project! Here are a few ideas:

  https://github.com/dealii/dealii/blob/master/CONTRIBUTING.md
  https://www.dealii.org/participate.html

--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] HP with parallel::shared::Triangulation

[Wolfgang Bangerth]

On 07/31/2016 12:06 PM, Pete Griffin wrote:


I get an error from an assertion in the call to:
   dof_handler.distribute_dofs (fe_collection);
The triangulation is a parallel::shared::Triangulation.

My C++ is not as good as it needs to be. It appears as if the exception
(below) means that HP is not implemented for a
parallel::shared::Triangulation. Is my interpretation correct ?

Is parallel::distributed::Triangulation implemented ?


Pete -- correct, hp finite elements are not implemented for either of the two 
parallel triangulations. I.e., they can currently only be used for sequential 
computations. This is pretty high on our list to get implemented in parallel.


Best
 W.


--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Re: PETScWrappers::SparseDirectMUMPS for symmetry system matrix and factorization resue

[Wolfgang Bangerth]

On 07/27/2016 01:43 PM, Alexander wrote:


step-17 claims to solve same problem as step-8. If one produces symmetry
matrix and the other does not or both do not -- this is suspicious.


The reason why the matrix in step-17 is not symmetric is actually 
explained in the code:


// The last argument to the call to
// MatrixTools::apply_boundary_values() below allows for some
// optimizations. It controls whether we should also delete
// entries (i.e., set them to zero) in the matrix columns
// corresponding to boundary nodes, or to keep them (and passing
// true means: yes, do eliminate the columns). If we
// do eliminate columns, then the resulting matrix will be
// symmetric again if it was before; if we don't, then it
// won't. The solution of the resulting system should be the same,
// though. The only reason why we may want to make the system
// symmetric again is that we would like to use the CG method,
// which only works with symmetric matrices. The reason why we may
// not want to make the matrix symmetric is because this
// would require us to write into column entries that actually
// reside on other processes, i.e., it involves communicating
// data. This is always expensive.
//
// Experience tells us that CG also works (and works almost as
// well) if we don't remove the columns associated with boundary
// nodes, which can be explained by the special structure of this
// particular non-symmetry. To avoid the expense of communication,
// we therefore do not eliminate the entries in the affected
// columns.
std::map boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
  0,
  ZeroFunction(dim),
  boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
solution,
system_rhs,
false);

In other words, it is expected that the matrix is not symmetric, and 
that consequently when you run a direct solver, you can't use the 
symmetric setting.


Best
 Wolfgang

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Observations on CPU usage on a single processor

[Wolfgang Bangerth]

Pete,


I was surprised to see all 4 cores (8 threads) being used in step-8 which does
not have PETSc/MPI support. In the solver, the main thread was using 100% of
it's resources. The others were using very consistently 24-25% independent of
#DOF. Is the normal SolverCG<> of step-8 multi-theaded?


Yes. Various parts of the library are multi-threaded and doing matrix-vector 
multiplications is one of those parts.




Would the overall
performance on a single processor improve if one were able to have multiple
controlling threads instead of one. Is it even possible? This might be even
more important for the 8 or 16 core CPUs I see available now. If I am correct
the code that exists was probably optimum with Single or Dual Core CPUs but it
does not appear to be true now. It may not be feasible but I though I'd
mention it.


It is possible to write some programs in a way so that there are multiple 
controlling threads, but this is not the case for typical finite element 
programs. If you look at a typical program, it almost always has the form

do this   (e.g., assemble the linear system)
  then
do that   (e.g., solve the linear system)
The "then" is a synchronization point. However, "do this" and "do that" may 
well be parallelizable, and we often do that in the library if possible. An 
example may be the loop over all cells in assembly. But it is difficult to do 
this with all statements, and that is where the inefficiency comes from.




Also, it appears that the efficiency with one CPU decreases as the #DOF
increase with PETSc/MPI. Does anyone know why?


How do you define "efficiency"?

Best
 W.

--
----
Wolfgang Bangerth   email:bange...@colostate.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] Daniel Arndt and Jean-Paul Pelteret appointed deal.II Developers

[Wolfgang Bangerth]

All,
I want to share the good news that Daniel Arndt and Jean-Paul Pelteret have 
accepted our invitation to become deal.II Developers. They are now listed at

  https://dealii.org/authors.html
Both of them have provided patches throughout the library for several years 
already, but you all probably know them best from their tireless efforts in 
keeping the mailing list well organized and making sure that most of the 
questions get answers.


Daniel & Jean-Paul, many thanks already for your past contributions, and we 
look forward to many more in the future!


Happy hacking
  Wolfgang Bangerth

--
----
Wolfgang Bangerth   email:bange...@colostate.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.

Re: [deal.II] Re: Entry of the own problem in dealii

[Wolfgang Bangerth]

On 07/12/2016 07:14 AM, Denis Davydov wrote:


thanks for your answer. Ok I will ask the question in a new post,but to your
information I compiled Trilinos-12.4.2 and I become my mpi version 5.1.2 as I
give the command 'mpiexec --version' bud during dealii compiling I see that
the MPI_VERSION=2.0.


This indicates that your MPI implementation (e.g. mpich, or openmpi, or Intel 
MPI) is 5.1.2, but that it supports the MPI Standard 2.0. There is no conflict 
I believe.


This is similar to saying that GCC 4.9.4 supports the C++11 standard.

Best
 W.

--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Re: Problem building dealii with PETSc

[Wolfgang Bangerth]

On 07/12/2016 09:59 AM, Pete Griffin wrote:


sudo apt-get install libopenmpi-dev
sudo apt-get install libmpich-dev(12.1.0)


This may not get all of the MPI libraries, but only the development 
parts (e.g., the header files, etc). Is there no package simply called 
"openmpi" or "mpich"?



> Is my inteperetation of the instructions correct, or do I have to
> remove openmpi-dev.

I'd try to start from as clean a slate as possible every time.

Best
 W.

--
----
Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Re: Problem building dealii with PETSc

[Wolfgang Bangerth]

deal.II does not use the C++ MPI interface (I do not know of any project that
does, for the record),


Which is, I guess, why the MPI committee has deprecated them :-)

But back to the issue at hand:



I used downloaded PETSc source and built with:

./configure --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-mpich


That last flag is generally a bad idea. It makes PETSc download and install an 
MPI version that may or may not conflict with any MPI package that is already 
installed on the system. In any case, it is going to be installed into some 
non-standard path where the remainder of the system will scratch its head 
trying to find out whether there is or is not an MPI installation. I wished 
the PETSc people just got rid of this option -- it has caused so many people 
so much heartache.


The way to install MPI is through the system's package manage system. Any 
distribution I can think of has either an openMPI or an mpich package. Either 
is fine (but not trying to install both). Once you've installed one of them on 
your system (via apt-get, zypper, or whatever package manager your system 
uses), you can install PETSc as you were doing, but without --download-mpich.


Best
 W.

--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Dealing with conflicting constraints

[Wolfgang Bangerth]

On 07/04/2016 08:52 AM, Ce Qin wrote:

Dear all,

I am now solving a problem using adaptive finite element method. The boundary
condition is Dirichlet condition. When the cells on boundary are locally
refined, the hanging node constraints are conflicting with inhomogeneous
constraints.

  In Constraints on degrees of freedom
<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dealii.org_8.4.1_doxygen_deal.II_group-5F-5Fconstraints.html=CwMFaQ=ODFT-G5SujMiGrKuoJJjVg=bn1clvQAiDQrfAC4yKbN0PlSr7RWRs-U3WJ0zRBB2qM=1SpqN_97meGfz2Xl2-E8Tbg0or5hD9WqUYXEKaSyPYc=g6gaCUE0kHiyFVTtZYEm3loEhpwmQKH7FkZCX1J1OIE=>,
 section
Dealing with conflicting constraints, it says:

If you want the hanging node constraints to win, then first build these
through the DoFTools::make_hanging_node_constraints()

<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dealii.org_8.4.1_doxygen_deal.II_group-5F-5Fconstraints.html-23ga3eaa31a679484e80c193e74e8a967dc8=CwMFaQ=ODFT-G5SujMiGrKuoJJjVg=bn1clvQAiDQrfAC4yKbN0PlSr7RWRs-U3WJ0zRBB2qM=1SpqN_97meGfz2Xl2-E8Tbg0or5hD9WqUYXEKaSyPYc=6UnsFr5iWSlrA4aiDNvfPd3uy0hbIygrP-fKWn6-mXQ=>
 function.
Then interpolate the boundary values using
VectorTools::interpolate_boundary_values()

<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dealii.org_8.4.1_doxygen_deal.II_namespaceVectorTools.html-23af6f700f193e9d5b52e9efe55e9b872d5=CwMFaQ=ODFT-G5SujMiGrKuoJJjVg=bn1clvQAiDQrfAC4yKbN0PlSr7RWRs-U3WJ0zRBB2qM=1SpqN_97meGfz2Xl2-E8Tbg0or5hD9WqUYXEKaSyPYc=5xj02EL9VScyZ6Q87iec6SjedHT33BYhVXmFNT8hVGY=>
 into
the same ConstraintMatrix

<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dealii.org_8.4.1_doxygen_deal.II_classConstraintMatrix.html=CwMFaQ=ODFT-G5SujMiGrKuoJJjVg=bn1clvQAiDQrfAC4yKbN0PlSr7RWRs-U3WJ0zRBB2qM=1SpqN_97meGfz2Xl2-E8Tbg0or5hD9WqUYXEKaSyPYc=QtUsjdtvZhO54EaZyqL2V2RDq87D7HIdGqo1DKLz2Zc=>
 object.
If the latter function encounters a boundary node that already is
constrained, it will simply ignore the boundary values at this node and
leave the constraint untouched.


But in documentation of VectorTools::project_boundary_values_curl_conforming
<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dealii.org_8.4.1_doxygen_deal.II_group-5F-5Fconstraints.html-23ga1c6685360c01c9c46eeb7575e8ef68ac=CwMFaQ=ODFT-G5SujMiGrKuoJJjVg=bn1clvQAiDQrfAC4yKbN0PlSr7RWRs-U3WJ0zRBB2qM=1SpqN_97meGfz2Xl2-E8Tbg0or5hD9WqUYXEKaSyPYc=bQl5C_-vKhJvAZZDD0NqBMW3JuklcaGy-9YSrXjdwuk=>,
it says:

If the ConstraintMatrix

<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dealii.org_8.4.1_doxygen_deal.II_classConstraintMatrix.html=CwMFaQ=ODFT-G5SujMiGrKuoJJjVg=bn1clvQAiDQrfAC4yKbN0PlSr7RWRs-U3WJ0zRBB2qM=1SpqN_97meGfz2Xl2-E8Tbg0or5hD9WqUYXEKaSyPYc=QtUsjdtvZhO54EaZyqL2V2RDq87D7HIdGqo1DKLz2Zc=>|constraints|
 contained
values or other constraints before, the new ones are added or the old ones
overwritten, if a node of the boundary part to be used was already in the
list of constraints. This is handled by using inhomogeneous constraints.
Please note that when combining adaptive meshes and this kind of
constraints, the Dirichlet conditions should be set first, and then
completed by hanging node constraints, in order to make sure that the
discretization remains consistent. See the discussion on conflicting
constraints in the module on Constraints on degrees of freedom

<https://urldefense.proofpoint.com/v2/url?u=https-3A__www.dealii.org_8.4.1_doxygen_deal.II_group-5F-5Fconstraints.html=CwMFaQ=ODFT-G5SujMiGrKuoJJjVg=bn1clvQAiDQrfAC4yKbN0PlSr7RWRs-U3WJ0zRBB2qM=1SpqN_97meGfz2Xl2-E8Tbg0or5hD9WqUYXEKaSyPYc=g6gaCUE0kHiyFVTtZYEm3loEhpwmQKH7FkZCX1J1OIE=>.

I was totally confused. Which one is right?


It seems to me that both of these texts describe what the functions do. But 
they do different things. I think the first one does the correct thing whereas 
the second may not. Want to fix this in the code and the documentation?


Best
 W.

--
----
Wolfgang Bangerth   email:bange...@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.

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

[Wolfgang Bangerth]

On 07/05/2016 10:08 AM, Ehsan Esfahani wrote:


thanks for your response. the problem is in preconditioner. I find out
that my code is running on my friend's machine but I get the error
related to SEGV. I have changed preconditioner to Jacobi and it's
running but with more iterations than AMG. Also, it's not running on my
computer and the problem on my computer remained the same as before :-(


Ehsan,
then let us come back to my original suggestion: when you run a program 
in a debugger and it crashes, the backtrace will show where this 
happens. If the problem was really in the line

  solver.solve (system_matrix, completely_distributed_solution_update,
system_rhs, preconditioner);
then that would mean that in creating one of the arguments to this call, 
the program did something illegal. For example, it could be that one of 
the arguments is a reference that is not bound to anything. You should 
be able to see this by looking at the various arguments in the debugger.


More likely, however, is that the program actually crashed *inside* the 
`solver.solve` function, or one of the functions that it calls. You can 
see this too in the backtrace, by seeing which functions are being 
called. Once you find out which function creates the problem, you can 
again look at the arguments and local variables and see what is happening.


Best
 W.

--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] (Multigrid) preconditioner for the PML method

[Wolfgang Bangerth]

On 07/05/2016 10:46 AM, Artur Safin wrote:


I am working on solving a Helmholtz-type problem, where I am using the
PML method to simulate absorbing boundary conditions. Originally (before
adding PML) I was using the complex GAMG preconditioner from the PETSc
package, which was working fine. But now I am having a difficulty where
the solver does not converge due to PML.

If anyone has any suggestions or experience, that would be greatly
appreciated.


Artur -- there are two options: either the matrix is not invertible, or 
the preconditioner is not a good preconditioner for the matrix. My first 
step in debugging such problems is usually to try a direct solver such 
as SparseDirectUMFPACK first, in order to make sure the matrix is really 
invertible and yields the correct solution. Once I'm sure of that, one 
can start thinking of the preconditioner.


Best
 W.


--

Wolfgang Bangerth   email:bange...@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.

Re: [deal.II] Re: Observations on CPU usage on a single processor

[Wolfgang Bangerth]

Pete,


Bruno, I assumed that the thread with 100% CPU usage was somehow feeding the
others in step-8,


It's more like for some functions, we split operations onto as many threads as 
there are CPUs. But then, the next function you call may not be parallelized, 
and so everything only works on one thread. On average, that one thread has a 
load of 100% whereas the others have a lesser load.




I just tested the step-8 program with PreconditionIdentity() witch showed 100%
CPU usage on all 8 CPUs. The results follow. Assuming having no preconditioner
only slows things down maybe getting 3 times the CPU power will make it up. I
haven't checked solve times yet. The preconditioner for step-8 was
PreconditionSSOR<> with relaxation parameter = 1.2. Is there an optimum
preconditioner/relaxation parameter for 3d elasticity problems that you know
of? Is their determination only by the trial and error?


1.2 seems to be what a lot of people use.

As for thread use: if you use PreconditionIdentity, *all* major operations 
that CG calls are parallelized. On the other hand, using PreconditionSSOR, you 
will be spend at least 50% of your time in the preconditioner, but SSOR is a 
sequential method where you need to compute the update for one vector element 
before you can more to the next. So it cannot be parallelized, and 
consequently your average thread load will be less than 100%.


Neither of these are good preconditioners in the big scheme of things, if you 
envision going to large problems. For those, you ought to use variations of 
the multigrid method.




Wolfgang, What I meant by efficiency was the CPU usage in the threads for
Step-17 NEW and OLD decreased with the larger #DOFs or cycle #'s.


If the load decreased for both codes, I would attribute this to memory 
traffic. If the problem is small enough, much of it will fit into the caches 
of the processor/cores, and so you get high throughput. If the problem becomes 
bigger, processors wait for data for longer. Waiting is, IIRC, still counted 
as processor load, but it may make some operations that are not parallelized 
take longer than those that are parallelized, and so overall lead to a lower 
average thread load.


But that's only a theory that would require a lot more digging to verify.

Best
 W.


--
----
Wolfgang Bangerth   email:bange...@colostate.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.

Re: [deal.II] Re: Observations on CPU usage on a single processor

[Wolfgang Bangerth]

Pete,


When I started to try to implement PETSc/MPI I didn't realize the normal
solverSG<> was running with all CPUs. As of now I have no need for PETSc/MPI.
I did a little more looking and found that on a single processor it was no
faster than the non-PETSc/MPI. For the 3D elastic problem I was working on I
found, also, by trial and error, that SSOR was fastest even with 25% CPU usage
with a relaxation parameter = 1.2. Jacobi which showed 80-90% CPU usage in the
threads was actually slower.


Yes, that's what you often find: the algorithms that are easy to parallelize 
are just not good enough to really compete with the more complex algorithms 
(which are then harder to parallelize efficiently).


Best
 W.

--
----
Wolfgang Bangerth   email:bange...@colostate.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.

Re: [deal.II] problems switching from one core to multiple cores - implementation of Chorin's scheme

[Wolfgang Bangerth]

On 08/02/2016 12:45 PM, Marek Čapek wrote:


I coded it in the step-40 like manner.
When I run it on one core (either ./main or mpirun -np 1 main ), the
results are physically reasonable.
However, when I switch to more cores, e.g. two - "mpirun -np 2 main", I
get meaningless solutions.
You know, that the Chorin method is three-step, wherein it is computed
with velocity in the first and third step.
I use the same degree of freedom handler for the solution vectors of the
first and third step. I initialize the vectors with
the same sets of dofs - locally relevant dofs  & locally owned dofs.
Could the problem originate from there?
I compiled the problem on my computer as well as on cluster and the
results are the same.
I suppose, that the scheme is coded correctly as it works on single core...


Marek -- "meaningless" leaves many options open. Since it's a three-step 
scheme, have you tried to output the solution after the first step and 
compare what you get with one and two processors? If you look at the 
solution, does it give an indication of what may be going wrong in 
parallel? (E.g., you could potentially see if something is wrong with 
boundary values, hanging node constraints, etc.)


Best
 W.

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Development activity over the past year

[Wolfgang Bangerth]

On 08/02/2016 02:33 PM, David Wells wrote:


About a year ago I did a bit of analysis on the commit rate: it has gone
up significantly over the least few years. Just for fun I attached a
plot of the total number of commits over time (as well as some important
points in our history) to visualize this. As far as I can tell, Matthias
joining the project increased the slope of the line remarkably :)


Matthias has done an incredible amount of work rewriting the 
configuration system, often late into the night based on the times when 
many of his patches were committed.


But him joining also happened around a time when the number of people 
contributing to the project grew significantly. Timo had made a graph a 
while back showing the number of individuals contributing patches in any 
given month, and it also showed a marked increase starting around that 
time. The whole project has been growing in recent years, and I'm 
excited about that!


Best
 W.

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] problems switching from one core to multiple cores - implementation of Chorin's scheme

[Wolfgang Bangerth]

On 08/02/2016 02:18 PM, Vinetou Incucuna wrote:


I have tried to output the solutions from the step one as
well as from the step three. There is some development, but
reasonable,i.e. solution from the first step and the third step look
similar (except the first iteration)
In the attachment you can see the results of my computations
and the partitions.
It appears to me, that the problematic edges in the solution (both
pressure and velocity) correspond the borders between the partitions.
Maybe you could give me a gist, what is going on.


It's hard to tell. The code doesn't look wrong, but without the context 
I also can't say that it's right (and I wouldn't have the time to look 
at it in sufficient detail even if I had the context). My best guess is 
that you may be assembling on the wrong cells. E.g., you may be 
assembling also on ghost cells, instead of only on locally owned ones. 
Of course I don't know whether that's true.


Are you running on a mesh without hanging nodes? What happens if you run 
just on a regular square instead of a complicated domain? My strategy 
for debugging is always to make the problem simpler and simpler.


Best
 W.

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
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

[Wolfgang Bangerth]

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

Re: [deal.II] error when trying to initialize PETSc preconditioners

[Wolfgang Bangerth]

On 08/15/2016 02:46 PM, Ignacio Tomas wrote:

Great! Many thanks for solving this problem. I owe at least a lunch.


Ah, now I'm jealous I left this easy to fix bug to Timo! ;-)

Cheers
 W.

--

Wolfgang Bangerth  email: bange...@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: [dealii-developers] dealii.org website down?

[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: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Extracting a component from Block FESystem and DofHandler

[Wolfgang Bangerth]

On 08/16/2016 09:34 AM, Spencer Patty wrote:


Supposing velocity_block represents the block number associated to the
velocity components, I noticed that solution.block{velocity_block} would be a
TrilinosWrapper::MPI::Vector of the velocity as desired and it is possible to
extract the sub FiniteElement fe_system.base_element{velocity_block}
corresponding only to the velocity, but the dof_handler has no such existing
capability and if I pass it as is with the extracted FiniteElement into an
FEValues, an assertion is triggered when we call fe_values.reinit(cell) saying
that the given finite element does not match the dofhandler's finite
element...


Correct -- there is no way to subset a DoFHandler by vector components. That's 
because a FESystem is really stored as a collection of simpler elements, 
whereas a DoFHandler is not.




My question is this:  What would be a best approach to dealing with this?  Is
there currently a way to extract a "diagonal block" from the system and reduce
the solution, fesystem and dofhandler to only represent that part?  Do I need
to create a new FESystem, DoFHandler and vector to provide the functionality I
need ?


The first two yes. If you do it right, i.e., ensure that the velocity 
DoFHandler numbers degrees of freedom in the same order as the system 
DoFHandler does, then you won't need a separate solution vector -- you can 
just use a block of your block vector.


The alternative -- and I think that is what Bruno was suggesting in his first 
email -- is to keep the one system DoFHandler and vector. In addition to 
those, you would then have to pass something to all consumers of the velocity 
field that identifies which components of the FESystem and of the block vector 
they need to know about.


Cheers
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Re_Error during installation of dealii.8.4.1

[Wolfgang Bangerth]

On 02/02/2017 03:41 PM, Sathish Kumar wrote:

//
/icpc: command line warning #10006: ignoring unknown option '-fuse-ld=gold'/
/lac/CMakeFiles/obj_lac.debug.dir/sparse_matrix_inst2.cc.o: file not
recognized: File format not recognized/


This typically happens if you compile in the same directory in which a 
previous compile job aborted because of a compiler error, because you 
were out of memory, or some other reason. In that case, the previous 
compiler run left an unfinished object file that the subsequent compiler 
run recognizes as "up to date" (because it's newer than its source 
files) but can't read successfully.


The solution is to do 'make clean' or just to blow away the build 
directory and start again.


Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] internal compiler error: Bus error

[Wolfgang Bangerth]

On 02/01/2017 06:23 PM, Lev Karatun wrote:


I'm trying to install the latest version of dealII on a cluster now, got
another error unfortunately.
The output is below, and I attached the logs too. I would appreciate any
help with this.


Are you compiling with -j4 or similar? Does the same happen if you 
compile with -j1? The kind of internal compiler error you have here 
seems to happen most often if you run out of memory.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] cannot call member function ‘void tbb::task::spawn(tbb::task&)’ without object

[Wolfgang Bangerth]

On 01/30/2017 08:47 PM, Lev Karatun wrote:


I was trying to install (to an empty directory) the newest version of dealII
but am getting the following error:
/include/deal.II/base/thread_management.h:2921:32: error: cannot call member
function ‘void tbb::task::spawn(tbb::task&)’ without object


Lev -- are you using the current development sources?

I think this was introduced by #3717:
  https://github.com/dealii/dealii/pull/3717
In the version of the threading building blocks Karl is using (and apparently 
all of us are using), the tbb::task::spawn() function is a static member function.


What version of the TBB are you using? And can you look up the declaration of 
the spawn() function in task.h in your version of the TBB?


Best
 W.


--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Fully distributed triangulation (level 0)

[Wolfgang Bangerth]

YC,


I have a project requiring to read in a large coarse mesh from gmsh to deal.ii

1M dofs. Most of cells have their own characteristics, which means I cannot

combine them and create a coarse mesh.
Currently, I implemented it by using shared memory triangulation for
parallelization. Since I want to scale it to a cluster system and target a
100M mesh (no need for mesh refinement), I have to use distributed tria via
MPI (is there any better solution?). I found out that the initial cost is
large because of the duplication of triangulation and p4est forest. I was
wondering if there is any method to remove part of triangulation or p4est
data.


No, unfortunately there is not. p4est is built on the assumption that the 
coarse mesh is replicated on all processors, and deal.II inherits this 
assumption. If your coarse mesh has 1M cells, that may just barely so be 
tolerable, although it will likely lead to inefficient code in some places 
where you loop over all cells and almost all of them turn out to be 
artificial. But I suspect that you will be in serious trouble if your coarse 
mesh has 100M cells.


You should really try to come up with a coarser coarse mesh that you can then 
refine.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] using VectorTools::interpolate_boundary_values :Simple verification of Stokes system with zero velocity (using step-22)

[Wolfgang Bangerth]

I'm trying to solve a stokes system with some funny boundary conditions
so I've stripped the code to make sure I was getting the correct answers
to a very simple system as follows:
governing equations are:
div grad u + grad p = f
div u = 0
i'm setting u=0 on the boundaries, giving u=0 (or close enough) everywhere.
setting f = -(0,1), i should get a linear profile with pressure
increasing with depth at the same rate going downwards.


Yes. But if you impose Dirichlet boundary values for the velocity all 
around the boundary, then the pressure is only determined up to a 
constant, i.e., the pressure should be a linear function but the offset 
is not determined by the equation.




I'm clearly not understanding how to use
VectorTools::interpolate_boundary_values correctly, and  I would
appreciate some clarification on it.

I'm imposing a zero pressure at the top of my rectangle with height of
10, which means I should get 10 at the bottom for pressure.


You can't do that. If you have Dirichlet boundary conditions for the 
velocity, then you can't impose anything on the pressure. In fact, for 
the incompressible Stokes equations, you can never impose any pressure 
boundary conditions.


There is a description of boundary conditions for the Stokes equations 
in the introduction of step-22. Take a look there to see what you can 
and cannot impose for the Stokes equations.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Possible Mistake and improvement in QMRS implementation

[Wolfgang Bangerth]

On 01/31/2017 11:24 AM, Ingo Kligge wrote:


during my work with the Krylov subspace method QMRS in deal.ii I met
some problems: the implementation of the method is, as stated in the
class description, adapted from Algorithm 5.1 in /Freund/Nachtigal:
Software for simplified Lanczos and QMR algorithms, Appl. Num. Math. 19
(1995), pp. 319-341/ (e.g http://dl.acm.org/citation.cfm?id=223396)
specialized for the solution of right-preconditioned symmetric linear
systems. The deal.ii version however calculates (at least up to version
8.4.1) the *preconditioned* initial residual - I think, referring to the
cited paper and any other right-preconditioned Krylov subspace method,
it should be the *unpreconditioned* one. I've noticed very bad
convergence behaviour and breakdowns as well. What do you say? Do you
agree with my conjecture?


I implemented this ~17 years and have no recollection. I would be 
surprised if anyone else has looked at it in the meantime. If you say 
that there is a bug, I have no reason to doubt that :-)



Beside of that: I would suggest a different implementation of the
symmetric QMR-method. In an earlier paper
(https://www.researchgate.net/publication/234171461_A_new_Krylov-subspace_method_for_symmetric_indefinite_linear_system
, 1994) Freund and Nachtigal propose the essentially same algorithm
especially for solving symmetric indefinite linear systems with a
symmetric precondition matrix, that is a) flexible for right, left and
split preconditioning and b) compared to the existing deal.ii code
requires only ONE application of the precondition matrix per iteration
in the case of right preconditioning (instead of two). The best, I
think, is a hybrid of both, avoiding squaring the scalars for updating
the iterates because of possible rounding errors.

Maybe I will rewrite the method on my own. if I find the time.


Since you seem to know far more about this than I suspect any of the 
other developers, I think it would be great if you could "fix" the 
implementation we have, or write a better one. Any help would definitely 
be much appreciated!


Best
 Wolfgang

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Shape function derivatives

[Wolfgang Bangerth]

On 02/01/2017 07:02 AM, 'Seyed Ali Mohseni' via deal.II User Group wrote:

My approach is the following:

|
template
FullMatrix get_b_operator (const FEValues _values, const
unsigned int dofs_per_cell, const unsigned int q)
{
FullMatrix tmp(dim, GeometryInfo::vertices_per_cell);

std::vector<DerivativeForm<1, dim, dim> > invJ =
fe_values.get_inverse_jacobians();

for (unsigned int i = 0; i < dofs_per_cell; i = i + dim)
{
const unsigned int index = i / dim;

// COMPUTE: B-operator (Remark: This version has to be extended for 3D!)
tmp[0][index] = invJ[q][0][0] * fe_values.shape_grad(i, q)[0] +
invJ[q][0][1] * fe_values.shape_grad(i, q)[1];
tmp[1][index] = invJ[q][1][0] * fe_values.shape_grad(i, q)[0] +
invJ[q][1][1] * fe_values.shape_grad(i, q)[1];
}

return tmp;
}
|


fe_values.shape_grad() returns the gradient *on the actual cell*, not on 
the reference cell. Consequently, you do not have to multiply it by invJ 
again. This is the source of your factor of two.


Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Wolfgang Bangerth]

On 02/08/2017 01:15 PM, Jaekwang Kim wrote:


I am using deal-ii to analyze Non-newtonian fluids, which varies viscosity.
For example, my viscosity eta is function of (u) field. I am solving
this with iteration method.. (Picard)
After coding, I am testing my code with method of manufacture.

For squared domain, I implemented all boundary condition with my assumed
solution and impose source term also to numerically generate my assumed
solution
however, the problem I got is also follow


While my velocity field is qualitatively similar to manufactured solution,
the pressure field is not following the exact form as seen below figure.
It is very strange that..I can approach u-soltution form while my
pressure does not


In addition , any mesh refinement would not make more accurate solution
now.

If you were me, what would you check more to find out what is problem?


I would make the problem simpler. Do you get the right solution if you 
have a constant viscosity? If you iterate, do you get the solution after 
one iteration?


I will also note that your manufactured solution is symmetric, but your 
computational solution is not. This provides you with a powerful way 
because you don't actually need to look at the convergence, it's enough 
to check that on a *coarse* mesh the solution is *symmetric*. It may be 
that solution is already non-symmetric on 1 cell, or on a 2x2 mesh -- in 
which case you already know that something is wrong. The question then 
is why that is so -- is your rhs correct, for example?


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Wolfgang Bangerth]

On 02/08/2017 12:35 PM, Jaekwang Kim wrote:

   */_local_rhs(i) += fe_values.shape_value(i,q) *_/*

*/_rhs_values[q] *_/*

*/_fe_values.JxW(q);_/*



rhs_values[q] is a Vector. You need to say which component of it 
you want to multiply with. (Hint: The correct component is 
fe.system_to_component_index(i).first).


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Re_Error during installation of dealii.8.4.1

[Wolfgang Bangerth]

On 02/04/2017 07:18 AM, Timo Heister wrote:

Now the installation says that there is no write permission on that
particular drive after 100% completion.


what is your CMAKE_INSTALL_PREFIX set to? It should be set to
something in your home directory.


In particular you set it to /path/to/install/dir/, which is a path you want to 
replace by the one where you want to install deal.II in :-)


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] using VectorTools::interpolate_boundary_values :Simple verification of Stokes system with zero velocity (using step-22)

[Wolfgang Bangerth]

On 02/01/2017 11:24 AM, Jane Lee wrote:


where the boundary values are just 1.0 so that u is a constant vector field.
this still gives me the bizarre 22 factor in my pressures solution
(which is just linear in z - the resulting pressure solution is atleast
linear), but the pressure solution aside, both the x and y components of
the velocity are 0.0 ON the boundaries and 1.0 everywhere else.


Are you calling constraints.distribute on the solution vector after 
solving your linear system?


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] cannot call member function ‘void tbb::task::spawn(tbb::task&)’ without object

[Wolfgang Bangerth]

On 01/31/2017 02:53 PM, Bruno Turcksin wrote:

According to your detailed.log, you are not.  You are using TBB 2.2
which is pretty old. Can you reconfigure and recompile deal.II using -D
|DEAL_II_FORCE_BUNDLED_THREADS=ON|


Or wait for the patch I'm going to submit in a minute :-)
W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] cannot call member function ‘void tbb::task::spawn(tbb::task&)’ without object

[Wolfgang Bangerth]

Or wait for the patch I'm going to submit in a minute :-)


Specifically, I would appreciate if you could see whether this helps:
  https://github.com/dealii/dealii/pull/3889

Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Fully distributed triangulation (level 0)

[Wolfgang Bangerth]

Yi-Chung,


Thank you for your prompt reply. I was wondeirng if I can integrate
other partition tools, such as metis or parmetis to handle the fully
distributed triangulation. I can develop that part by myself (or with
some help from the community). Do you have any suggestions?


I believe that it would not be impossible to develop this, but it will 
probably not be a small exercise. You would have to develop a fair share 
of code, and gain an understanding of parts of the library you don't 
usually get to see if you use deal.II.


If you're willing to do this, we can guide you along the process, but 
it's going to be a bit of work for sure.




My following
project also relies on this since I will try to manage number of cells
in each processor. With p4est, it is hard to manage number of cells
based on a in-house algorithm.


That is actually not quite true. p4est (and the deal.II classes that use 
it) allow you to describe a "weight" for each cell, and p4est will 
partition in a way that the sum of weights is roughly equal among all 
processors.




My application is about IC designs that
may have million to billion cells. A fully distributed triangulation
helps to reduce memory usage. The current shared_memory system can
handle 20M (single core) in system of 32GB main memory.


That's already quite impressive :-) What kind of meshes do you have that 
require so many cells? Are they geometrically incredibly complicated to 
require that many cells already at the coarse level?




Any design of 1M cells on the distributed triangulation have problem of
computation time because of the reorder step. This is why I bypassed it
and provided a sorted mesh to grid_in (read_msh). For a problem of 5M
cells, I can save 200sec at the step of create_triangulation.


Yes, I'm willing to believe this. The algorithm wasn't intended for 
meshes of this size, though we did test it with ~300k cells in 2d and we 
know that it scales like O(N). So 200 seconds seems like a long time. Is 
this in debug mode?


Best
 Wolfgang


--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Add heterogeneous parameters controlled by function on cells

[Wolfgang Bangerth]

On 01/26/2017 01:20 AM, Lam DANG wrote:


But I still a concern that the function class (the Coefficient class in step
6) will be called in many time in the main code.


Yes, but so are many other functions (e.g., FEValues::shape_grad). Is this 
function *particularly* expensive in your case?


I'm asking because unless you have *concrete* evidence that something is slow, 
then it's not worth changing -- in particular, if the result is code that is 
more complex. There is a saying in computer science that "premature 
optimization is the root of all evil", and I think there is truth to that. So, 
*unless you have measured that it is a problem*, then the simplest approach is 
the way to go.


Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Assembly of material forces

[Wolfgang Bangerth]

On 01/28/2017 04:56 AM, 'Seyed Ali Mohseni' via deal.II User Group wrote:


I somehow figured the assembly out by listening to Prof. Bangerth's advice
that I should try to solve it the deal.II way. And after Jean-Paul told me to
take a look at step-42 I realized how things are working in deal.II.
Unfortunately, I have still some issues with the values. They are not the
same. I believe this is dependent on the B-operator, namely the jacobians.
They differ by factor 2 probably since our unit cell in deal.II is within
[0,1]^dim and the C++ code uses [-1,1]^dim.
But shouldn't the assembled configurational forces still be the same
independent from the mapping?

Here is what I got so far :)


I really have no idea what the pictures show, nor which ones correspond to the 
deal.II solution and your own code, respectively.


But I do see that one seems to be computed on a finer mesh -- maybe a 
triangular mesh, whereas the other one is computed on a quadrilateral mesh.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Add heterogeneous parameters controlled by function on cells

[Wolfgang Bangerth]

On 01/25/2017 05:37 AM, Jean-Paul Pelteret wrote:

Dear Lam,

There are a number of ways that you can do this. One approach would be to
define a class that derives from the function class
<https://www.dealii.org/developer/doxygen/deal.II/classFunction.html>; this is
alluded to in step-8
<https://www.dealii.org/developer/doxygen/deal.II/step_8.html#ElasticProblemassemble_system>,
but here the parameters are considered constant. For this approach you would
need to describe the material parameters by position, rather than material_id.


Indeed, this is what step-6 does. (For the Laplace equation, but if you see it 
happen there, you will know how it works for the elasticity equation as well.)


Best
 W.


--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assembly of material forces

[Wolfgang Bangerth]

Now I try to compute the material forces for the whole structure. Hence, I
tried using the same procedure of assembly for our LHS or RHS matrices in
deal.II, but it fails.


"It fails" does not carry enough information to really say what may be the 
cause. Please take a look at


https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#im-stuck



Extra question: Would my work be of interest for your code gallery maybe in
the near future?


If you have good questions, always!

Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Approximate Derivative Tensor functionality

[Wolfgang Bangerth]

On 02/15/2017 08:27 PM, Sumedh Yadav wrote:


I just had a query, when supplying <DoFHandler<2>,2,2> these template
arguments I suppose the highlighted 2 is for the rank of tensor?


Not quite. The template list for the function looks like this:

template


If yes then
as I also need to call the function for cell_grad<1,2> I should then be
supplying <DoFHandler<2>,1,2> as template arguments. But supplying 1 in place
of 2 doesn't work!


Because it refers to the dimension, not the order of the tensor. It turns out 
that you can safely let the compiler determine the order itself, depending on 
the argument you provide.


Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] A bug of visualization with VISIT

[Wolfgang Bangerth]

On 02/21/2017 09:32 PM, Jiaqi Zhang wrote:


I am not sure if this is the right place to report this.

In one of the following figures, the color plot doesn't math the contour plot,
is this the problem with VISIT?


Can you explain how the way you created the two figures differ?

Visit visualizes things by splitting each quad into two along one of the two 
diagonals, and then using a linear approximation on each triangle. These do 
not necessarily correspond to the bilinear function that you are outputting, 
and may be the difference between the pseudocolor plot and the isocontour plot 
(neither of which are equal to the P2 function in deal.II).


You can probably get a better picture by providing a positive argument to the 
build_patches() function.


Best
 W.


--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Example of using MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks?

[Wolfgang Bangerth]

On 02/17/2017 09:14 PM, Fred wrote:



I currently have a few codes written in the meshworker framework that solve
systems of PDE's. (Elasticity, Laplace/Poisson, Stokes) I am interested in
experimenting with Multigrid preconditioner for these.  I have been unable to
get this to work using the
MeshWorker::Assembler::MGMatrixLocalBlocksToGlobalBlocks function.  I am
interested in seeing an example of a code where this strategy is used.  If
there is any interest in seeing my code, I will be happy to produce a minimal
non-working example.  For now, the compilation error that I get from clang is
as follows (g++ is similar).


git/dealiiBase/dealiiInstall/include/deal.II/lac/matrix_block.h:932:10: error:
  no viable conversion from returned value of type 'const
  MGLevelObject<dealii::SparseMatrix>' to function return type 
'const
  MGLevelObject<MatrixBlock<dealii::SparseMatrix >>'
  return *matrices.read* >(i);

The line that triggers it is the following

dealii::MeshWorker::integration_loop<dim,dim>
  ( dof_handler.begin_mg(),
dof_handler.end_mg(),
dof_Info,
info_box,
LDGintegrator,
assembleSystem);

LDGintegrator is my custom integrator of a type inheriated from
MeshWorker::LocalIntegrator;
and assembleSystem is of type

MeshWorker
  ::Assembler
  ::MGMatrixLocalBlocksToGlobalBlock< SparseMatrix,
double>;


Any help would be greatly appreciated.  I would be happy to contribute any
changes required to get this to work.


Can you open an issue on the github forum at
  https://github.com/dealii/dealii/issues
It would be great to have a minimal example that shows the problem.

Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Reading user_indices from a vtk file.

[Wolfgang Bangerth]

On 02/24/2017 04:15 AM, 987093...@qq.com wrote:


In the file,I have set up user_indices for each grid cell, but when I use
GridIn to generate a mesh, there is no user indices that I have saved on each
grid cell. What should i do?


I don't actually see anything in the vtk file that would describe user 
indices. All I see are sections for points (vertices), cells, cell types, and 
a variable "rho". The latter is data defined *on a mesh*, and is ignored when 
reading the mesh.


So it's not clear to me what data you want to import. In general, though, we 
don't read user indices. We do read material_ids, however.


Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Access specific element within a distributed triangulation

[Wolfgang Bangerth]

Seyed,


Generally, in FEM you number elements like in the below picture:

<https://lh3.googleusercontent.com/-0n-HlF1EsRs/WLBEmN34aaI/AGU/GlK9R9OqDJ8FSTLKGXeaI85eU6pdV3MkwCLcB/s1600/System.png>


No, that is patently wrong. First, you are only considering rectangular 
domains subdivided into rectangular cells. Second, you think that 
because there is a "simple" numbering that that is what codes "should 
do". Both of these assumptions are wrong. How, for example, should a 
code number these cells:


https://raw.githubusercontent.com/dealii/dealii/master/doc/doxygen/images/distributed_mesh_1.png

https://raw.githubusercontent.com/dealii/dealii/master/doc/doxygen/images/hyper_shell_96_cut.png

Starting a sentence with "Generally, in FEM you number elements..." is 
generally wrong.




I am a bit confused since it is a rather easy task which should not be
such a big "deal" in deal.II.


Because it's *not* easy. It's only easy if you consider simple situation 
like yours. But if you want to do complicated things, they generally 
turn out to be complicated.


Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: FE_Q constructor

[Wolfgang Bangerth]

On 02/23/2017 08:47 AM, Aycil Cesmelioglu wrote:

I have an elasticity type problem; the horizontal displacement depends
on (x,y) but the vertical component of the displacement doesn't depend
on y. It's like it is 1.5 dimensional instead of 2. Because of this the
variational formulation also has some simplifications.


I think that's going to be a bit more complicated to implement. You can 
either try to work with both a Triangulation<1> *and* a Triangulation<2> 
(and associated finite elements and DoFHandler objects), or do 
everything in 2d and use constraints (or a penalty term) to force one 
function to be constant in y-direction.


But why care? 2d computations are so cheap these days. Can't you just 
assume that everything depends on x and y?


Best
 W.


--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] bug in program only with adaptive mesh refinement

[Wolfgang Bangerth]

Hi Dan,


* In the formula above, P(.) is a functional, I assume, i.e., it
takes a
function and returns a number, right?
* If so, what exactly does
 f(u) . v
actually mean? How do you compute this?
* Same for the second derivatives?


Sorry if that was unclear -- if only we could write LaTeX in a newsgroup
:) Yes, P is a functional, so P : H^1 to R; f : H^1 to (H^1)^*; and df :
H^1 to the space of linear operators from (H^1) to (H^1)^*. Maybe a
better way to write it would be

P(u + h * v) = P(u) + h * <f(u), v> + h^2 * <df(u)v, v> / 2 + O(h^3)

where < * , * > is the duality pairing.


...which you compute via quadrature? Or do you compute a vector F that 
corresponds to f(u) somehow?


If you go the route via vectors you have to pay attention to *what kind 
of vector you have*, namely one that does or does not incorporate 
constraints. Dealing with dual space vectors (e.g., rhs vectors) is 
tricky because they obey different rules than primal space vectors 
(e.g., solution vectors).




I assume you mean by "error" the size of the second and third term?


Yes exactly, I checked that the errors in the linear approximation
looked quadratic and that the errors in the quadratic approximation
looked cubic as a function of h.


> In addition, the
> value of the action is roughly the same for each h for both the
adaptively
> refined and uniform meshes.

I don't think I understand this statement.


By that, I meant that, if u and v are the velocity fields on the uniform
mesh and u_a, v_a are the velocity fields on the adaptively-refined
mesh, then P(u + h * v) = P(u_a + h * v_a) to within some negligible error.


So u_a, v_a are the same as u, v in the pointwise sense, just 
interpolated from a uniform mesh to an adaptively refined one obtained 
from the uniform one?




> The only other clue I have is that the error of
> the numerical solution against the analytic solution is actually
somewhat
> worse on the adaptively-refined mesh than for the uniform mesh,
but before I
> assumed this just had something to do with the linear solver.

How do you measure "worse"? As the error as a function of the number
of unknowns?


If u_true is the projection of the exact solution onto the uniform mesh
and u_a_true is the projection of the exact solution onto the
adaptively-refined mesh, then |u - u_true| is less than |u_a -
u_a_true|, where |*| is the L^2 norm of the functions (not the l^2 norm
of the vector of coefficients). This is one of the weirder things to me;
I would think that the error would be lower for the adaptive case since
overall the cells are finer. But that could just be a quirk of my
hand-rolled Newton solver.


Are you using the L2 projection? For example, for the Laplace equation, 
I think you can only guarantee that the H1 seminorm error is smaller for 
a finer mesh, but not necessarily the L2 error. It all depends on your 
equation.


And how do the two meshes differ?

Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] A bug of visualization with VISIT

[Wolfgang Bangerth]

On 02/22/2017 06:41 AM, Jiaqi ZHANG wrote:

Soluiton2 is fine. For Solution1, I checked the result with MATLAB, the
contour plot is right.
I also follow your instruction and tried data_out.build_patches(2).


Try build_patches(6) or build_patches(10).

As I stated, (i) what you output is not what is used internally in deal.II, 
and (ii) what Visit plots is not what you output. But with higher arguments to 
build_patches, you get closer in both regards.


Best
 W.


--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] FE_Q constructor

[Wolfgang Bangerth]

Aycil,


Hi, I am a new user and I want to find out whether I can use deal.ii to solve
a 2D problem where there are two unknowns; a scalar function of x,y and a
scalar function of only x. I would like to initialize these as a vector and my
test functions will have the same structure.

Is it possible to declare this using FESystem and FE_Q?


I want to use a very simple mesh

GridGenerator::hyper_rectangle (triangulation,

Point<2>(0,-0.05),

Point<2>(1,0.05));


and I thought

fe (FE_Q(1), 1,

FE_Q<dim-1, dim>(1), 1),


may work but I am not sure whether I understood dim and spacedim correctly.
Can this be done with the same triangulation for both unknowns?


No. You can only use FE_Q<N,M> on a DoFHandler<N,M> that needs to live on a 
Triangulation<N,M>.


I don't think I know what you want to do. Can you explain what PDE you're 
trying to solve and why you want to mix 2d and 1d things?


Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] bug in program only with adaptive mesh refinement

[Wolfgang Bangerth]

Daniel,


Hi all -- I'm writing a library that involves solving a nonlinear elliptic
PDE, which I'll write as f(u) = 0. There is an exact solution for this PDE for
a certain simplified geometry. To test everything, I check that my numerical
solution is tolerably close to the analytic solution, with both bilinear and
biquadratic elements, and on a mesh that has been adaptively refined in part
of the domain. This all works fine.

The PDE can be derived from an action principle, i.e. there's a convex
functional P such that f = dP. I've decided to reimplement some of my code to
explicitly use the action principle, so I wrote a routine to calculate P. The
action is analytically computable for the exact solution I'm already using to
test the code, so this is already one unit test. However, I decided to take it
a step further, and check that the nonlinear differential operator f is the
derivative of P, and that its linearization df is the Hessian operator of P.
To check that, I take some vector field u, another vector field v, and check 
that

  P(u + h * v) = P(u) + h*f(u) . v + h^2 * (df(u)v) . v / 2 + O(h^3)


> [...]

I agree it makes not very much sense, but that often points to a 
misunderstanding rather than a bug. So let me poke:


* In the formula above, P(.) is a functional, I assume, i.e., it takes a 
function and returns a number, right?

* If so, what exactly does
f(u) . v
actually mean? How do you compute this?
* Same for the second derivatives?



This test passes in the simple geometry with both bilinear and biquadratic
elements, but *not* for an adaptively refined mesh. The errors in both the
local linear approximation to the action functional and the local quadratic
approximation do not go to 0 as the increment h is reduced.


I assume you mean by "error" the size of the second and third term?



In addition, the
value of the action is roughly the same for each h for both the adaptively
refined and uniform meshes.


I don't think I understand this statement.



The only other clue I have is that the error of
the numerical solution against the analytic solution is actually somewhat
worse on the adaptively-refined mesh than for the uniform mesh, but before I
assumed this just had something to do with the linear solver.


How do you measure "worse"? As the error as a function of the number of 
unknowns?

Best
 W.

PS: I do like this approach to testing. It shows great maturity in designing 
code and testing numerical methods to think this way. Well done!



--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] second derivatives of solution

[Wolfgang Bangerth]

On 02/12/2017 04:50 AM, hanks0...@gmail.com wrote:


It seems that my question is very easy one, but I can't find the similar 
thing...

I know how I can use solution and first derivative solution that comes from
fe_values.get_function_values and fe_values.get_function_gradients

So I thought that second derivative can be obtained by

/
std::vector<Tensor<2,dim> > sol_hess(n_q_points);

fe_values.reinit (cell);
fe_values.get_function_hessians(solution, sol_hess);
...

sol_hess[q_index][0][0] //it is the same as d^2(solution)/dx^2
sol_hess[q_index][1][1] //it is the same as d^2(solution)/dy^2


However, even though the solution is continuous, the plot of
d^2(solution)/dx^2 and d^2(solution)/dy^2 look very weird.

it looks like set of delta function on the each element.

So my question is..

Is it wrong way to get second derivatives of solution??


How do you define "looks like set of delta function"? :-)

One thing to remember is that for finite elements, the solution is continuous, 
the gradient is continuous on every cell but discontinuous between cells, and 
the second derivatives is, in general, a delta function at the interfaces 
between cells (because you are taking the derivative of a discontinuous function).


That doesn't mean that the second derivatives can't be used for some purposes, 
but you need to know what you are doing.


Best
 W.


--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Recommendation: BCs in config

[Wolfgang Bangerth]

On 02/13/2017 07:09 AM, Bruno Turcksin wrote:



you could write everything into one string and then split it. So it
would like this:

subsection bc
  set value = 0.5*x+y; x*y^2
end

Take a look here
http://stackoverflow.com/questions/236129/split-a-string-in-c for some
ideas how to split a string.


Yes, this is the way to go. In fact, we have a function that splits strings:

https://www.dealii.org/8.4.0/doxygen/deal.II/namespaceUtilities.html#a8d799bb35ac16d206818c88e82afbfae

We use the scheme Bruno proposes to split strings first by semicolon, 
and then again by comma all the time in our own projects. If you don't 
like semicolons/commas, any other choice of course also works.


As far as working with people who are not used to programming -- this is 
a tricky question. It is true that you need to provide a way that's 
readable, but using work-arounds like have multiple sections is also not 
entirely obvious and will invite questions whose answers you also have 
to document somewhere. It is a balance that's sometimes difficult to 
strike. I typically strive for a solution that is easiest to describe 
and easiest to understand *in its entirety*. In your case, I would 
probably choose splitting on semicolons and colons, and in the 
documentation show how you can split different, semicolon-separated 
parts of the formula onto separate lines using a backslash at the end of 
the line to show that that makes formulas easier to read.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Parameter map default value

[Wolfgang Bangerth]

On 02/13/2017 06:15 AM, Franco Milicchio wrote:

Hello,

I know this is a weird question, but is there a way of telling a *Map*
that there is a default value?

For instance I'd like to express something like "*stress:4, power*"
instead of what I have to do right now, such as "*stress:4, power:8*"
(or whatever number).

The following is the code I am using now:

params->declare_entry("csv", "strain:3, work:1",
dealii::Patterns::Map(dealii::Patterns::Selection("time|displacement|stress|strain|work|power"),
dealii::Patterns::Integer(0), 1), "Dump to CSV file, comma separated
values (time, displacement, stress, strain, work, power)");


Not a weird question. But also not something that's implemented.

I've been thinking for a minute on how to implement this, but the answer 
is that it's not going to be easy because the patterns are really only 
used to describe if a particular value for a parameter *satisfies 
certain requirements*, but the patterns have nothing to do with the code 
that later on parses and interprets the value.


So one thing you can do is to describe the parameter as
  Patterns::Map(Patterns::Selection("..."), Patterns::Anything(), 1)
which works because Anything also matches the empty string (along with 
anything else, of course). You would then just have to ensure later on 
that what is in this spot is either the empty string (which you can 
interpret as a default value), or indeed a positive integer.


I'm going to note that figuring out that a string is indeed a positive 
integer (and nothing else) is not trivial.


Best
 W.


--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Query using hdf5 with deal.ii

[Wolfgang Bangerth]

Rajat,


I am using deal.ii to solve a 3D solid mechanics problem. I am using
Petsc and P4est along with deal.ii

There is a previous post
<https://groups.google.com/forum/#!searchin/dealii/hdf5%7Csort:relevance/dealii/_51WLg4jDGI/yvFDxSLeTswJ>
 related
to my query but it does not talk in detail about using hdf5.

My code currently follows step 40 and outputs one vtu file from each
processor and one pvtu file from a rank 0 processor.
However, as the number of processors grow, this is generating huge
number of files as output which inturn takes a lot of time to copy /
move to my desktop for postprocessing.
Also, if I output data from multiple equations (being solved
concurrently), the mesh data is replicated in each of these files.

I came across using hdf5 which might solve this issue. Since there is no
example that demonstrates the use of hdf5, I have a a question before
I start exploring and using it. Will this output just 1 file independent
of the number of processors on which the code is run ?

This is asource file of aspect
<https://github.com/geodynamics/aspect/blob/6db2fa985d24f7d4e60dc8cde9b9a2d63817113b/source/particle/output/hdf5.cc>
that demonstrates using hdf5 in detail.


This is the only example I know of.

But you may want to look at the VTU part there again: it allows writing 
the data from whole *groups* of processors into one file so that you 
don't end up with one file per processor. We typically use this for 
larger simulations.




Also, for snapshot creation and restart purposes, I output 1 file from
each processor containing gauss point data at a current time step. This
is also not scalable for the same reasons.
Is there any way that I can get only 1 file (or maximum 1 file per node)
that contains the data from all the gauss points spread over multiple
processors ? Or can this data be included while
doing triangulation.save() which saves triangulation and data that is
associated with it using Parallel solution transfer class ?


Both is possible. You can attach the data to the triangulation. ASPECT 
again does this with particle data, for example -- I don't know where 
exactly this happens, but you can take a look at the particles directory.


It is also possible to let each processor collect its own data, put 
things into one long array, send that to processor 0, and let that 
processor put it into a central file. I believe that that, too, is what 
ASPECT does. You may want to look into the 
source/simulator/checkpoint_restart.cc file. (I'm going to say that 
we've implemented these pieces so long ago that I don't recall all of 
the details. So I may be wrong in my claims above, but that's where I 
would look things up.)


Best
 W.



--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] second derivatives of solution

[Wolfgang Bangerth]

On 02/12/2017 06:15 PM, hanks0...@gmail.com wrote:


Do you mean that it is right way to get second derivatives but it could not be
continuous?


Yes, you are using the correct way of computing something that is of 
questionable usefulness. As I mentioned, this is because for finite elements, 
the function is continuous, the derivative is discontinuous, and then the 
second derivative is questionable.


Think of the situation just in 1d where the solution is a piecewise linear 
function. Imagine how the derivative looks (piecewise constant) and how the 
second derivative would look like (a set of delta functions).




I also attach the plot for d^2(solution)/dx^2 that look weird to me...


I think it is plausible. It is at least not obviously wrong in my mind.

Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Assemble Righthand Side for vector-valued problem

[Wolfgang Bangerth]

On 02/09/2017 04:45 AM, Jaekwang Kim wrote:


The manufactured solution of velocities has not satisfies the continuity
equation.
(i.e. manufactured solution does not satisfy, div.u=0)...


Yes, that would do it :-)
Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Nonhomogeneous Dirichlet Boundary conditions using a Dirichlet lift

[Wolfgang Bangerth]

as far as I have understood (but I might be wrong), the functions
VectorTools::interpolate_boundary_values
<https://www.dealii.org/8.4.0/doxygen/deal.II/namespaceVectorTools.html#af6f700f193e9d5b52e9efe55e9b872d5>
MatrixTools::apply_boundary_values
<https://www.dealii.org/8.4.0/doxygen/deal.II/namespaceMatrixTools.html#a41a069894610445f84840d712d4f891e>
find the nodes where Dirichlet BC's are applied and then there impose the
corrensponding boundary value, after having built the system matrix and
right-hand side.

Another possibility would be to use a Dirichlet lift, change the weak
formulation and solve for homogeneous Dirichlet boundary conditions. I am
wondering if someone already did this or if it somewhere implemented in deal.ii


It may not look like it, but that's really what the functions do that you cite 
above.


The algorithm is a bit complicated, but take a look at lectures 21.6 and 21.65 
here:

  http://www.math.colostate.edu/~bangerth/videos.html

Best
 Wolfgang

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Efficient implementation of varying Dirichlet BCs when everything else in the system is constant

[Wolfgang Bangerth]

On 08/19/2016 04:23 PM, David F wrote:


P.S.: I tried to simply repeat the process to re-apply the BCs with a new
value using the same matrix, but I get always the first BCs I apply  (I think
because once I apply it, the matrix knows it has been already condensed and
ignores the next calls).


Yes. What you need to do is in every time step copy the matrix you have 
previously assembled into a new matrix, call apply_boundary_values() to this 
matrix and your current right hand side, and then solve the linear system. 
This way, you always start from the same matrix.


I think one of steps 23, 24, 25, 26 does this.

Best
 Wolfgang

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] weakly enforcing an identity - vector mean curvature as the surface laplacian of the identity on a codim-1 manifold

[Wolfgang Bangerth]

On 08/19/2016 01:15 PM, thomas stephens wrote:

Wolfgang, I appreciate your assistance so far.

I'm getting closer, but I still need some help.  I would like to look at this
problem as it appears in [1], Equation 27 on p. 11, as the bilinear forms are
written as integrals, making things more clear.

Below I have set up what I think to be the system for (k_{bar,h})_X  =  -
(nabla_X v_h, nabla_X id_X)_X - which, in light of Eqn 27 in [1], is really

\integral  k_bar_h \dot v_h = \integral nabla_X id_X : nabla_X v_h


|
 typename DoFHandler<dim,spacedim>::active_cell_iterator cell =
dof_handler.begin_active();
  for (; cell!=dof_handler.end(); ++cell)
  {
cell_mass_matrix = 0;
cell_rhs = 0;
fe_values.reinit (cell);

// create lhs mass matrix, (k*nu, v)_ij
for (unsigned int i=0; i<dofs_per_cell; ++i)
  for (unsigned int j=0; j<dofs_per_cell; ++j)
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
  cell_mass_matrix(i,j) += fe_values.shape_value(i,q_point) *
   fe_values.shape_value(j,q_point) *
   fe_values.JxW(q_point);

// create rhs vector, (nabla_X id_X, nabla_X v)_i := \int nabla_X id_X :
nabla_X v
for (unsigned int i=0; i<dofs_per_cell; ++i)
  for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_rhs(i) += (fe_values.shape_grad_component(i,q_point,0)*

identity_on_manifold.shape_grad_component(fe_values.normal_vector(q_point),0) +
fe_values.shape_grad_component(i,q_point,1)*

identity_on_manifold.shape_grad_component(fe_values.normal_vector(q_point),1) +
fe_values.shape_grad_component(i,q_point,2)*

identity_on_manifold.shape_grad_component(fe_values.normal_vector(q_point),2))*
fe_values.JxW(q_point);

cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
  system_rhs(local_dof_indices[i]) += cell_rhs(i);

  for (unsigned int j=0; j<dofs_per_cell; ++j)
mass_matrix.add (local_dof_indices[i],
 local_dof_indices[j],
 cell_mass_matrix(i,j));
}
  }

|

Already I am confused by this, since the double contraction : (which is
implemented long-hand above) ultimately puts a scalar into cell_rhs(i),
whereas I think I am looking for each element of cell_rhs(i) to have
spacedim=3 components.  I am quite confused about this actually.  (The reason
why I think elements at cell_rhs(i) should be vector-valued is because each
element in k_bar is vector-valued, so in mass_matrix*k_bar = system_rhs, the
components of system_rhs should be vectors.


No, that doesn't make any sense. For each i, the weak form of the equation 
leads to just a scalar equation on both the left and right hand side. Take a 
look how this works for other vector-valued problems in the documentation 
module on vector-valued problems.


In your code above, it is not clear to me what 'identity_in_manifold' is. Does 
this simply represent X? If so, aren't you just computing

   grad phi_i * grad X
? You need to somehow use the grad_Gamma instead, though.



If the above is correct, then I need help getting the mean curvature vector
out of this system.  If it is the case that I have the pieces to write
mass_matrix*mean_curvature_vector = system_rhs, and I need to invert
mass_matrix, then how to extract mean_curvature_vector?  The following code
yields all zeros for mean_curvature_vector:

|
/// solve  mean_curvature_vector = inverse_mass_matrix * system_rhs
  SolverControl solver_control (mean_curvature_vector.size(), 1e-7 );
  SolverCG<> cg_solver (solver_control);

  const auto inverse_mass_matrix = 
inverse_operator(linear_operator(mass_matrix),
cg_solver,
PreconditionIdentity());
  mean_curvature_vector.reinit(dof_handler.n_dofs());
  mean_curvature_squared.reinit(dof_handler.n_dofs());
  system_rhs *= -1.0;
  inverse_mass_matrix.vmult(mean_curvature_vector,system_rhs);


  std::cout << "mean curvature vector" << mean_curvature_vector << std::endl;
  mean_curvature_squared = mean_curvature_vector*mean_curvature_vector;
  std::cout << "mean curvature squared " << mean_curvature_squared << std::endl;

|


That is, I do not think I am handling these data structures correctly -
perhaps that LinearOperator inverse_mass_matrix is not just an ordinary
dealii::Matrix that I can then apply to dealii::Vectors.


Well, is system_rhs nonzero?

Best
  W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For maili

[deal.II] Are you funded to write applications based on deal.II?

[Wolfgang Bangerth]

All,
Timo and I are writing a proposal that will (hopefully) fund further 
development of deal.II in the coming years. It would be helpful if we had an 
overview of what other projects there are where people are funded to write 
applications based on deal.II (whether or not they will be made available to 
anyone else is not important).


So: If you or your adviser has a grant that specifically mentions that 
software based on deal.II will be developed, please drop me an email (not as a 
reply to the list, though).


Thanks
 Wolfgang

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Does MappingFEField apply to all cells ?

[Wolfgang Bangerth]

On 08/22/2016 07:30 AM, Praveen C wrote:

Thanks Wolfgang. I should have checked the built_patches function. It works
fine if I give the third argument. But now I have to find the real reason why
my code is crashing with MappingFEField :-(


There is always one more bug. Always.
Best
 W.


--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] write_vtk fail when using hp functionality

[Wolfgang Bangerth]

On 02/27/2017 12:24 PM, Weixiong Zheng wrote:

Dear all,

Might be a trivial questions to you guys. I get a scalar problem using hp
class. When I tried to output the results using

|
  std::vectorsolution_names;
|

  solution_names.push_back ("phi");



  std::vector

  data_component_interpretation(1,
DataComponentInterpretation::component_is_scalar);



  DataOut<dim,hp::DoFHandler > data_out;



  data_out.attach_dof_handler (dof_handler);

  data_out.add_data_vector (scalar_flux,

solution_names,

DataOut<dim,hp::DoFHandler >::type_dof_data,

data_component_interpretation);



  data_out.build_patches ();

  std::string name = "solution.vtk";

  std::ofstream output (names.c_str());

  data_out.write_vtk (output);

  std::cout << "see see 3" << std::endl;


I got the following error:

CDLS-DSA(49490,0x7c0cb000) malloc: *** error for object 0x7f7f87707c98:
incorrect checksum for freed object - object was probably modified after being
freed.

*** set a breakpoint in malloc_error_break to debug


Weixiong,
Can you run this in a debugger, set the breakpoint at the function suggested 
here, and see if you can get a backtrace? Alternatively, can you create a 
small, self-contained testcase that shows the problem?


Best
 W.


--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Configuring deal.ii error - can't pinpoint what the issue is

[Wolfgang Bangerth]

Phil,


I'm a new user trying to ultimately install ASPECT, currently setting up 
deal.ii.

I've managed to almost get deal.ii to configure, but errors have occurred.


Which errors specifically do you see? detailed.log doesn't show anything 
obvious, but I may also be looking in the wrong place.


Best
 W.


--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] weakly enforcing an identity - vector mean curvature as the surface laplacian of the identity on a codim-1 manifold

[Wolfgang Bangerth]

Thomas,


I was able to solve for the vector mean curvature using the weak form of
the identity k_bar = laplacian_X id_X, where X is a codimension 1
manifold without boundary.  The image above is a plot of the square of
the mean curvature on an ellipsoid with semi principle axes a=1,b=2,c=3.


Great, this looks roughly correct I would say. Have you checked that the 
numerical values are correct as well?


I think this would actually make for a really nice addition to the set 
of tutorials, or the code gallery at

  http://dealii.org/code-gallery.html
Would you be interested in getting it there? (Maybe after rewriting it 
in such a way that you solve for all three components at once?)


Cheers
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Installation/Tutorial 1 Issue

[Wolfgang Bangerth]

In the case of trying to compile the tbb.cc code using the command
sudo gcc tbb.cc -std=c++11

I get a long string of errors that look something like this


tbb.cc:(.text._ZN9__gnu_cxx13new_allocatorISt10_List_nodeIN6dealii10WorkStream8internal15Implementation225IteratorRangeToItemStreamINS_17__normal_iteratorIPiSt6vectorIiSaIi12scratch_data9copy_dataE8ItemType17ScratchDataObjectEEE8allocateEmPKv[_ZN9__gnu_cxx13new_allocatorISt10_List_nodeIN6dealii10WorkStream8internal15Implementation225IteratorRangeToItemStreamINS_17__normal_iteratorIPiSt6vectorIiSaIi12scratch_data9copy_dataE8ItemType17ScratchDataObjectEEE8allocateEmPKv]+0x46):
undefined reference to `operator new(unsigned long)'
/tmp/ccAgCMHw.o: In function


Peter -- two comments:

1/ Using sudo is almost always a terrible idea. What you're doing is 
calling the compiler as *root* (i.e., as the superuser). There really 
shouldn't be any reason to do this, but so many reasons *not* to do this 
because you may accidentally overwrite files that shouldn't be 
overwritten, you could be writing into directories that need to stay as 
they are, etc.


2/ The command 'gcc tbb.cc -std=c++11' trie to compile and link the 
tbb.cc file into an executable. But tbb.cc depends on many other files, 
so while you should be able to compile it into an object file, you won't 
be able to link it into an executable without linking in a large number 
of other object files. That's in essence what the error message is saying.



My suggestion is to just stick with how the installation process is 
intended to be used:

* no sudo, no superuser
* compile and install into a directory somewhere within your homedir
* go to one of the example dirs in the installation directory
* call 'cmake .'
* compile
This should work.

For your own project, you can copy one of the existing example 
directories somewhere else, then call

  cmake -DDEAL_II_DIR=/where/you/installed/dealii .
and it should still work.

Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Creating my own constraints matrix for applying dirichlet bc

[Wolfgang Bangerth]

On 09/05/2016 09:32 AM, Praveen C wrote:

So if I loop over all non-artificial cells instead of locally owned cells to
build the dirichlet map, then it seems to work correctly.


Yes, that is the correct solution -- you need all locally owned or ghost cells.

Best
 W.
--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] merging geometries(geo files) in gmsh for computational geometry definition

[Wolfgang Bangerth]

Marek,


I have prepared slightly better mesh - without "doubling" in the intersection,
see attachment.
However I got error on reading in the mesh file.

n error occurred in line <1447> of file

in function
void dealii::GridIn<dim, spacedim>::read_msh(std::istream&) [with int dim
= 3; int spacedim = 3; std::istream = std::basic_istream]
The violated condition was:
false
The name and call sequence of the exception was:
ExcInvalidVertexIndex(cell, subcelldata.boundary_lines.back().vertices[i])
Additional Information:
While creating cell 55, you are referencing a vertex with index 0 but no
vertex with this index has been described in the input file.


Do You have any advice, how to cope with this problem.


The error message is correct. If you look at the file, it starts with

$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
9708
1 0 0 0
2 1 0 0
3 0 1 0
...

That is, it starts numbering vertices at 1 (first column). But then, later on, 
it describes cells as follows:


$Elements
11703
1 15 2 0 1 1
2 15 2 0 2 2
3 15 2 0 3 3
...
56 1 2 0 5 0 52
...

Here, each line corresponds to an point, edge, face, or cell with properties. 
Column 1 has the number of the object, column 2 the type of the object 
(line=1, quadrangle=2, hexahedron=5, point=15). After the type is a number 
that indicates how many properties there are for this object (here, 2 
properties for points and lines), followed by this many property values, 
followed by the vertex indices. Indeed, object 56 (a line, with two properties 
with values 0 and 5) consists of vertices 0 and 52. But we don't have a vertex 
with number 0, since they were numbered starting at 1.


I don't know how this mesh was created, but the content of the file indeed 
does not seem correct to me. You probably need to talk to the people who wrote 
gmsh what this is supposed to mean.


Best
 Wolfgang


--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Doubt regarding constraints

[Wolfgang Bangerth]

Suppose I want to constraint dof1 owned by processor 1 to be equal to dof2
owned by processor 2.


I can do that by calling


|
ConstraintMatrixconstraints;
constraints.reinit(locally_relevant_dofs);
constraints.add_line(dof1);
constraints.add_entry(dof2);
|

Q1. Should it only be called on one processor ?


If one of dof1 or dof2 are among the locally relevant ones (i.e., on active or 
ghost cells), then both processors need to know about this constraint.



Q2. My question is, what happens when this is called on both the processors
(processor 1 and processor 2) ?


Nothing bad. In fact, both need to know about it.



Q3. What happens when the processor 2 calls the same thing but in the other
order ?


Bad things. Don't do it :-)

Best
 Wolfgang

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] modify solution vector

[Wolfgang Bangerth]

On 09/02/2016 06:07 AM, Bryukhanov Ilya wrote:


I want to change the solution vector (for example displacement vector from
step-8 and step-18) by the vector
that depends on the point. I want to iterate over all vertexes and add to each
solution component some value
that is a function of the vertex coordinate.
How can I do it in dealii?


It sounds like you have a solution u_h(x) and you want to add to it a function 
v_h(x) that you know at each vertex location, i.e., it *interpolates* a 
continuous function v(x) that you know.


If put this way, I would suggest you call
  VectorTools::interpolate
on v(x) to get the nodal values V of v_h(x), and then add V to your solution 
vector U.


I think that ought to be easier than what you suggest (though internally it 
does exactly that).


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: applying 3D mesh in step-20

[Wolfgang Bangerth]

On 09/04/2016 05:38 AM, hanks0...@gmail.com wrote:


I have to show my result to my boss(variable p). And, To show the inside of
torus, It seems that I have to use contour plot using Visit.(In the Wolfgang
lecture, he used contour plot to see the inside of the calculated result in 3D
problem)

 But, When I used  Contour plot in Visit for variable p, It is said that " The
Contour plot of variable "p" yielded no data."


That's because you use piecewise constant pressures. Any given isosurface 
values does not typically equal the value the pressure has on any of the cells.




Is there anyway for me to show the inside of the calculated result to my boss?


Use cutplanes, cutspheres, thresholds, etc., to cut away some cells while 
showing the others.


Best
 Wolfgang

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Applying Multi Point Constraints

[Wolfgang Bangerth]

Yes you are right. It turns out this is the problem.
I was creating a hyper_rectangle using grid generation class. Then 4 and 5
corresponded to the -/+ z surface.
(I used hard coded numbers as it saved time for looping over all the elements.)

>
> Now when I am using the same mesh but now reading via read_ucd(), upon reading
> the documentation carefully, I saw that the ordering is changed. So 4 and 5
> now no more correspond to -/+ z surface.

It's always a safe bet in debugging other people's codes to challenge their 
underlying assumptions :-)




However, If I call cell->face(4)->unit_cell_direction, it still gives output
as 2, which means that this face is normal to z surface.
This normal direction is for the cell in the natural coordinates, right ?


In the "reference coordinate system" (which we often call the "unit cell 
coordinate system". But that's not the coordinate system you care about.




So, in general, how should I know if the real cell face is normal to -/+ z
direction ?


You could query the z-coordinate of a face's center. Loop over all faces of a 
cell, if you find one whose center's z-coordinate is at the top or bottom of 
the domain, then you have a cell you care about.




Will having a constraint such that both the dofs are locally relevant & not
locally owned result in anything bad ? I think, it should have no effect
atleast for that processor on which both these two dofs lie.
Am I right here ?


Correct.
W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Help with Vector template instantiations

[Wolfgang Bangerth]

Dragan,


It is ok, you can apply patch yourself.


I've committed under your name (which hopefully I got right):
  https://github.com/dealii/dealii/pull/3114



I actually also had to specialize many template functions to work with the new
type (those related to matrices and vectors, BlockMatrixBase::add,
SparseMatrix::add and also local_apply_boundary_conditions etc.). There are
some functions that do not work well with the user-defined types like
Vector<>::allocate(), Vector<>::deallocate() and Vector::operator
=(). They cause seg. faults. I am not sure why, perhaps because of mem.align
you use (Utilities::System::posix_memalign) and it is fine for primitive
types. Anyway that is something internal for my case.


Hm, good question what is happening there. Feel free to propose individual 
patches for each of these cases if you can identify an underlying reason.




I am finishing what I started doing and will write you about an interesting
application of deal.II. It is mostly useful for multi scale modelling in
chemical engineering but is general in nature. It is coupling of deal.II with
an equation-based simulator: basically, using deal.II only to assemble
matrices/rhs (including non-linear FE cases), generating equations and then
solving one or more of these systems together with other differential and
algebraic equations in a large DAE system. All that is done from python
(although implemented in c++). An example could be Lithium-ion batteries:
electrodes are typically made out of a porous material composed of large
numbers of solid particles, and particles are in a electrolyte. Particles
(transport modelled using the FE method) are coupled at the larger electrode
length scale via electrolyte transport (using the finite volume method). Two
software are fully coupled and the boundary conditions can be set using the
equations from the other software etc.


That would make for a fantastic code gallery project!

Best
 Wolfgang

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] How to declare a constant that is used in a template?

[Wolfgang Bangerth]

On 09/10/2016 03:49 AM, dealii.gr...@gmail.com wrote:


I have a dimension independent code which can work for 1D, 2D and 3D. The
dimension of the problem is defined in the input file. Is there a way to
declare a constant integer variable (dim) for template?

for 1D
dim = 1;
for 2D
dim = 2;
else
dim = 3;

then,

FEM solid("parameters.prm");
solid.run();


The problem is that the  template parameter is something the compiler 
needs to know at compile time, but you want to only determine it at run time. 
This doesn't easily work.


But you can do this with tricks -- take a look at the main() function of 
ASPECT, for example:

 https://github.com/geodynamics/aspect/blob/master/source/main.cc
lines 452-490.

Best
 Wolfgang

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Problem creating a pull request on code-gallery

[Wolfgang Bangerth]

Hi Manuel,


I want to contribute with an small project to the code-gallery.


Fantastic!


> I followed the

instructions at (http://www.dealii.org/code-gallery.html) but I think I might
have had a problem on the very last step; i.e., creating a "pull request". I
went to my github account and clicked on "New pull request" but got


  There isn’t anything to compare.

*dealii:master* and *manuel-quezada:master* are identical.


This looks like you are trying to create a pull request from your deal.II fork 
against the deal.II master repository. But you want to create a pull request 
from your code-gallery fork against the code-gallery master repository.


Does this help?

Cheers
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

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

[Wolfgang Bangerth]

On 09/08/2016 10:17 PM, hanks0...@gmail.com wrote:

Thanks for your answer. But, I still don't know how I can deal with this error.


Then you need to describe more clearly what you want to do. As I mentioned 
here...



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


...your original formula for K_Inv[0][0] is so that it can lead to 1/0 for 
some values of x,y. There is simply nothing you can do about this -- it's not 
a question how you *compute* it, it's that that's what the *model* says. If 
you don't like it, you need to start thinking about whether that's the 
*correct model*!


Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Help with Vector template instantiations

[Wolfgang Bangerth]

On 09/10/2016 11:22 AM, dragniko...@gmail.com wrote:

Yes, no problem. I see those files changed a lot since the 8.4.1 release.
Anyway, I see the two spots are the same as in 8.4.1. I'll send you the patch.


Fantastic, thanks!
W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Doubt with data output in parallel

[Wolfgang Bangerth]

My program uses Parallel::distributed::triangulation.

If I want to output a data array associated with cell, which one of the
following is the correct way of doing it ? And why ?
(This cell data type vector is to be viewed in Paraview.)

|
intindex =0
for(;cell!=endc,++cell,++index)
if(cell->is_locally_owned)
data_array(index)=value;
|


or


|




for(;cell!=endc,++cell)
if(cell->is_locally_owned())
data_array(cell->active_cell_index())=value;
|


Both will work. They are identical.

Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Help with Vector template instantiations

[Wolfgang Bangerth]

Dragan,


Here is the patch. I added one more change in functions Vector::reinit:
   if (omit_zeroing_entries == false)
*this = static_cast(0);
The line above fails if the Number type is not primitive. If it is a class it
can't cast. So it could be:
   if (omit_zeroing_entries == false)
*this = Number(0);
Therefore it will call a constructor or simply create a primitive type. I do
not know if there are some ramifications of changing that line.
Please have a look.


Yes, all of this looks correct. I think the *declaration* of these explicit 
specializations should be moved to the bottom of the .h file, though.


Do you know how to use github? This would make sure that you get credit for 
the patch, as you should! If you want to learn how to use it, take a look at 
lecture 32.8 here:

  http://www.math.colostate.edu/~bangerth/videos.html
If you want me to apply the patch, just say so, and that's ok with me as well.

Best
 Wolfgang

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] After the calculation, Can the result be used??

[Wolfgang Bangerth]

Hi Kyusik,


I just wonder if it is possible that The result of the calculation can be used
for the other calculation.

For example, In step5, The Eq is -div(alp(x)*grad(u(x)))=f with appropriate
Boundary Condition.

After this calculation, is it possible for me to use the solution that is u,
for the other values?(for example grad(u), |grad(u)|^2/x etc...)

If it is possible, Could you please explain to me...?


There are many ways how one can use a solution. But your question is too vague 
because you don't describe what exactly you want to do with the solution. For 
example, it may be that you want to use it as a coefficient for another 
equation. Or you may want to create graphical output. Or you may want to 
compute some output quantity.


All of these are possible, but depending on what you want to do, the 
approaches will be very different. Consequently, your question will have to be 
more specific.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Help with Vector template instantiations

[Wolfgang Bangerth]

On 09/10/2016 10:33 AM, dragniko...@gmail.com wrote:

Hello Wolfgang,

Yes, I do explicitly instantiate these templates in my code. The datatype is
different from what is already defined in deal.II. Something like:
#include 
#include 
template class Vector;

But these two guys:
dealii::Vector<std::complex >::operator=(std::complex)
and
dealii::Vector::lp_norm(int) const
are explicitly instantiated in vector.templates.h (lines 855 and 1476).


Ah, I see. Good point. Yes, moving these two functions into the .cc file and 
leaving a *declaration* of the explicit instantiation in the .h file is the 
right solution.


Want to try writing a patch for this? We'd be happy to merge it!

Best
 Wolfgang

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

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

[Wolfgang Bangerth]

Metin,


I need help to figure out some functions in fe_values.h and fe_values.cc
sources.

Could you please help me understand what's happening in the divergence
function (SymmetricTensor<2, dim, spacedim>::divergence) and the
symmetric_gradient function for the vectors
(Vector<dim,spacedim>::symmetric_gradient), or point out any material
(manual, video etc.)?


Can you be more specific which functions and which parts within them you 
have trouble with?


Best
 Wolfgang

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

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

[Wolfgang Bangerth]

On 09/28/2016 11:32 PM, Erik Svensson wrote:



I’m trying to implement homogeneous Neumann bc on part of the boundary in
tutorial example step-20.

The Neumann bc should be imposed strongly (as for Dirichlet bc in normal,
non-mixed, formulation). Is this correct?

In any case, I’m trying to impose the Neumann bc strongly using the technique
from step-22.

Running program I get the error:

---
An error occurred in line <1878> of file

in function
void
dealii::VectorTools::{anonymous}::do_interpolate_boundary_values(const
M_or_MC&, const
DoFHandlerType&, const typename dealii::FunctionMap::type&, std::map&, const
dealii::ComponentMask&, dealii::internal::int2type) [with DoFHandlerType
= dealii::DoFHandler<2>; M_or_MC = dealii::Mapping; int dim_ = 2; typename
dealii::FunctionMap::type =
std::map*, std::less, std::allocator<std::pair*> > >]
The violated condition was:
cell->get_fe().is_primitive (i)
The name and call sequence of the exception was:
ExcMessage ("This function can only deal with requested boundary " "values
that correspond to primitive (scalar) base " "elements")
Additional Information:
This function can only deal with requested boundary values that correspond to
primitive (scalar) base elements
---

Can someone please help me to interpret the error message, or even better
guide me how to proceed implementing the Neumann bc?


Erik -- imposing zero Neumann boundary values for the mixed Laplace implies 
that the normal velocity is zero on the boundary. This is no problem if you 
discretize the velocity variable through the ordinary Lagrange elements, but 
in step-20 we use the div-conforming Raviart-Thomas element for which this 
function is not equipped to handle the situation because the Raviart-Thomas 
element mixes both the x- and y-components of the velocity vector.


I think you should be able to use
  VectorTools::project_boundary_values_div_conforming
instead to do what you want to do.

Best
 W.


--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Issue with Trilinos installation (issue with blas?)

[Wolfgang Bangerth]

On 10/04/2016 07:12 PM, Matt Weller wrote:


Okay I tried as you suggested, and same issue. It can't
find -l/usr/lib64/libblas.so.3.6.1 . There is something off with the library
but I have not the foggiest what it is.


I don't know what to suggest other than to play with your package manager 
until you figure this out. At least you now have a small testcase.


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Running deal.ii with METIS

[Wolfgang Bangerth]

On 10/05/2016 08:41 AM, Bruno Turcksin wrote:

Do you know whether automatic
> partitioning uses METIS? But still, I think the old step-17 and 18 should
> work well with version 8.3...

Automatic partitioning does not use METIS and does not work across
nodes.


I believe we only do the automatic partitioning in step-18. But in any case, 
even the automatic partitioning relies on METIS if you run it in parallel.


Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Using p-refinement with high order elements

[Wolfgang Bangerth]

On 10/05/2016 09:03 AM, Deepak Gupta wrote:



Since it works for the combination of p=1 and 2, I do not think hanging nodes
are a problem. I looked at the hp-paper as well and do not see something that
might be going wrong related to the hanging nodes. I am still trying to narrow
down the problem and come up with a simplified example, however, I still
thought of raising this question. Based on the above scenario, can someone
figure out what might be going wrong?


I am now aware of anything that would cause this, but it certainly sounds like 
a bug. A small testcase would definitely be useful!


In the past, I have debugged problems such as this by interpolating a known 
and simple function onto the mesh, and then applying the constraints via 
ConstraintMatrix::distribute() to the vector with the interpolant. If the 
constraints are wrong, you will see this if you just output this as a field 
(with sufficiently high n_subdivisions passed to DataOut::build_patches()).


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Using p-refinement with high order elements

[Wolfgang Bangerth]

Deepak,


By the way, I am constructing the cell matrices using B'DB, where B is
the strain-displacement matrix and D is the constitutive relation. Then
I add it to the system matrix using distribute_local_to_gobal function.
Since B and D are the same for global and adaptive case, I do not expect
any difference in the cell matrices. The system_rhs is correct since I
have only a point load and I can check it by printing the vector. Thus,
the possible issue can be only in the assembly (for the local refinement
one). My problem is vector-values. In case this information can lead to
any other suggestion from someone, it would be great.


I don't have any experience with this approach, but at the end of the 
day, you get a local matrix  A_K = B_K' D_K B_K  for each cell K, and 
that should be all that matters. I assume you mean 
ConstraintMatrix::distribute_local_to_gobal, not 
cell->distribute_local_to_gobal, right?




Meanwhile I am going to extend the simple example for a basic version of
the elastic problem I am trying to solve. Hope then I can figure out the
error.


Good plan.
Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

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

[Wolfgang Bangerth]

On 10/04/2016 09:12 AM, hanks0...@gmail.com wrote:


template 
void Step6::get_profile ()
{
  const QGauss  quadrature_formula(3);
  FEValues fe_values (fe, quadrature_formula,
   update_values|  update_gradients |
   update_quadrature_points  |  update_JxW_values);
  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
  const unsigned int   n_q_points= quadrature_formula.size();

  std::vector local_dof_indices (dofs_per_cell);


  typename DoFHandler::active_cell_iterator
  cell = dof_handler.begin_active(),
  endc = dof_handler.end();
  f.reinit (dof_handler.n_dofs());//f is already declared as private
like solution

  for (; cell!=endc; ++cell)
{
  fe_values.reinit (cell);

  for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
{
   const Point =fe_values.quadrature_point(q_index);

f[q_index]=p(0)*solution[q_index]+p(1);


'solution' is a global vector of the values for each degree of freedom. 
But you index it with the quadrature point number. That makes no sense.


What you need to do is evaluate the solution at the quadrature points. 
This is done using the

  FEValues::get_function_values()
function. Several example programs use it, so you will find use cases there.

Best
 W.

--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

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

[Wolfgang Bangerth]

On 09/20/2016 09:04 AM, JAEKWANG KIM wrote:

yes...I tried that before, but I get an error message as follow..


As a general workflow, you need to
1/ make the code compile
2/ run it without any runtime error
3/ make sure the results are correct.

In your original message, you were stopped in 1 because the compiler 
complained about Solution<dim+1> where it should have been 
Solution. You then fixed this, and got a run time error, i.e., you 
were stopped in step 2.


This is *progress*. But you saw an error and instead went back to 
Solution<dim+1>. This makes no sense. A program typically has many bugs. 
You need to work on them one after the other, until you pass step 3. It 
may mean that you need to fix multiple compile errors first before you 
can run the code, and fix multiple runtime errors before you get any 
solution.


Best
 W.


--
----
Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Why is this asserted to be impossible?

[Wolfgang Bangerth]

On 09/17/2016 05:20 AM, Alex Zimmerman wrote:


Not understanding why this is impossible, I began digging through the code a
bit and thought I would write my own 1D implementation. I noticed that the
code in the existing template should work correctly without any changes. So
instead I simply copied and renamed the template, without the implementation
asserting impossibility in 1D ( Assert(false, ExcImpossibleInDim(1)); ). My
first test ran exactly as expected.

deal.II often asserts that something doesn't make sense or is impossible in a
certain dimension. Sometimes I agree and sometimes I don't. This is the first
time that it's just seemed completely wrong to me. What am I missing?


Quite possibly nothing -- we may simply be wrong. Can you point out the file 
and line where the offending Assert is located?


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Why is this asserted to be impossible?

[Wolfgang Bangerth]

On 09/18/2016 03:45 AM, Alex Zimmerman wrote:

I question all of the asserts in this
file: 
https://github.com/dealii/dealii/blob/master/source/numerics/vector_tools_rhs.cc

My particular case is the one on line 36, for dim = 1 and spacedim = 1.

Thanks for looking into this.


This looks like we disabled these functions back in the day when computing 
boundary integrals in 1d did not work. I think this ought to work now. What 
happens if you just remove these 4 functions from the file, recompile the 
deal.II library, and test whether the result works for your code?


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

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

[Wolfgang Bangerth]

On 09/18/2016 10:01 AM, 'Franck Kalala' via deal.II User Group wrote:


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


Do you actually need to form the product matrix? Or are you only interested in 
*applying* this matrix to vectors?


In the latter case, the operation
  diag(V) * vec
can be achieved using the Vector::scale() function.

Best
 W.

--

Wolfgang Bangerth  email: bange...@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.