[deal.II] Re: Reading in higher order meshes from GMSH

2020-09-15 Thread 'peterrum' via deal.II User Group
Dear Sepehr,

If I understand you correctly you would like to have support for quad9 and 
hex27.

We have an open pull request for this issue (see: 
https://github.com/dealii/dealii/pull/10163). Hopefully we get that PR (or 
some version of it) merged soon.

Regards,
Peter

On Tuesday, 15 September 2020 at 17:26:45 UTC+2 bananacreami...@gmail.com 
wrote:

> Hello,
>
> My name is Sepehr and im an undergraduate mechanical engineering student 
> at McGill.
>
> I am currently looking to see if Dealii has support for reading in higher 
> order meshes.
>
> I am aware we can read in a course linear mesh and us built in functions 
> to add nodes in between to match out desired shape. However, making a 
> higher order mesh in GMSH is more desirable.
>
> The only similar discussion board I could find was from 9 years ago:
> https://dealii.narkive.com/93ybreg6/gmsh-higher-order-elements
>
> I was wondering maybe the new version of dealii has addressed this,
>
> Thank you for your time,
>
> Sepehr
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"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/d773248c-1f5c-4b0c-96fc-1d93a37b7b2bn%40googlegroups.com.


[deal.II] Re: Error while installing

2020-07-01 Thread 'peterrum' via deal.II User Group
The error message already tells you what to do. You need to run the commend 
with admi privileges, i.e., `sudo make install`.

Hope that helps,
Peter

On Wednesday, 1 July 2020 13:21:12 UTC+2, ME20S001 Bardawal Prakash wrote:
>
> Hello,
> Someone help me solving this issue, here I'm attaching snapshots of 
> that error message. And also the installation files was not created after 
> running 'make install' command.
>
>

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


[deal.II] Re: Is it possible to copy_triangulation for fullydistributed with periodic face?

2020-06-11 Thread 'peterrum' via deal.II User Group
Dear Heena,

may I ask you to be more specific regarding to 
parallel::fullydistributed::Triangualation (p:f:t) error. In the case of 
p:f:t you can copy indeed refined meshes, however users need to deal with 
periodicity on their own by applying the periodicy once again. See the 
following test: 
https://github.com/dealii/dealii/blob/master/tests/fullydistributed_grids/copy_serial_tria_04.cc#L102

This is kinda annoying, but I was not able to come up with a more 
transparent solution during the development of p:f:t.

Hope that helps,
Peter

On Thursday, 11 June 2020 10:06:49 UTC+2, heena patel wrote:
>
> Dear all,
>  I am trying to solve semilagrangian advection 
> problem?
> I have traingulation with periodic face.
> I wanted to copy_triangulation and trace back mesh. The thing is with 
> *parallel 
> distributed triangulation* it refuse to copy for refine mesh, that is 
> written in documentation as well. When  I try to use fullydistributed 
> triangulation it gives me error with periodicity. Is there a way that I can 
> copy triangulation with periodicity or I am missing something? The code 
> works if I create new_triangulation with parallel distributed triangulation 
> and work on it. But it is expensive step. Kindly find code in attachment.
>
>
> Regards,
> Heena
>

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


[deal.II] Re: Broadcasting packed objects

2020-06-08 Thread 'peterrum' via deal.II User Group
What you could also do is to turn compression off.

Peter

On Monday, 8 June 2020 14:19:25 UTC+2, peterrum wrote:
>
> Dear Maurice,
>
> The problem is that the size of `auto buffer = 
> dealii::Utilities::pack(r1);` is not the same on all processes, which is a 
> requirement if you use `MPI_Bcast`. My suggestion would to split the 
> procedure into two steps: 1) bcast the size on rank 1; 2) bcast the actual 
> data.
>
> Peter
>
> On Monday, 8 June 2020 12:46:53 UTC+2, Maurice Rohracker wrote:
>>
>> Dear deal.II community,
>>
>> For a distributed implementation, we would like to broadcast packed data 
>> objects.
>>
>> For the beginning, we would like to understand how the serialize function 
>> for packing an object is working in the sense of the deal.II way, before it 
>> is getting more complicated.
>>
>> I created a small test with a class, which should be broadcasted. As a 
>> comparison, I took the dealii:Point, which can obviously be packed and 
>> broadcasted without any issues.
>>
>> Unfortunately, the broadcast for our class does not work. We assume that 
>> our serialize function is not in the correct manner. Is it possible to pack 
>> own objects using the dealii::Utilities:pack() function? What is the proper 
>> way defining the serialize function?
>>
>> Attached you'll find our small example code. We are using deal.II 9.0.1.
>>
>> Thank you very much in advance.
>>
>> Best regards,
>> Maurice Rohracker
>> Master Student Computational Engineering
>> FAU Erlange-Nürnberg
>>
>

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


