Re: [deal.II] Periodic Boundary conditions on Step-3 tutorial with dim != spacedim

[Timo Heister]

Hi Malhar,

if you look at
https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/source/dofs/dof_tools_constraints.inst.in#L69

you can see that it is only instantiated for
DoFHandler
meaning only for dim=spacedim (the second argument is the default right now).

Instead, we would need to move the instantiation down into a loop that
has two loops over dim and instantiate it for all pairs like it is
done for
https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/source/dofs/dof_tools_constraints.inst.in#L123

Do you want to try if you can make this work? I am not sure if the
implementation will work correctly, but we can check that in a second
step.
If this works, we would appreciate a pull request from you to fix this!

Best,
Timo



On Sun, Nov 22, 2020 at 4:11 PM Malhar T.  wrote:
>
> Hello All !
>
> I was trying the step-3 with periodic boundary conditions and dim = 2 and 
> spacedim = 3. The main reason behind this is I want to test periodic boundary 
> conditions for dim != spacedim. I have tested the code for dim = spacedim = 2 
> and the periodic boundary conditions and it works (never mind the results for 
> now !). I have also tested the code for dim = 2 and spacedim = 3 without 
> periodic boundary conditions and it works too. But only when I add periodic 
> conditions for that case, the issue arises.
>
> I got the following error,
>
> error: undefined reference to 'void 
> dealii::DoFTools::make_periodicity_constraints<2, 3, 
> double>(dealii::DoFHandler<2, 3> const&, unsigned int, unsigned int, unsigned 
> int, dealii::AffineConstraints&, dealii::ComponentMask const&, 
> double)'
> collect2: error: ld returned 1 exit status
> make[2]: *** [CMakeFiles/step1.dir/build.make:139: step1] Error 1
> make[1]: *** [CMakeFiles/Makefile2:76: CMakeFiles/step1.dir/all] Error 2
> make: *** [Makefile:84: all] Error 2
>
> The version that I am using is 9.3.0-pre. Although I am not really great at 
> C++ but according to the template and source, it should have worked, but I 
> could not find anything here. I tried running it without template arguments 
> or with only '2' as a template argument, but couldn't make it work. I have 
> attached the code which I have tried, and also the code dim = 2 and spacedim 
> = 3 without periodic boundary conditions. Any help will be greatly 
> appreciated.
>
> Thank You
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see 
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups 
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an 
> email to dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/e11a9f41-4750-459d-bedc-ef0b3ae83dddn%40googlegroups.com.



-- 
Timo Heister
http://www.math.clemson.edu/~heister/

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

[deal.II] Periodic Boundary conditions on Step-3 tutorial with dim != spacedim

[Malhar T.]

Hello All !

I was trying the step-3 
 with periodic 
boundary conditions and dim = 2 and spacedim = 3. The main reason behind 
this is I want to test periodic boundary conditions for dim != spacedim. I 
have tested the code for dim = spacedim = 2 and the periodic boundary 
conditions and it works (never mind the results for now !). I have also 
tested the code for dim = 2 and spacedim = 3 without periodic boundary 
conditions and it works too. But only when I add periodic conditions for 
that case, the issue arises. 

I got the following error,

error: undefined reference to 'void 
dealii::DoFTools::make_periodicity_constraints<2, 3, 
double>(dealii::DoFHandler<2, 3> const&, unsigned int, unsigned int, unsigned 
int, dealii::AffineConstraints&, dealii::ComponentMask const&, double)'
collect2: error: ld returned 1 exit status
make[2]: *** [CMakeFiles/step1.dir/build.make:139: step1] Error 1
make[1]: *** [CMakeFiles/Makefile2:76: CMakeFiles/step1.dir/all] Error 2
make: *** [Makefile:84: all] Error 2

The version that I am using is 9.3.0-pre. Although I am not really great at 
C++ but according to the template 

 
and source 
,
 
it should have worked, but I could not find anything here 
.
 
I tried running it without template arguments or with only '2' as a 
template argument, but couldn't make it work. I have attached the code 
which I have tried, and also the code dim = 2 and spacedim = 3 without 
periodic boundary conditions. Any help will be greatly appreciated.

