Re: [deal.II] make run error: "dyld: Library not loaded"

[Jean-Paul Pelteret]

Dear Zhibo,

I guess that this could be fixed in Candi (? ping Uwe, Timo). As of the El 
Capitan release (I think), its not possible to use DYLD_LIBRARY_PATH anymore if 
SIP (System Integrity Protection) is enabled. For an immediate solution I would 
suggest trying to use install_name_tool to add the full path of (or an @rpath 
to) these problematic libraries. Here’s a link to a StackOverflow post 
 that has more details. I hope that using 
this tool to direct to the location of the linked in library will resolve the 
issue.

Best,
Jean-Paul

> On 22 Aug 2018, at 21:30, Zhibo Luo  wrote:
> 
> 
> Hi all,
> 
> I recently update my MacOs system to latest macOS High Sierra (10.13.6) and 
> crash my dealii. so I decide to reinstall dealii through candi script. The 
> installation works fine except some warnings. The problem is that when I cd 
> to step-1, cmake and make command get passed, however, when I make run, error 
> comes:
> 
> [ 66%] Built target step-1
> [100%] Run step-1 with Debug configuration
> dyld: Library not loaded: libsuperlu_dist.5.dylib
>   Referenced from: 
> /Users/zhiboluo/deal.ii-candi/deal.II-v9.0.0/examples/step-1/./step-1
>   Reason: image not found
> make[3]: *** [CMakeFiles/run] Abort trap: 6
> make[2]: *** [CMakeFiles/run.dir/all] Error 2
> make[1]: *** [CMakeFiles/run.dir/rule] Error 2
> 
> I found similar issue here: 
> https://groups.google.com/forum/#!searchin/dealii/dyld$3A$20Library$20not$20loaded$20%7Csort:date/dealii/fIXfMnaFdEc/zgHjkY3_CAAJ
>  
> 
> 
> And I tried the suggested method by export 
> DYLD_LIBRARY_PATH=$DYLD_LIBRARY_PATH:~/deal.ii-candi/superlu_dist_5.1.2/lib, 
> but it doesn't work for me. 
> 
> I used otool command to output the link libraries:
> 
> zhiboluo$ otool -L step-1
> step-1:
>   
> /Users/zhiboluo/deal.ii-candi/deal.II-v9.0.0/lib/libdeal_II.g.9.0.0.dylib 
> (compatibility version 9.0.0, current version 9.0.0)
>   /Users/zhiboluo/deal.ii-candi/p4est-2.0/DEBUG/lib/libp4est-2.0.dylib 
> (compatibility version 0.0.0, current version 0.0.0)
>   /Users/zhiboluo/deal.ii-candi/p4est-2.0/DEBUG/lib/libsc-2.0.dylib 
> (compatibility version 0.0.0, current version 0.0.0)
>   /usr/lib/libz.1.dylib (compatibility version 1.0.0, current version 
> 1.2.11)
>   @rpath/libmuelu-adapters.12.dylib (compatibility version 12.0.0, 
> current version 12.10.1)
>   @rpath/libmuelu-interface.12.dylib (compatibility version 12.0.0, 
> current version 12.10.1)
>   @rpath/libmuelu.12.dylib (compatibility version 12.0.0, current version 
> 12.10.1)
>   @rpath/libteko.12.dylib (compatibility version 12.0.0, current version 
> 12.10.1)
>   @rpath/libstratimikos.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libstratimikosbelos.12.dylib (compatibility version 12.0.0, 
> current version 12.10.1)
>   @rpath/libstratimikosaztecoo.12.dylib (compatibility version 12.0.0, 
> current version 12.10.1)
>   @rpath/libstratimikosamesos.12.dylib (compatibility version 12.0.0, 
> current version 12.10.1)
>   @rpath/libstratimikosml.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libstratimikosifpack.12.dylib (compatibility version 12.0.0, 
> current version 12.10.1)
>   @rpath/libifpack2-adapters.12.dylib (compatibility version 12.0.0, 
> current version 12.10.1)
>   @rpath/libifpack2.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libanasazitpetra.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libModeLaplace.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libanasaziepetra.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libanasazi.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libamesos2.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libbelostpetra.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libbelosepetra.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libbelos.12.dylib (compatibility version 12.0.0, current version 
> 12.10.1)
>   @rpath/libml.12.dylib (compatibility version 12.0.0, current version 
> 12.10.1)
>   @rpath/libifpack.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libzoltan2.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libpamgen_extras.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libpamgen.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libamesos.12.dylib (compatibility version 12.0.0, current 
> version 12.10.1)
>   @rpath/libgaleri-xpetra.12.dylib (

Re: [deal.II] 3D hp dof handler with distributed triangulation generates error

[Xiukun Hu]

Thank you so much for your explanation, professor. 

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

[Jie Cheng]

Dear Wolfgang

Thanks for your attention. The physical background is that I want to 
interpolate fluid velocity (v^f) onto an arbitrary point (which is in fact 
a solid support point). Instead of doing a standard finite element 
interpolation, I want to do a "delta function interpolation":

v^s = \int_\Omega v^f \delta(x - x^s) d\Omega

where \delta is the "delta function" defined in the previous post. 

I've tried to read your question a couple of times, but I can't seem to 
> figure 
> out what exactly you want to do. I understand the definition of the delta 
> function located around a given point 'p', but what is (in mathematical 
> formulas) what you want to do in the interpolation process? You have an 
> input 
> and an output vector, but I don't think I understand how they are related. 
>

The logic of the Interpolator class is: for a specific target point, 
iterate over all the nodes on the source grid, compute the distance between 
the node and the target point, then the weight. Next, given a source 
vector, the Interpolator does the interpolation using the 
previously-computed weights.

1. However, there is no convenient way to iterate over all the "nodes" on a 
grid and query their corresponding dofs. Therefore I have to iterate over 
all cells, then iterate over all the support points and find which support 
points are contributing to the target point.

2. In addition, I may want to interpolate more than one vector onto a 
specific point, so I want to cache the influential nodes, compute their 
weights, and reuse them.

3. The tuple is used to mark the "node" that is influential to the target 
point. Intuitively one would use a "global node id" for that, but deal.II 
doesn't work this way. So I use {cell_pointer, support_point_id,  weight}, 
which essentially identifies a "node". While doing interpolation, I must be 
careful that I don't use the same node for multiple times.

Yes, system_to_component_index() returns a std::pair, that was a bug.

What I did is as following:

  template 
  class DiracDeltaInterpolator
  {
  public:
DiracDeltaInterpolator(const DoFHandler &, const Point &, 
double);
void interpolate(const VectorType &,
 Vector &);

  private:
/// Source DoFHandler
const DoFHandler &dof_handler;
/// Target point to interpolate to
const Point ⌖
/// Background mesh size
double h;
/**
 * Source nodes that contribute to the target point,
 * denoted as tuples of cell iterator, supporting point id,
 * and weight
 * as there is no convenient way to express "global node id".
 */
std::vector::active_cell_iterator,
   unsigned int,
   double>>
  sources;
  };

  template 
  DiracDeltaInterpolator::DiracDeltaInterpolator(
const DoFHandler &dof_handler, const Point &point, double h)
: dof_handler(dof_handler), target(point), h(h)
  {
const FiniteElement &fe = dof_handler.get_fe();
const std::vector> &unit_points = 
fe.get_unit_support_points();
Quadrature dummy_q(unit_points.size());
MappingQGeneric mapping(1);
FEValues dummy_fe_values(
  mapping, fe, dummy_q, update_quadrature_points);
std::vector dof_indices(fe.dofs_per_cell);
for (auto cell : dof_handler.active_cell_iterators())
  {
dummy_fe_values.reinit(cell);
cell->get_dof_indices(dof_indices);
// Real coordinates of the current cell's support points
auto support_points = dummy_fe_values.get_quadrature_points();
for (unsigned int v = 0; v < unit_points.size(); ++v)
  {
// Compute the weight of each support point
double weight = 1;
for (unsigned int d = 0; d < dim; ++d)
  {
double xbar = (unit_points[v][d] - target[d]) / h;
weight *= (std::abs(xbar) <= 2
 ? xbar * 0.25 * (1 + std::cos(M_PI * xbar / 2))
 : 0);
  }
// If weight is nonzero then it is an influential point
if (weight > 0)
  {
sources.push_back(std::make_tuple(cell, v, weight));
  }
  }
  }
  }

  template 
  void DiracDeltaInterpolator::interpolate(
const VectorType &source_vector,
Vector &target_vector)
  {
std::vector bs(dof_handler.n_dofs(), false);
target_vector = 0;
const FiniteElement &fe = dof_handler.get_fe();
std::vector dof_indices(fe.dofs_per_cell);
for (auto p : sources)
  {
std::get<0>(p)->get_dof_indices(dof_indices);
// Vector component of this support point
auto d = fe.system_to_component_index(std::get<1>(p)).first;
AssertThrow(d < dim, ExcMessage("Vector component not less than 
dim!"));
// Global dof index of this support point
auto index = dof_indices[std::get<1>(p)];
AssertThrow(index < dof_handler.n_dofs(),
Ex

Re: [deal.II] 3D hp dof handler with distributed triangulation generates error

[Wolfgang Bangerth]

On 08/23/2018 05:51 PM, Xiukun Hu wrote:


I've been trying to parallelize step-46 for a while (mainly by combining with 
step-40). Now I can make the code run in release mode and get the correct 
output, but still get (different kind of) errors in debug mode. After some 
efforts I finally reproduced one in a quite short code.


There will likely be many more errors to come, and I would doubt that your 
results will be correct even in release mode. The issue, fundamentally, is 
that in deal.II 9.0, hp::DoFHandler does not yet work with 
parallel::distributed::Triangulation. With some luck, this combination will be 
available in deal.II 9.1, though.


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.

[deal.II] 3D hp dof handler with distributed triangulation generates error

[Xiukun Hu]

Hi everybody,

I've been trying to parallelize step-46 for a while (mainly by combining 
with step-40). Now I can make the code run in release mode and get the 
correct output, but still get (different kind of) errors in debug mode. 
After some efforts I finally reproduced one in a quite short code.

*Code To Reproduce the Error:*

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 


int main(int argc, char *argv[]) {
 using namespace dealii;

 Utilities::MPI::MPI_InitFinalize mpi_init_fin(argc, argv, 1);
 MPI_Comm mpi_comm(MPI_COMM_WORLD);

 constexpr int dim = 3;

 parallel::distributed::Triangulation mesh(mpi_comm);
 // Triangulation mesh;
 GridGenerator::subdivided_hyper_cube(mesh, 8, -1, 1);
 hp::DoFHandler dof_handler(mesh);

 hp::FECollection fe_collection(FE_Q(2), FE_Q(2));
 for (auto && cell:dof_handler.active_cell_iterators()){
   if (cell->is_locally_owned()) {
 cell->set_active_fe_index(0);
   }
 }

 dof_handler.distribute_dofs(fe_collection);

 std::cout << "Done.\n";

 return 0;
}


*Error Message:*


An error occurred in line <3654> of file 

 
in function
unsigned int dealii::DoFCellAccessor, 
false>::active_fe_index() const [DoFHandlerType = dealii::hp::DoFHandler<3, 
3>, lda = false]
The violated condition was:
(dynamic_cast*> 
(this->dof_handler) != nullptr) || (this->is_locally_owned() || 
this->is_ghost())
Additional information:
You can only query active_fe_index information on cells that are either 
locally owned or (after distributing degrees of freedom) are ghost cells.

Stacktrace:
---
#0  2   libdeal_II.g.9.0.0.dylib0x00011842d31c 
_ZNK6dealii15DoFCellAccessorINS_2hp10DoFHandlerILi3ELi3EEELb0EE15active_fe_indexEv
 
+ 380: 2   libdeal_II.g.9.0.0.dylib0x00011842d31c 
_ZNK6dealii15DoFCellAccessorINS_2hp10DoFHandlerILi3ELi3EEELb0EE15active_fe_indexEv
#1  3   libdeal_II.g.9.0.0.dylib0x0001195af850 
_ZN6dealii8internal2hp24DoFHandlerImplementation14Implementation13reserve_spaceILi3EEEvRNS_2hp10DoFHandlerILi3EXT_EEE
 
+ 304: 3   libdeal_II.g.9.0.0.dylib0x0001195af850 
_ZN6dealii8internal2hp24DoFHandlerImplementation14Implementation13reserve_spaceILi3EEEvRNS_2hp10DoFHandlerILi3EXT_EEE
#2  4   libdeal_II.g.9.0.0.dylib0x0001195aeb06 
_ZN6dealii2hp10DoFHandlerILi3ELi3EE15distribute_dofsERKNS0_12FECollectionILi3ELi3EEE
 
+ 630: 4   libdeal_II.g.9.0.0.dylib0x0001195aeb06 
_ZN6dealii2hp10DoFHandlerILi3ELi3EE15distribute_dofsERKNS0_12FECollectionILi3ELi3EEE
#3  5   rep 0x00010bb5e2fd main + 413: 
5   rep 0x00010bb5e2fd main
#4  6   libdyld.dylib   0x7fff6004f015 start + 1: 6 
  libdyld.dylib   0x7fff6004f015 start


Calling MPI_Abort now.
...

*Environment:*

macOS High Sierra + llvm 9.1.0 + cmake 3.12.1 + deal.II 9.0.0

*Observations:*

Won't throw error if in any following situation:

   - build type is release,
   - dim = 2,
   - use (non-distributed) Triangulation class, or
   - run with only one processor.

The error happens inside distribute_dofs(). I'm not sure if I used the 
wrong classes/functions for parallel or there is a bug in dealii. I'd be 
very glad to see some examples/tutorials of using hp::* in distributed 
system, but couldn't find any yet.

Thanks ahead for any help.

Xiukun

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

[Wolfgang Bangerth]

Jie,
I've tried to read your question a couple of times, but I can't seem to figure 
out what exactly you want to do. I understand the definition of the delta 
function located around a given point 'p', but what is (in mathematical 
formulas) what you want to do in the interpolation process? You have an input 
and an output vector, but I don't think I understand how they are related.



I changed the type of DiracDeltaInterpolator::sources to 
std::vector::active_cell_iterator,

    
   unsigned int,
  
 double>>
which stores the cell iterator, id of support point, weight of an "influential 
node". So that DiracDeltaInterpolator::interpolate can be implemented as:


   template 
   void DiracDeltaInterpolator::interpolate(
   const VectorType &source_vector,
   Vector &target_vector)
   {
     target_vector = 0;
     const FiniteElement &fe = dof_handler.get_fe();
     std::vector dof_indices(fe.dofs_per_cell);
     for (auto p : sources)
   {
     std::get<0>(p)->get_dof_indices(dof_indices);
     // Vector component of this support point
     auto d = fe.system_to_component_index(std::get<1>(p));
     Assert(d < dim, ExcMessage("Vector component not less than dim!"));
     // Global dof index of this support point
     auto index = dof_indices[std::get<1>(p)];
     // Interpolate
     target_vector[d] += std::get<2>(p) * source_vector[index];


Can you explain what your tuple now contains? In your first implementation, it 
was
  * cell
  * index of a dof within a cell
But it's not clear to me whether that is still true here. Also, 
system_to_component_index() returns a std::pair, so that's going to be the 
data type of 'd'. But then you can write 'd'target_vector[d]' in the last line.


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] error while running step-40 on crayXC40 with petsc and mpi

[Wolfgang Bangerth]

Finally I get the following error while running step-40 with 24 
processes(single node) or 48 processes(two nodes).

|

Anerror occurred inline <250>of file 
infunction

voiddealii::PETScWrappers::MPI::Vector::create_vector(unsignedint,unsignedint)
Theviolated condition was:
     size()==n
Additionalinformation:
Dimension4294967295notequal to 0.

Stacktrace:
---
#0  /mnt/lustre/serc/secbteja/bhanu3/dealiiinst/lib/libdeal_II.g.so.9.1.0-pre: 
dealii::PETScWrappers::MPI::Vector::create_vector(unsigned int, unsigned int)
#1  /mnt/lustre/serc/secbteja/bhanu3/dealiiinst/lib/libdeal_II.g.so.9.1.0-pre: 
dealii::PETScWrappers::MPI::Vector::Vector()

#2  ./step-40: Step40::LaplaceProblem<2>::LaplaceProblem()
#3  ./step-40: main

|
Incidentally 4294967295 = 2^32 - 1.
Please suggest.


Bhanu,
it's hard to tell whether this is specific to Cray, is a compiler bug, a bug 
in deal.II, or anything else. Does it work with 1 processor? With 2 
processors? With 4? Etc? The program does run on my machine with 24 
processors, so the question is what exactly is different on your system.


By the way: The error happens in this function
```
void
Vector::create_vector(const size_type n, const size_type local_size)
{
  (void)n;
  Assert(local_size <= n, ExcIndexRange(local_size, 0, n));
  ghosted = false;

  const PetscErrorCode ierr =
VecCreateMPI(communicator, local_size, PETSC_DETERMINE, &vector);
  AssertThrow(ierr == 0, ExcPETScError(ierr));

  Assert(size() == n, ExcDimensionMismatch(size(), n));// line 250
}
```
So the error suggests that after the function creates the vector, size() of 
the vector returns (unsigned int)-1. It's not clear to me how this can happen. 
If you run this in a debugger, can you check what the arguments 'n' and 
'local_size' to this function are?


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: How do I choose my spaces to enforce gradient continuity across an interface?

[Wolfgang Bangerth]

Thank you for your reply. My interface condition is \integral_{Interface} 
p*I*n_stokes - jump(velocity_gradient) ds,


Do you really mean that the *integral* is zero, or should this hold pointwise? 
Either way, I am pretty sure that you can't expect this to be true pointwise 
for the discrete solution. Furthermore, what it really says is *not* that the 
velocity gradient must be continuous, but that the gradient has a jump equal 
to the pressure times the normal vector. So a method that guarantees a C^1 
solution would not actually help you.



where I is the interface and 
jump(velocity_gradient) = vel_grad_stokes*n_stokes + 
vel_grad_laplace*n_laplace.  I am using this with zero boundary conditions on 
pressure.  The interface  in the code I attached is y=0.5, with y>=0.5 stokes 
and below y=0.5 laplace.  When I let the mesh size go to zero, even when I 
prescribe zero pressure everywhere (by defining zero contribution from 
pressure on my right hand side) I am still not getting zero pressure as the 
mesh size goes to zero.  Furthermore, as my mesh size goes to zero my velocity 
converges to something that is not related to the right-hand side I chose.  
Specifically, it feautures a non zero jump in the gradient across the 
interface, even though I am prescribing the same right-hand side to both sides 
(i.e. with zero pressure on both sides).  This does not happen if I solve just 
stokes everywhere or just laplace (with the same right hand side :i.e. zero 
contribution from pressure).


I don't have a suggestion here other than simplifying the problem and spending 
some time gathering other cues. Does the code work if one or the other medium 
makes up the whole domain?


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.

[deal.II] Re: error while running step-40 on crayXC40 with petsc and mpi

[Bhanu Teja]

I checked with other versions of MPICH, petsc and tpsl also, but the same 
error occurs or the installation fails.

Thanks,
Bhanu.


>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group 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 to deal.ii's matrix structures after assembly for solution using an external solver

[Bruno Turcksin]

Pratik,

Le jeu. 23 août 2018 à 09:34,  a écrit :

>
> Thank you. So, I can get the matrix values by &system_matrix.begin() and
> the column indices by doing &sparsity_pattern.begin() but I am a bit
> confused as to how I get the row pointers. Or did I misunderstand something
> ?
>
Take a look here

for the method I have explained. Note that you do not get access to a
pointer to the rows and columns but to an iterator to a structure which you
can ask for the row and column indices. If you need to have access to the
raw pointers, you need to modify deal.II by making your matrix class a
friend of SparseMatrix (see here

)

Best,

Bruno

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

[deal.II] Re: Access to deal.ii's matrix structures after assembly for solution using an external solver

[pratik . hpc]

Bruno,

Thank you. So, I can get the matrix values by &system_matrix.begin() and 
the column indices by doing &sparsity_pattern.begin() but I am a bit 
confused as to how I get the row pointers. Or did I misunderstand something 
?

Thanks,
Pratik.

On Thursday, August 23, 2018 at 3:05:07 PM UTC+2, Bruno Turcksin wrote:
>
> Pratik,
>
> On Thursday, August 23, 2018 at 8:38:34 AM UTC-4, prati...@gmail.com 
>  wrote:
>>
>> Essentially, it says the the member variable cols within system_matrix is 
>> private. Is there any other way or a remedy to this ? How do I access this 
>> matrix and rhs data ?
>>
> You can get the underlying data of the matrix using this function 
> .
>  
> To get the column indices, you first need to get the SparsityPattern of the 
> matrix, then you can get the underlying data using this 
> .
>  
> This will give an iterator 
> 
>  
> that you can dereference to an accessor. 
> 
>  
> This gives you everything that you need. Another solution is to modify 
> deal.II directly. You could just make the variables public (the easiest but 
> also the dirtiest method) or make your own Matrix class a friend of 
> SparseMatrix.
>
> Best,
>
> Bruno
>

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

[deal.II] Re: Access to deal.ii's matrix structures after assembly for solution using an external solver

[Bruno Turcksin]

Pratik,

On Thursday, August 23, 2018 at 8:38:34 AM UTC-4, pratik@gmail.com 
wrote:
>
> Essentially, it says the the member variable cols within system_matrix is 
> private. Is there any other way or a remedy to this ? How do I access this 
> matrix and rhs data ?
>
You can get the underlying data of the matrix using this function 
.
 
To get the column indices, you first need to get the SparsityPattern of the 
matrix, then you can get the underlying data using this 
.
 
This will give an iterator 

 
that you can dereference to an accessor. 

 
This gives you everything that you need. Another solution is to modify 
deal.II directly. You could just make the variables public (the easiest but 
also the dirtiest method) or make your own Matrix class a friend of 
SparseMatrix.

Best,

Bruno

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

Re: [deal.II] Small bug in ConstraintMatrix::add_entries () ?!

[Dustin Kumor]

Daniel,

thanks for the hint. I will have a look at this on the weekend.

Best,
Dustin

Am Donnerstag, 23. August 2018 13:45:41 UTC+2 schrieb Daniel Arndt:
>
> Dustin,
>
> I can put that into a pull request but I've never done this before so I'm 
>> going to need someone guiding me through the whole process.
>> Same for the test writing issue.
>>
> Have a look at https://github.com/dealii/dealii/wiki/Contributing and 
> don't hesitate to ask if you have questions.
>
> Best,
> Daniel
>

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

Re: [deal.II] Re: How do I choose my spaces to enforce gradient continuity across an interface?

[georgios . sialounas]

Dear Prof. Bangerth,

Thank you for your reply. My interface condition is \integral_{Interface} 
p*I*n_stokes - jump(velocity_gradient) ds, where I is the interface and 
jump(velocity_gradient) = vel_grad_stokes*n_stokes + 
vel_grad_laplace*n_laplace.  I am using this with zero boundary conditions 
on pressure.  The interface  in the code I attached is y=0.5, with y>=0.5 
stokes and below y=0.5 laplace.  When I let the mesh size go to zero, even 
when I prescribe zero pressure everywhere (by defining zero contribution 
from pressure on my right hand side) I am still not getting zero pressure 
as the mesh size goes to zero.  Furthermore, as my mesh size goes to zero 
my velocity converges to something that is not related to the right-hand 
side I chose.  Specifically, it feautures a non zero jump in the gradient 
across the interface, even though I am prescribing the same right-hand side 
to both sides (i.e. with zero pressure on both sides).  This does not 
happen if I solve just stokes everywhere or just laplace (with the same 
right hand side :i.e. zero contribution from pressure).

Kind regards,
Georgios Sialounas

On Wednesday, August 22, 2018 at 5:42:34 AM UTC+1, Wolfgang Bangerth wrote:
>
>
> Georgios, 
>
> > Thanks for your reply! Essentially what I was trying to do is patch 
> together 
> > two product spaces:  ([H1_0]^2, L^2_0) in one part of the domain and 
> > ([H1_0]^2, L^2_0) in another, with an interface condition to ensure 
> > solvability for my problem (Stokes in one part of the domain, laplace in 
> > another).  The interface condition was chosen such that there would be 
> no jump 
> > of the velocity gradient across the interface (this was done by choosing 
> an 
> > L2_0 pressure that goes in the right-hand side). 
>
> Can you state this interface condition in mathematical terms? In general, 
> functions in H^1 can be expected to be continuous along interfaces, but 
> their 
> gradients can not. Indeed, their normal derivatives along the interface 
> will 
> only be in H^{-1/2}, and so the difference between the normal derivatives 
> from 
> both sides of the interface can at best be expected to be weakly 
> continuous 
> (i.e., when tested with a sufficiently smooth test function). 
>
>
> >   Unfortunately this did not 
> > work in the implementation phase and I ended up getting a jump of the 
> velocity 
> > gradient across the interface. 
>
> But is this different than the jump in the gradient between each cell and 
> the 
> next? There, the gradient is not continuous either. The question is what 
> happens as you let the mesh size go to zero. 
>
>
> >  At the moment I don't know if this is a coding 
> > bug, or an erroneous choice of spaces.  To be honest I expected weak 
> > differentiability to be enough to ensure a zero velocity gradient jump 
> across 
> > the interface too. As such, I was thinking how I should choose my spaces 
> to 
> > ensure gradient continuity across the interface to see if this would fix 
> my 
> > problem but I do not know how to enforce that in deal.ii.  As such, any 
> advice 
> > on how to do that would be great. 
>
> It's actually quite complicated to ensure continuity of the gradient 
> strongly. 
> This is why the biharmonic equation is so difficult to solve! 
>
> Best 
>   Wolfgang 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

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

[deal.II] Access to deal.ii's matrix structures after assembly for solution using an external solver

[pratik . hpc]

Hello,

I am looking to use deal.ii to generate an assembled sparse matrix from a 
physical problem, but would like to use my own solver software. For 
example, if you consider the step-9, 

In AdvectionProblem::solve() function I would get the system-matrix, 
system_rhs from deal.ii's data structures and use them to get the solution 
vector and convert them back to deal.ii's type for 
hanging_node_constraints.distribute(solution) function.

Therefore, I would like to replace


 SolverControl   solver_control (1000, 1e-12);
 SolverBicgstab<>bicgstab (solver_control);

 PreconditionJacobi<> preconditioner;
 preconditioner.initialize(system_matrix, 1.0);

 bicgstab.solve (system_matrix, solution, system_rhs,
 preconditioner);

with my own solver. The purpose is to test our solver and be able to 
interface it with deal.ii. As I see it, I need:

1. the system_matrix to be converted to CSR format. I do this by using the 
member variables of the system_matrix, namely cols, which again has 
rowstart, colnums, so I should be able to access the column indices and row 
pointers by system_matrix.cols->colnums.get() and 
system_matrix.cols->rowstart.get() and the values by 
system_matrix.val.get(). The size of the matrix by system_matrix.m() and 
system_matrix.n() and the rhs values by system_rhs.values.get(). But I run 
into problems when I do this. I get the following error:

error: ‘dealii::SmartPointer > dealii::SparseMatrix::cols’ is 
private within 
this context
 auto A = mtx::create(exec,gko::dim<2>(system_matrix.cols->cols),
   ~~^~~~
/home/pratik/Documents/10Softwares/dealii/install_dealii/include/deal.II/lac/sparse_matrix.h:1615:61:
 
note: declared private here
   SmartPointer > cols;

Essentially, it says the the member variable cols within system_matrix is 
private. Is there any other way or a remedy to this ? How do I access this 
matrix and rhs data ?

I am sorry for the long post. I wanted to give as much information as 
possible.

Thank you,
Pratik.

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group 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] Small bug in ConstraintMatrix::add_entries () ?!

[Daniel Arndt]

Dustin,

I can put that into a pull request but I've never done this before so I'm 
> going to need someone guiding me through the whole process.
> Same for the test writing issue.
>
Have a look at https://github.com/dealii/dealii/wiki/Contributing and don't 
hesitate to ask if you have questions.

Best,
Daniel

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

Re: [deal.II] Small bug in ConstraintMatrix::add_entries () ?!

[Dustin Kumor]

Dear Wolfgang

I can put that into a pull request but I've never done this before so I'm 
going to need someone guiding me through the whole process.
Same for the test writing issue.

Best regards
Dustin

Am Mittwoch, 22. August 2018 23:25:11 UTC+2 schrieb Wolfgang Bangerth:
>
>
> Dustin, 
> good observation. I think you have found a bug indeed. I've put this 
> here for now: 
>https://github.com/dealii/dealii/issues/7101 
> Your fix seems reasonable to me. Do you want to put that into a pull 
> request? Let us know if you need help with that (and with writing a test). 
>
> (For reference, post-9.0, the class is now called AffineConstraints.) 
>
> Best and thanks for pointing this out! 
>   Wolfgang 
>
>
> On 08/22/2018 08:30 AM, Dustin Kumor wrote: 
> > Dear all, 
> > 
> > in my opinion there is a small bug in the implementation of the 
> > ConstraintMatrix member function add_entries in file 
> > constraint_matrix.cc (release Version 9.0, lines 224-261). The current 
> > piece of code is: 
> > 
> > | 
> > 224void 
> > 225ConstraintMatrix::add_entries 
> > 226(constsize_type  line, 
> > 227conststd::vector>&col_val_pairs) 
> > 228{ 
> > 229Assert(sorted==false,ExcMatrixIsClosed()); 
> > 230Assert(is_constrained(line),ExcLineInexistant(line)); 
> > 231 
> > 232ConstraintLine*line_ptr 
> =&lines[lines_cache[calculate_line_index(line)]]; 
> > 233Assert(line_ptr->index ==line,ExcInternalError()); 
> > 234 
> > 235// if in debug mode, check whether an entry for this column already 
> > 236// exists and if its the same as the one entered at present 
> > 237// 
> > 238// in any case: skip this entry if an entry for this column already 
> > 239// exists, since we don't want to enter it twice 
> > 240for(std::vector>::const_iterator 
> > 241   col_val_pair =col_val_pairs.begin(); 
> > 242   col_val_pair!=col_val_pairs.end();++col_val_pair) 
> > 243{ 
> > 244Assert(line !=col_val_pair->first, 
> > 245ExcMessage("Can't constrain a degree of freedom to itself")); 
> > 246 
> > 247for(ConstraintLine::Entries::const_iterator 
> > 248   p=line_ptr->entries.begin(); 
> > 249   p !=line_ptr->entries.end();++p) 
> > 250if(p->first ==col_val_pair->first) 
> > 251{ 
> > 252// entry exists, break innermost loop 
> > 253Assert(p->second ==col_val_pair->second, 
> > 254ExcEntryAlreadyExists(line,col_val_pair->first, 
> > 255   
> > p->second,col_val_pair->second)); 
> > 256break; 
> > 257} 
> > 258 
> > 259  line_ptr->entries.push_back (*col_val_pair); 
> > 260} 
> > 261} 
> > | 
> > 
> > Although, it is stated in the comments (lines 235-239) there is no skip 
> > of the entry col_val_pairin the case this entry already exist. The 
> > breakcommand terminates the innermost loop but the push_backoperation is 
> > performed nevertheless, i.e. for every element of col_val_pairs. This 
> > leads to duplicates in the constraint matrix. 
> > I solved the problem by just addind the three lines 246a, 255a, 258a. 
> > 
> > | 
> > 224void 
> > 225ConstraintMatrix::add_entries 
> > 226(constsize_type line, 
> > 227conststd::vector>&col_val_pairs) 
> > 228{ 
> > 229Assert(sorted==false,ExcMatrixIsClosed()); 
> > 230Assert(is_constrained(line),ExcLineInexistant(line)); 
> > 231 
> > 232ConstraintLine*line_ptr 
> =&lines[lines_cache[calculate_line_index(line)]]; 
> > 233Assert(line_ptr->index==line,ExcInternalError()); 
> > 234 
> > 235// if in debug mode, check whether an entry for this column already 
> >   236// exists and if its the same as the one entered at present 
> >   237// 
> >   238// in any case: skip this entry if an entry for this column already 
> >   239// exists, since we don't want to enter it twice 
> >   240for(std::vector >::const_iterator 
> >   241   col_val_pair = col_val_pairs.begin(); 
> >   242   col_val_pair!=col_val_pairs.end(); ++col_val_pair) 
> >   243{ 
> >   244Assert(line != col_val_pair->first, 
> >   245ExcMessage("Can't constrain a degree of freedom to itself")); 
> >   246 
> > 246aboolentry_exists =false; 
> > 247for(ConstraintLine::Entries::const_iterator 
> > 248   p=line_ptr->entries.begin(); 
> > 249   p !=line_ptr->entries.end();++p) 
> > 250if(p->first ==col_val_pair->first) 
> > 251{ 
> > 252// entry exists, break innermost loop 
> >   253Assert(p->second == col_val_pair->second, 
> >   254ExcEntryAlreadyExists(line, col_val_pair->first, 
> >   255  p->second, 
> > col_val_pair->second)); 
> > 255a   entry_exists =true; 
> > 256break; 
> > 257} 
> > 258 
> > 258aif(entry_exists ==fase) 
> > 259  line_ptr->entries.push_back (*col_val_pair); 
> > 260} 
> > 261} 
> > | 
> > 
> > May be there is a more elegant way to do this. 
> > 
> > Best regards 
> > Dustin 
> > 
> > -- 
> > The deal.II project is located at http://www.dealii.org/ 
> > For mailing list/forum options, see 
> > http