Re: [deal.II] Step-4 1D problem

2020-06-08 Thread 'peterrum' via deal.II User Group
Indeed! Christoph, you seem to be right!

Feel free to create a pull request on GitHub for this inconsistency! We 
will help you if you need some assistance! Amazing that there are still 
errors in the first tutorials although - probably - all deal.II user have 
had a look at these...

Thanks,
Peter

On Monday, 8 June 2020 17:50:25 UTC+2, Christoph Kammer wrote:
>
> Thank you Praveen,
>
> that solved my problem. After fixing the issue with the boundary id things 
> made sense. Is there a particular reason that in 1D the boundary points are 
> labeled in a different manner than in 2+ D?
>
> PS. deal ii team. In the documentation on step 3 + 4 it states that the 
> Poisson problem is solved on a domain [0,1]^dim. The actual solver however 
> is implemented to solve this problem on the domain [-1,1]^dim. Perhaps the 
> introduction to this tutorial should be made consistent with the commented 
> program.
>
>
> Thanks again,
>
>
> Christoph
>
>
> On 6/7/20 1:14 AM, Praveen C wrote:
>
> In 1-D, 
>
> GridGenerator::hyper_cube(triangulation, -1, 1);
>
> will assign different boundary indicators to left and right side, even 
> though colorize=false is default, see
>
>
> https://dealii.org/developer/doxygen/deal.II/namespaceGridGenerator.html#acea0cbcd68e52ce8113d1134b87de403
>
>  I think left=0 and right=1. Due to this reason, bc may not be applied on 
> both sides.
>
> Do the apply bc twice, once with indicator=0 and once with indicator=1.
>
>   std::map boundary_values;
>   VectorTools::interpolate_boundary_values(dof_handler,
>0, // left boundary
>BoundaryValues(),
>boundary_values);
>   MatrixTools::apply_boundary_values(boundary_values,
>  system_matrix,
>  solution,
>  system_rhs);
>
>   VectorTools::interpolate_boundary_values(dof_handler,
>1, // right boundary
>BoundaryValues(),
>boundary_values);
>   MatrixTools::apply_boundary_values(boundary_values,
>  system_matrix,
>  solution,
>  system_rhs);
>
> Best
> praveen
>
> On 07-Jun-2020, at 9:50 AM, 'Christoph Kammer' via deal.II User Group <
> dea...@googlegroups.com > wrote:
>
> Hi Deal.ii team,
>
> I went through tutorial step-4 and I tried to use the solver for Poisson's 
> problem in 2D and rewrite it for a 1D problem. This is basically step-4 and 
> setting RHS f =1 and solving for homogeneous BCs. The program is compiling 
> just fine and executing, however, the solution I am getting from that 
> doesn't make any sense. The BCs aren't even satisfied at x = 1.
>
> Could anyone let me know what's going on here? Setting f = 1 and changing 
> everything to solve for homogeneous BC in 1D seems like a very simple 
> problem.
>
> Thank you,
>
> Christoph
>
>
> -- 
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see 
> https://groups.google.com/d/forum/dealii?hl=en
> --- 
> You received this message because you are subscribed to the Google Groups 
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an 
> email to dea...@googlegroups.com .
> To view this discussion on the web visit 
> https://groups.google.com/d/msgid/dealii/C66D5394-FFCF-426E-A66F-9F3CA5A0F7B0%40gmail.com
>  
> 
> .
>
>

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


[deal.II] Re: Broadcasting packed objects

2020-06-08 Thread 'peterrum' via deal.II User Group
Dear Maurice,

The problem is that the size of `auto buffer = 
dealii::Utilities::pack(r1);` is not the same on all processes, which is a 
requirement if you use `MPI_Bcast`. My suggestion would to split the 
procedure into two steps: 1) bcast the size on rank 1; 2) bcast the actual 
data.