Thank You

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/e11a9f41-4750-459d-bedc-ef0b3ae83dddn%40googlegroups.com.
/* -
 *
 * Copyright (C) 1999 - 2020 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE.md at
 * the top level directory of deal.II.
 *
 * -
 *
 * Authors: Wolfgang Bangerth, 1999,
 *  Guido Kanschat, 2011
 */
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
using namespace dealii;
class Step3
{
public:
  Step3();
  void run();
private:
  void make_grid();
  void setup_system();
  void assemble_system();
  void solve();
  void output_results() const;
  Triangulation<2,3> triangulation;
  FE_Q<2,3>  fe;
  DoFHandler<2,3>dof_handler;
  MappingQ<2,3>  mapping;
  SparsityPattern  sparsity_pattern;
  SparseMatrix system_matrix;
  Vector solution;
  Vector system_rhs;
};
Step3::Step3()
  : fe(1)
  , dof_handler(triangulation)
  , mapping(2)
{}
void Step3::make_grid()
{
  GridGenerator::subdivided_hyper_cube(triangulation, 10, -10.0, 10.0, 1);
  // triangulation.refine_global(5);
  std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
void Step3::setup_system()
{
  dof_handler.distribute_dofs(fe);
  std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
  DynamicSparsityPattern dsp(dof_handler.n_dofs());
  DoFTools::make_sparsity_pattern(dof_handler, dsp);
  sparsity_pattern.copy_from(dsp);
  system_matrix.reinit(sparsity_pattern);
  solution.reinit(dof_handler.n_dofs());
  system_rhs.reinit(dof_handler.n_dofs());
}
void Step3::assemble_system()
{
  QGauss<2> quadrature_formula(fe.degree + 1);
  FEValues<2,3> fe_values(mapping, fe,
quadrature_formula,
update_values | update_gradients | update_JxW_values);
  const unsigned int dofs_per_cell = fe.dofs_per_cell;
  FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell);
  Vector cell_rhs(dofs_per_cell);
  std::vector 

[deal.II] Re: Creating a SOA with nice OOP like interface

['peterrum' via deal.II User Group]

Hi Zachary,

Regarding your two TODOs. I am explaining how we do it in deal.II and then 
you can decide if this approach also fits your needs.

*Index handling*
The simple case is that the indices (level, index within the level) are 
passed to `TriaAccessorBase` in the constructor (like here 
https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/include/deal.II/grid/tria_accessor.h#L341-L344).
 
The more common approach is that `TriaAccessorBase` is used in the context 
of an iterator; the iterator calls the `TriaAccessorBase::operator++()` 
function, which increments the indices internally appropriately (see: 
https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/include/deal.II/grid/tria_accessor.templates.h#L198-L236).

*Allocation*
In deal.II we have dedicated allocation functions, called `reserve_space` 
see:  
https://github.com/dealii/dealii/blob/a56c8585863b774de9dce5f7dd5c334ef56f3a51/source/grid/tria.cc#L2638-L2744.
 
But doing the reservation in the constructor is also a way if you now the 
size at time of construction, which we don't know in the Triangulation 
class, and it does not change over time (e.g., due to adaptivity).

A general remark: it might be also useful - depending on the regularity of 
data-access pattern -  to have a look into explicit SIMD vectorization, 
i.e., in your case to pack 8 orbitals in the case of AVX512 together. For 
this purpose, deal.II has the VectorizedArray struct, which is a wrapper 
class around intrinsic instruction so that one does not have mess with 
these.

I hope this helps. 

Peter

On Sunday, 22 November 2020 at 01:39:40 UTC+1 zachary...@gmail.com wrote:

> Hi Peter,
>
> Thank you for your email.  Here is what I have so far:
>
> ```c++
> // This will be filled in by looping over number of orbitals, with each 
> orbital containing an array of coefs.
> struct Orbital_Data{
> \\TODO: need a constructor that allocates memory for all the orbital's 
> info here?
> private:
>   std::shared_ptr[ ] > orbital_coefs;
>   std::shared_ptr orbital_coefs_index;
> }
>
> struct Orbital{
> \\TODO: How to handle the ith_orbital index?  
> \\ Something like: Orbital(const uint32_t _ith_orbital)
>   inline std:complex
>   coef(uint32_t i) const {
> uint32_t coef_index =
>   orbital_data->orbital_coefs_index[ith_orbital*size_of_coefs_array + 
> i]; // Note the ith_orbital here used as a stride
>   return orbital_data->orbital_coefs[ coefs_index];
>   }
> private:
>   Orbital_Data *orbital_data;
>   uint32_t ith_orbital;  // This is the heart of SOA
> }
> ```
> I would rather allocate enough space instead of using std::vectors in a 
> constructor possible.  Does this look reasonable to you?
>
> My two questions are my TODOs.  I'm thinking about passing in the size the 
> arrays need to be into a constructor for Orbital_Data.  I am looking to 
> fill in this data inside a loop over orbitals.  So I need to know the 
> ith_orbital and pass it into the Orbital struct somehow.
>
> This follows Prof. Wolfgang's example in his lecture.  I am just missing a 
> few little things conceptually.
>
> Cheers,
>
> Zachary
> On Saturday, November 21, 2020 at 2:21:29 PM UTC-6 peterrum wrote:
>
>> Hi Zachary,
>>
>> I haven't watched the lecture, but I think I can, nevertheless, answer 
>> this question since I have been refactoring and generalizing the relevant 
>> data structures in the last months. 
>>
>> The key aspect is that the data is separated from the class that gives 
>> coordinated access to the data. In our case, the `Triangulation` classes 
>> (and similarity the `DoFHandler`) return iterators, which are more or less 
>> wrapper around `TriaAccessor` or `CellAccessor` (which is also a 
>> `TriaAccessor` with some additional functionalities). The core of these two 
>> classes are the internal fields `present_level` and `present_index`, which 
>> function as pointers into some arrays (within the `Triangulation` class).
>>
>> Let's for example have a look at the implementation of 
>> `CellAccessor::subdomain_id()`:
>> ```cpp
>> template 
>> inline types::subdomain_id
>> CellAccessor::subdomain_id() const
>> {
>>   return 
>> this->tria->levels[this->present_level]->subdomain_ids[this->present_index];
>> }
>> ```
>> which access via the `Triangulation` class 
>> `internal::TriangulationImplementation::TriaLevel::subdomain_ids` array. 
>> The functions `child()` and `vertex()` are setup with a similar logic but 
>> are a bit more complicated: the first one has a stride of 2^dim and the 
>> second one uses many indirections (e.g. in the case of 3D cell, querying 
>> the vertex of a face, which queries the vertex of a line - however each 
>> query is similar to the one of `subdomain_id()`).
>>
>> The assumption of this storage is that in one place of your code you loop 
>> over all cells and only need one or a few properties of a cell (e.g. 
>> `subdomain_id` or `boundary_id`). As a consequence, how to store data is 
>> very