Peter

On Monday, 8 June 2020 12:46:53 UTC+2, Maurice Rohracker wrote:
>
> Dear deal.II community,
>
> For a distributed implementation, we would like to broadcast packed data 
> objects.
>
> For the beginning, we would like to understand how the serialize function 
> for packing an object is working in the sense of the deal.II way, before it 
> is getting more complicated.
>
> I created a small test with a class, which should be broadcasted. As a 
> comparison, I took the dealii:Point, which can obviously be packed and 
> broadcasted without any issues.
>
> Unfortunately, the broadcast for our class does not work. We assume that 
> our serialize function is not in the correct manner. Is it possible to pack 
> own objects using the dealii::Utilities:pack() function? What is the proper 
> way defining the serialize function?
>
> Attached you'll find our small example code. We are using deal.II 9.0.1.
>
> Thank you very much in advance.
>
> Best regards,
> Maurice Rohracker
> Master Student Computational Engineering
> FAU Erlange-Nürnberg
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"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/2b2fb018-0dfd-497a-a1df-9372f74e9e23o%40googlegroups.com.


[deal.II] Re: hp fem error assigning Fourier

2020-06-08 Thread 'peterrum' via deal.II User Group
Dear Ihsan,

is the issue solved now? I have compiled your code with the current version 
of deal.II and it works.

Peter 

On Monday, 8 June 2020 09:56:21 UTC+2, A.Z Ihsan wrote:
>
> Oops, i was wrong. I followed the deal.ii 9.2.0 tutorial meanwhile in my 
> local deal.ii version is 9.1.
> There is a couple different implementation in terms of FESeries::Fourier.
>
> On Friday, June 5, 2020 at 12:25:47 PM UTC+2, A.Z Ihsan wrote:
>>
>> Hi Peter, 
>> thank you for the answer. Actually i did put the fe_series.h. 
>> I forgot to mention that the problem arise when i use template 
>> specialization by the end the implementation
>>
>> #include 
>> #include 
>> #include 
>> #include 
>>
>> using namespace dealii;
>>
>> namespace hpfe{
>> template 
>> class HPSolver
>> {
>> public:
>>  HPSolver(
>>  const unsigned int max_fe_degree);
>>  //virtual ~HPSolver();
>>
>> const hp::FECollection& get_fe_collection();
>>  const hp::QCollection& get_face_quadrature_collection();   
>>  
>> protected:
>>  hp::FECollection fe_collection;
>>  hp::QCollection quadrature_collection;
>>  hp::QCollection face_quadrature_collection;
>>  hp::QCollection  fourier_q_collection;
>>  std::unique_ptr> fourier;
>>  std::vector ln_k;
>>  Table>fourier_coefficients;
>> };
>>
>> template 
>> HPSolver::HPSolver(
>>  const unsigned int max_degree)
>> {
>>  for (unsigned int degree=2; degree <= max_degree; ++degree)
>>{
>>  fe_collection.push_back(FE_Q(degree));
>>  quadrature_collection.push_back(QGauss(degree+1));
>>  face_quadrature_collection.push_back(QGauss(degree+1));
>>}
>>  const unsigned int N = max_degree;
>>  QGauss<1>  base_quadrature(2);
>>  QIterated quadrature(base_quadrature, N);
>>  for (unsigned int i = 0; i < fe_collection.size(); i++)
>>fourier_q_collection.push_back(quadrature);
>>  std::vector n_coefficients_per_direction(fe_collection.
>> size(), N);
>>  fourier = std::make_unique>(n_coefficients_per_
>> direction, fe_collection, fourier_q_collection);
>>  //resize(fourier_coefficients, N);
>> }
>> }
>> template class hpfe::HPSolver<3, Vector> ;
>>
>> can you try once more?
>>
>> BR, 
>> ihsan
>>
>>
>>> On Friday, 5 June 2020 10:40:29 UTC+2, A.Z Ihsan wrote:


 Hi All, 

 I am trying to implement hp-fem into my problem according to the 
 step-27. But, i have an error when i am trying to compile, 

 error: no matching function for call to 'dealii::FESeries::Fourier<3, 
> 3>::Fourier(std::vector&, dealii::hp::FECollection<3, 3>&, 
> dealii::hp::QCollection<3>&)'
>  { return unique_ptr<_Tp>(new 
> _Tp(std::forward<_Args>(__args)...)); }


 i believe there is a mistake in assigning fourier, but i copied the 
 step-27 exactly. 
 Here is the code snippet... could someone help me?

 template 
 class HPSolver
 {
 public:
 HPSolver(
 const unsigned int max_fe_degree);
 virtual ~HPSolver();

 const hp::FECollection& get_fe_collection();
 const hp::QCollection& get_face_quadrature_collection();   
  
 protected:
 hp::FECollection fe_collection;
 hp::QCollection quadrature_collection;
 hp::QCollection face_quadrature_collection;
 hp::QCollection  fourier_q_collection;
 std::unique_ptr> fourier;
 std::vector ln_k;
 Table>fourier_coefficients;
 };

 template 
 HPSolver::HPSolver(
 const unsigned int max_degree,
 :
 max_fe_degree(max_degree)
 {
 for (unsigned int degree=2; degree <= max_fe_degree; ++degree)
   {
 fe_collection.push_back(FE_Q(degree));
 quadrature_collection.push_back(QGauss(degree+1));
 face_quadrature_collection.push_back(QGauss(degree+1));
   }
 const unsigned int N = max_fe_degree;
 QGauss<1>  base_quadrature(2);
 QIterated quadrature(base_quadrature, N);
 for (unsigned int i = 0; i < fe_collection.size(); i++)
   fourier_q_collection.push_back(quadrature);
 std::vector 
 n_coefficients_per_direction(fe_collection.size(), N);
 fourier = 
 std_cxx14::make_unique>(n_coefficients_per_direction,
  
 fe_collection, fourier_q_collection);
 resize(fourier_coefficients, N);
 }

 BR, 
 Ihsan

>>>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"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/59bc6f2b-6fbc-4fd9-840f-6f491004450bo%40googlegroups.com.


[deal.II] Re: triangulation save not working for 1D domain

2020-06-06 Thread 'peterrum' via deal.II User Group
What type of triangulation are you using?

Peter

On Saturday, 6 June 2020 17:52:53 UTC+2, Amaresh B wrote:
>
>   Dear all,
>
> I am trying to save my triangulation and using the lines below. 
> 'triangulation.save' commands works and saves the mesh if my domain is 
> either 2D or 3D (i.e. dim=2 or 3). However, it does not save the 
> triangulation if dim=1, nor it gives any error message. I will be thankful 
> if someone can comment on it and give suggestions.
>
>
>
>std::vector 
> sol_transfer_vectors;
>sol_transfer_vectors.push_back();
> 
> 
> parallel::distributed::SolutionTransfer 
> sol_transfer(dof_handler);
>  sol_transfer.prepare_for_serialization 
> (sol_transfer_vectors);
>
>   triangulation.save("restart.mesh");
>
> Thank you,
>
> Regards
> Amaresh
>  
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"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/00da9e61-551a-4a19-851a-3179fc4df131o%40googlegroups.com.


[deal.II] Re: hp fem error assigning Fourier

2020-06-05 Thread 'peterrum' via deal.II User Group
Dear Ihsan,

I have no problem to compile the following code (your code with minor 
adjustments):


#include 
#include 
#include 
#include 

using namespace dealii;

template 
class HPSolver
{
public:
   HPSolver(
   const unsigned int max_fe_degree);
   //virtual ~HPSolver();

const hp::FECollection& get_fe_collection();
   const hp::QCollection& get_face_quadrature_collection();   
 
protected:
   hp::FECollection fe_collection;
   hp::QCollection quadrature_collection;
   hp::QCollection face_quadrature_collection;
   hp::QCollection  fourier_q_collection;
   std::unique_ptr> fourier;
   std::vector ln_k;
   Table>fourier_coefficients;
};

template 
HPSolver::HPSolver(
   const unsigned int max_degree)
{
   for (unsigned int degree=2; degree <= max_degree; ++degree)
 {
   fe_collection.push_back(FE_Q(degree));
   quadrature_collection.push_back(QGauss(degree+1));
   face_quadrature_collection.push_back(QGauss(degree+1));
 }
   const unsigned int N = max_degree;
   QGauss<1>  base_quadrature(2);
   QIterated quadrature(base_quadrature, N);
   for (unsigned int i = 0; i < fe_collection.size(); i++)
 fourier_q_collection.push_back(quadrature);
   std::vector n_coefficients_per_direction(fe_collection.size
(), N);
   fourier = std::make_unique>(
n_coefficients_per_direction, fe_collection, fourier_q_collection);
   //resize(fourier_coefficients, N);
}

int main()
{
   HPSolver<3,Vector > solver(3);
}


So my guess is that you have forgotten: `#include https://www.dealii.org/developer/doxygen/deal.II/fe__series_8h_source.html>
>`?

Hope this helps!

Peter

On Friday, 5 June 2020 10:40:29 UTC+2, A.Z Ihsan wrote:
>
>
> Hi All, 
>
> I am trying to implement hp-fem into my problem according to the step-27. 
> But, i have an error when i am trying to compile, 
>
> error: no matching function for call to 'dealii::FESeries::Fourier<3, 
>> 3>::Fourier(std::vector&, dealii::hp::FECollection<3, 3>&, 
>> dealii::hp::QCollection<3>&)'
>>  { return unique_ptr<_Tp>(new _Tp(std::forward<_Args>(__args)...)); }
>
>
> i believe there is a mistake in assigning fourier, but i copied the 
> step-27 exactly. 
> Here is the code snippet... could someone help me?
>
> template 
> class HPSolver
> {
> public:
> HPSolver(
> const unsigned int max_fe_degree);
> virtual ~HPSolver();
>
> const hp::FECollection& get_fe_collection();
> const hp::QCollection& get_face_quadrature_collection();   
>  
> protected:
> hp::FECollection fe_collection;
> hp::QCollection quadrature_collection;
> hp::QCollection face_quadrature_collection;
> hp::QCollection  fourier_q_collection;
> std::unique_ptr> fourier;
> std::vector ln_k;
> Table>fourier_coefficients;
> };
>
> template 
> HPSolver::HPSolver(
> const unsigned int max_degree,
> :
> max_fe_degree(max_degree)
> {
> for (unsigned int degree=2; degree <= max_fe_degree; ++degree)
>   {
> fe_collection.push_back(FE_Q(degree));
> quadrature_collection.push_back(QGauss(degree+1));
> face_quadrature_collection.push_back(QGauss(degree+1));
>   }
> const unsigned int N = max_fe_degree;
> QGauss<1>  base_quadrature(2);
> QIterated quadrature(base_quadrature, N);
> for (unsigned int i = 0; i < fe_collection.size(); i++)
>   fourier_q_collection.push_back(quadrature);
> std::vector 
> n_coefficients_per_direction(fe_collection.size(), N);
> fourier = 
> std_cxx14::make_unique>(n_coefficients_per_direction, 
> fe_collection, fourier_q_collection);
> resize(fourier_coefficients, N);
> }
>
> BR, 
> Ihsan
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"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/8312e8f9-6237-4c24-a452-6c445a609b73o%40googlegroups.com.


[deal.II] Re: LinearOperator MPI Parallel Vector

2020-04-23 Thread 'peterrum' via deal.II User Group
Dear Doug,

Could you post a short code how you want to use the LinearOperator so that 
I know what actually is not working.

Regarding Trilinos + LA::dist::Vectror: there is an open PR (
https://github.com/dealii/dealii/pull/9925) which adds the instantiations 
(hope I did not miss any).

Regarding GPU: currently there is only support for matrix-free and non for 
matrix-based algorithms. These GPU implementations use currently 
LA::dist::Vector. 

Personally, I am always using LA::dist::Vector (not just in the matrix-free 
code but also) in combination with TrilinosWrappers::SparseMatrix. It works 
very well! I have no experience how well it works with PETSc.

Peter

On Thursday, 23 April 2020 04:00:10 UTC+2, Doug Shi-Dong wrote:
>
> Hello,
>
> I have been using the LinearAlgebra::distributed::Vector class for MPI 
> parallelization since the way it works is more familiar to what I had 
> worked with and seemed more flexible.
>
> However, for parallelization, I have to either use a Trilinos or PETSc 
> matrix since the native deal.II SparseMatrix is only serial (correct me if 
> I'm wrong). Seems like I can do matrix-vector multiplications just fine 
> between LA::dist::Vector and the wrapped matrices. However, when it gets to 
> LinearOperator, it looks like a TrilinosWrappers::SparseMatrix wrapped 
> within a LinearOperator only works with a TrilinosWrappers::MPI::Vector, 
> and same thing for PETSc.
>
> I am wondering what the community is using as their go-to parallel 
> matrices and vectors, and if you've been mixing them. E.g. matrix-free with 
> Trilinos/PETSc vectors, or PETSc matrices with LA::dist::Vector. From what 
> I've seen from some tutorials, there is a way to code it up such that 
> either Trilinos or PETSc wrappers are used interchangeably, but the 
> LA::dist::Vector does not seem be nicely interchangeable with the 
> Trilinos/PETSc ones. 
>
> I was kind of hoping to be able to use LA::dist::Vector for everything, am 
> I expecting too much from it? Maybe I just need to fix the LinearOperator 
> implementation to mix-and-match the data structure? If I do commit to 
> Trilinos matrices/vectors, will I have trouble doing some matrix-free or 
> GPU stuff in the far future?
>
> Best regards,
>
> Doug
>

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


[deal.II] Re: Assembling sparse matrix from Matrix-free vmult with constrains

2020-03-10 Thread 'peterrum' via deal.II User Group
Hi Michal,

any chance that you post or send me a small runnable test program.

By the way, there is an open PR (https://github.com/dealii/dealii/pull/9343) 
for the computation of the diagonal in a matrix-free manner. Once this is 
merged, I will work the matrix-free assembly of sparse matrices.

Thanks,
Peter

On Monday, 9 March 2020 15:26:10 UTC+1, Michał Wichrowski wrote:
>
> Dear  all,
> I've got matrix-free multigrid solver for Stokes problem. The main 
> bottleneck is solution of coarse problem, so I tried to assemble the 
> regular sparse matrix and use direct solver. Since the coarse problem is 
> (relatively) small, I used vmults by unit vector to obtain columns of the 
> matrix.  This is my code:
>
>
> setup_matrix(); // generates sparsity patters and reinits matrices
>
> LinearAlgebra::distributed::BlockVector  dst;
> LinearAlgebra::distributed::BlockVector  src;
>
> src.reinit(2);
> for (unsigned int b = 0; b < 2; ++b)
> stokes_matrix.get_matrix_free()->initialize_dof_vector(
> src.block(b), b);
> src.collect_sizes();
> src =0;
>
>
> dst.reinit(2);
> dst.block(0).reinit(owned_dofs_u, relevant_dofs_u, MPI_COMM_WORLD);
> dst.block(1).reinit(owned_dofs_p, relevant_dofs_p, MPI_COMM_WORLD);
> dst.collect_sizes();
>
> for(types::global_dof_index i =0; i< owned_dofs_u.size(); ++i){
> src=0;
> dst=0;
> if(owned_dofs_u.is_element(i) )
> src.block(0)(i)=1;
>
> src.compress(VectorOperation::insert);
> stokes_matrix.vmult(dst, src);
>
> dst.update_ghost_values();
>
> for(IndexSet::ElementIterator  index_iter =owned_dofs_u.begin();
> index_iter != owned_dofs_u.end();
> ++index_iter){
> if(dst.block(0)(*index_iter)!=0 )
> block_A->set(*index_iter,i, dst.block(0)(*index_iter) );
> }
>
> }
>
> block_A->compress(VectorOperation::unknown);
>
> Without constrains the matrix matches the matrix-free operator, but with 
> constrains present it does not. What is the proper way to assemble the 
> matrix with vmult?
>
> Best,
> Michał Wichrowski
>
>
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"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/11de6867-e589-439c-8aa1-e56d808a0345%40googlegroups.com.


[deal.II] Re: Application of inverse mass matrix in matrix-free context using cell_loop instead of scale()

2020-01-18 Thread 'peterrum' via deal.II User Group
Yes, like 
here 
https://github.com/dealii/dealii/blob/b84270a1d4099292be5b3d43c2ea65f3ee005919/tests/matrix_free/pre_and_post_loops_01.cc#L100-L121

On Saturday, 18 January 2020 12:57:24 UTC+1, Maxi Miller wrote:
>
> In step-48 the inverse mass matrix is applied by moving the inverse data 
> into a vector and applying the function scale(), i.e. as in the following 
> code:
> data.cell_loop(::local_apply_cell,
>this,
>dst,
>src,
>true);
> computing_times[0] += timer.wall_time();
> timer.restart();
>
>dst.scale(inv_mass_matrix);
>
>
>
> Now I would like to do the same, but use a cell_loop instead of scale(). 
> Thus, I created an additional function called 
> "local_apply_inverse_mass_matrix" as 
> template 
> void LaplaceOperator::
> local_apply_inverse_mass_matrix(
> const MatrixFree &   data,
> LinearAlgebra::distributed::Vector &  dst,
> const LinearAlgebra::distributed::Vector ,
> const std::pair & cell_range) 
> const
> {
> (void) data;
> (void) cell_range;
> dst = src;
> dst.scale(inv_mass_matrix);
> }
>
> When running that code, though, using
> LinearAlgebra::distributed::Vector tmp(src);
>
> data.initialize_dof_vector(tmp);
> data.initialize_dof_vector(dst);
> data.cell_loop(::local_apply_cell,
>this,
>tmp,
>src,
>true);
> computing_times[0] += timer.wall_time();
> timer.restart();
>
> data.cell_loop(::local_apply_inverse_mass_matrix,
>this,
>dst,
>tmp,
>true);
> computing_times[1] += timer.wall_time();
>
> computing_times[3] += 1.;
>
>
> I get the error 
> An error occurred in line <3338> of file  matrix_free/matrix_free.h> in function
> void dealii::internal::VectorDataExchange VectorizedArrayType>::compress_start(unsigned int, VectorType&) [with 
> VectorType = dealii::LinearAlgebra::distributed::Vector; typename 
> std::enable_if<(dealii::internal::has_compress_start::value && 
> dealii::internal::has_exchange_on_subset::value), VectorType
> >::type*  = 0; int dim = 2; Number = double; 
> VectorizedArrayType = dealii::VectorizedArray]
> The violated condition was: 
> vec.has_ghost_elements() == false
> Additional information: 
> You are trying to use functionality in deal.II that is currently not 
> implemented. In many cases, this indicates that there simply didn't 
> appear much of a need for it, or that the author of the original code did 
> not have the time to implement a particular case. If you hit this 
> exception, it is therefore worth the time to look into the code to find out 
> whether you may be able to implement the missing functionality. If you do, 
> please consider providing a patch to the deal.II development sources (see 
> the deal.II website on how to contribute).
>
> Stacktrace:
> ---
> #0  ./MF_FES_RK4-Test: void dealii::internal::VectorDataExchange<2, 
> double, dealii::VectorizedArray 
> >::compress_start dealii::MemorySpace::Host>, 
> (dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>*)0>(unsigned int, 
> dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>&)
> #1  ./MF_FES_RK4-Test: void dealii::internal::compress_start<2, 
> dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>, double, dealii::VectorizedArray, 
> (dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>*)0>(dealii::LinearAlgebra::distributed::Vector  
> dealii::MemorySpace::Host>&, dealii::internal::VectorDataExchange<2, 
> double, dealii::VectorizedArray >&, unsigned int)
> #2  ./MF_FES_RK4-Test: dealii::internal::MFWorker double, dealii::VectorizedArray >, 
> dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>, 
> dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>, Step40::LaplaceOperator<2, 2, 4>, 
> true>::vector_compress_start()
> #3  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::internal::MatrixFreeFunctions::MPICommunication::execute()
> #4 
>  
> /opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
>  
>
> #5 
>  
> /opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
>  
>
> #6 
>  
> /opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
>  
>
> #7 
>  
> /opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
>  
>
> #8 
>  
> /opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
>  
>
> #9  /lib64/libpthread.so.0: 
> #10  /lib64/libc.so.6: 

[deal.II] Re: Application of inverse mass matrix in matrix-free context using cell_loop instead of scale()

2020-01-18 Thread 'peterrum' via deal.II User Group
Hi Maxi,

I guess I am not the correct person to explain you the reason for that 
assert. But what you are doing is that while calling scale you are messing 
with the ghost values (which prevents the compress step). 

You should do it only locally. What you might want to check it out are the 
new `operation_after_loop` functions provided by the cell_loops (on the 
master branch), which gives you the information local vector entry has been 
finished and you are allowed to apply e.g. the mass matrix diagonal on.

Peter

On Saturday, 18 January 2020 12:57:24 UTC+1, Maxi Miller wrote:
>
> In step-48 the inverse mass matrix is applied by moving the inverse data 
> into a vector and applying the function scale(), i.e. as in the following 
> code:
> data.cell_loop(::local_apply_cell,
>this,
>dst,
>src,
>true);
> computing_times[0] += timer.wall_time();
> timer.restart();
>
>dst.scale(inv_mass_matrix);
>
>
>
> Now I would like to do the same, but use a cell_loop instead of scale(). 
> Thus, I created an additional function called 
> "local_apply_inverse_mass_matrix" as 
> template 
> void LaplaceOperator::
> local_apply_inverse_mass_matrix(
> const MatrixFree &   data,
> LinearAlgebra::distributed::Vector &  dst,
> const LinearAlgebra::distributed::Vector ,
> const std::pair & cell_range) 
> const
> {
> (void) data;
> (void) cell_range;
> dst = src;
> dst.scale(inv_mass_matrix);
> }
>
> When running that code, though, using
> LinearAlgebra::distributed::Vector tmp(src);
>
> data.initialize_dof_vector(tmp);
> data.initialize_dof_vector(dst);
> data.cell_loop(::local_apply_cell,
>this,
>tmp,
>src,
>true);
> computing_times[0] += timer.wall_time();
> timer.restart();
>
> data.cell_loop(::local_apply_inverse_mass_matrix,
>this,
>dst,
>tmp,
>true);
> computing_times[1] += timer.wall_time();
>
> computing_times[3] += 1.;
>
>
> I get the error 
> An error occurred in line <3338> of file  matrix_free/matrix_free.h> in function
> void dealii::internal::VectorDataExchange VectorizedArrayType>::compress_start(unsigned int, VectorType&) [with 
> VectorType = dealii::LinearAlgebra::distributed::Vector; typename 
> std::enable_if<(dealii::internal::has_compress_start::value && 
> dealii::internal::has_exchange_on_subset::value), VectorType
> >::type*  = 0; int dim = 2; Number = double; 
> VectorizedArrayType = dealii::VectorizedArray]
> The violated condition was: 
> vec.has_ghost_elements() == false
> Additional information: 
> You are trying to use functionality in deal.II that is currently not 
> implemented. In many cases, this indicates that there simply didn't 
> appear much of a need for it, or that the author of the original code did 
> not have the time to implement a particular case. If you hit this 
> exception, it is therefore worth the time to look into the code to find out 
> whether you may be able to implement the missing functionality. If you do, 
> please consider providing a patch to the deal.II development sources (see 
> the deal.II website on how to contribute).
>
> Stacktrace:
> ---
> #0  ./MF_FES_RK4-Test: void dealii::internal::VectorDataExchange<2, 
> double, dealii::VectorizedArray 
> >::compress_start dealii::MemorySpace::Host>, 
> (dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>*)0>(unsigned int, 
> dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>&)
> #1  ./MF_FES_RK4-Test: void dealii::internal::compress_start<2, 
> dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>, double, dealii::VectorizedArray, 
> (dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>*)0>(dealii::LinearAlgebra::distributed::Vector  
> dealii::MemorySpace::Host>&, dealii::internal::VectorDataExchange<2, 
> double, dealii::VectorizedArray >&, unsigned int)
> #2  ./MF_FES_RK4-Test: dealii::internal::MFWorker double, dealii::VectorizedArray >, 
> dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>, 
> dealii::LinearAlgebra::distributed::Vector dealii::MemorySpace::Host>, Step40::LaplaceOperator<2, 2, 4>, 
> true>::vector_compress_start()
> #3  /opt/dealii/lib/libdeal_II.g.so.9.2.0-pre: 
> dealii::internal::MatrixFreeFunctions::MPICommunication::execute()
> #4 
>  
> /opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
>  
>
> #5 
>  
> /opt/intel/compilers_and_libraries_2019.5.281/linux/tbb/lib/intel64_lin/gcc4.7/libtbb_debug.so.2:
>  
>
> #6 
>  
>