date:20170323

[deal.II] making a package of PETSc vectors

2017-03-23 Thread itomas

Dear everybody:

I don't have any particular problem, just trying to streamline a code. I am 
have a mess with the interfaces of my code, in particular with the 
functions and classes I have created, primarily because I have to handle a 
''packet'' of vectors, say 

PETScWrappers::MPI::Vector  Vect1;
PETScWrappers::MPI::Vector  Vect2;
PETScWrappers::MPI::Vector  Vect3;
PETScWrappers::MPI::Vector  Vect4;
…
PETScWrappers::MPI::Vector  Vect7;

and use them either as input or output arguments of some functions, 
solvers, complementary computations, etc. This is quite cumbersome, and as 
solution I naively thought I could do the following

std::vector  VectorPacket ;
VectorPacket.resize(7) ;
for (int i=0; i<7 ; ++i) 
  VectorPacket[i].reinit (locally_owned_dofs, mpi_communicator ) ;

So, now VectorPacket is really a packet, rather than seven isolated 
vectors. It is almost needless to say that this is bad idea, std::vector 
will demand quite a few things from the class PETScWrappers::MPI::Vector 
which are just not there (in particular some operators). In reality the 
size of the packet of vector is not always seven, which means that I really 
need to use a proper container for which I can change it's size. Do you 
have any suggestion (a.k.a. creative approach) on how to deal with an issue 
like this one? Which container from the STL would be appropriate?

Many thanks.

-- 
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: Vertices constraints inside the domain

2017-03-23 Thread seven

Guess I figured it out. I didn't do 
constraints.distribute();


在 2017年3月23日星期四 UTC-4上午1:38:32，seven写道：
>
>
> Hello all,
>
> I am trying to project a solution in FE_DGP(sol_dg) to FE_Q (sol_cg), but 
> I only know the gradients of sol_dg and its values at some vertices, so I 
> want to fix values of sol_cg at the same vertices after the projection.
>
> The first method I tried is quite straightforward, but the resulting 
> global matirx is very ill-conditioned. I set the corresponding row to zero 
> except the diagonal entry, and then change the right hand size.
>
> for( ; cell_c!= endc;++cell_dg, ++cell_c ){
> const unit active_index = cell_dg->active_cell_index();
> if(interface_flag[active_index] == 1){
> fe_values_continuous.reinit(cell_c);
> fe_values_v_dg.reinit(cell_dg);
> fe_values_v_dg.get_function_values(solution, solu_dg);
> cell_c->get_dof_indices(local_dof_indices_c);
> for(unit vert=0; vert<4; ++vert){
> // local_dof_indices_c[vert] global index
> const unit index_i = local_dof_indices_c[vert];
> ls_rhs_continuous(index_i) = solu_dg[vert]* aver_diag;
> for(unit i=0; i ls_matrix_continuous.set(local_dof_indices_c[vert], i, 0.);
> ls_matrix_continuous.set(index_i, index_i, aver_diag);
> }
> }
>
> } 
> so I want to use constraint matrix to do the same thing, but cannot get 
> the right answer. 
>
> for( ; cell_c!= endc;++cell_dg, ++cell_c ){
> const unit active_index = cell_dg->active_cell_index();
> if(interface_flag[active_index] == 1){
> fe_values_continuous.reinit(cell_c);
> fe_values_v_dg.reinit(cell_dg);
> fe_values_v_dg.get_function_values(solution, solu_dg);
> cell_c->get_dof_indices(local_dof_indices_c);
> for(unit vert=0; vert<4; ++vert){
> // local_dof_indices_c[vert] global index
> const unit index_i = local_dof_indices_c[vert];
> constraints.add_line(index_i);
> constraints.set_inhomogeneity(index_i, solu_dg[vert]);
> }
> }
> }
> constraints.close();
>
> cell_c = dof_handler_continuous.begin_active();
> cell_dg = dof_handler.begin_active();
> for( ; cell_c!=endc; ++cell_c, ++cell_dg){
> fe_values_dg.reinit(cell_dg);
> fe_values_dg.get_function_gradients(solution, grad_dg);
>
> ls_matrix = 0;
> ls_rhs = 0;
> fe_values_continuous.reinit(cell_c);
>
> for(unit q=0; q for(unit i=0; i ls_rhs(i) += fe_values_continuous.shape_grad(i,q)
>   * fe_values_continuous.JxW(q)
>   * grad_dg[q];
> for(unit j=0; j ls_matrix(i,j) += fe_values_continuous.shape_grad(i,q)
>* fe_values_continuous.shape_grad(j,q)
>* fe_values_continuous.JxW(q);
> }
> }
> cell_c->get_dof_indices(local_dof_indices_c);
> constraints.distribute_local_to_global (ls_matrix, ls_rhs,
> local_dof_indices_c,
> ls_matrix_continuous, 
> ls_rhs_continuous);
>
> Any help would be appreciated
>
> Thanks,
> Jiaqi 
>

-- 
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] fixing one component of solution to the same value

2017-03-23 Thread RAJAT ARORA

Thanks Professor.


On Wednesday, March 15, 2017 at 3:35:42 PM UTC-4, Wolfgang Bangerth wrote:
>
>
> > I am using deal.ii to solve a 3D solid mechanics problem which uses 
> > p::d::triangulation. 
> > 
> > I am solving equilibrium equations to solve for the displacement in the 
> domain. 
> > The domain of the problem is a cylinder with the z-axis aligned along 
> the axis 
> > of the cylinder. 
> > 
> > To avoid the warping, I want constraints such that 
> > displacement in the z direction is same for all nodes on the +z surface. 
>
> Choose one degree of freedom for the vertical displacement on one 
> processor 
> (say, it is u_42), and then on all other processors you want to go through 
> all 
> vertical displacements at the boundary in question and set 
>u_i = u_42 
> This is a straightforward constraint to add to the ConstraintMatrix class. 
> Take a look at step-11 to see how constraints are added. 
>
> Best 
>   W. 
>
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

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

Re: [deal.II] How can I resize a std::vector<Functions::ParsedFunction>?

2017-03-23 Thread luca.heltai

Resizing is fine. It is the equality that won’t work:

function_pointers[f] = 
 std::shared_ptr(new 
dealii::Functions::ParsedFunction()));

is the correct way to do it. 

L.


> On 23 Mar 2017, at 15:09, Alex Zimmerman  wrote:
> 
> Also for the record, here's a copy of issue_10_push_back.cc from that 
> repository, which is the example of how to solve my issue:
> 
> #include 
> 
> #include 
> #include  
> 
> int main(int /*argc*/, char** /*argv*/)
> {
> const unsigned int dim = 2;
> const unsigned int function_count = 4;
> 
> std::vector>
> function_pointers;
> 
> dealii::ParameterHandler prm;
> dealii::Point point;
> double value;
> 
> for (unsigned int f = 0; f < function_count; ++f)
> {
> function_pointers.push_back(
> std::shared_ptr(new 
> dealii::Functions::ParsedFunction()));
> 
> prm.enter_subsection("parsed_function_"+std::to_string(f));
> {
> function_pointers[f]->declare_parameters(prm);
> function_pointers[f]->parse_parameters(prm);
> }
> prm.leave_subsection();
> 
> function_pointers[f]->value(point, value);
> 
> std::cout << "f(" << point << ") = " << value << std::endl;
> }
> 
> return 0;
> }
> 
> 
> On Thursday, March 23, 2017 at 3:05:53 PM UTC+1, Alex Zimmerman wrote:
> Timo, thanks for the extra clarification. As I mentioned in my reply to Luca, 
> I don't understand why declare_parameters appeared to succeed in my test, 
> which failed at parse_parameters.
> 
> My issue is solved using the vector.push_back(new 
> ParsedFunction) approach. My attempt at the resize approach isn't working; 
> but push_back is fine. 
> 
> Related test code is here: 
> https://github.com/alexanderzimmerman/nsb-pcm/tree/bugs/tests
> ctest passes all tests except for issue_10_resize.cc. I'm guessing that I 
> have the syntax wrong somewhere.
> 
> Thanks for the help!
> 
> 
> On Wednesday, March 22, 2017 at 10:01:21 PM UTC+1, Timo Heister wrote:
> Your call to resize() will create 4 shared_ptrs, but they are pointing 
> to nothing (NULL) because you are never allocating any objects and 
> assigning them. 
> 
> Think of this example: 
> std::vector v; 
> v.resize(4); 
> now v[0] is a pointer to an int, but it is NULL unless you do something like 
> v[0] = new int; 
> 
> 
> 
> On Wed, Mar 22, 2017 at 3:09 PM, Alex Zimmerman 
>  wrote: 
> > Also for the record, I verified that the test passes if I use a 
> > std::vector without resizing: 
> > https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_alexanderzimmerman_nsb-2Dpcm_blob_bugs_tests_bug-5F10-5Fpass.cc=DwIFaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=bMNkhuzVlxuFjrw_sjeYQVtM1sMLS6wmKIkP9ILH7d4=v1-30ZnlLX9dscEMGIPAxToOfyXrgOcTjUUS54eVdJA=
> >  
> > 
> > 
> > 
> > On Wednesday, March 22, 2017 at 7:17:25 PM UTC+1, Alex Zimmerman wrote: 
> >> 
> >> Luca, 
> >> 
> >> Thanks for the idea. Conceptually I don't understand how a shared_ptr 
> >> solves my problem. I can resize 
> >> std::vector>; but 
> >> all of the shared pointers are empty. I imagine that they are supposed to 
> >> point to ParsedFunction objects that I've instantiated somewhere. 
> >> 
> >> That being said, I think you are leading me in the correct direction. In 
> >> the following test code, "default.prm" is successfully generated with the 
> >> four parsed function sections. It seg faults when I call parse_parameters; 
> >> so I'm still debugging: 
> >> 
> >> #include  
> >> 
> >> #include  
> >> #include  
> >> 
> >> int main(int /*argc*/, char** /*argv*/) 
> >> { 
> >> const unsigned int dim = 2; 
> >> const unsigned int function_count = 4; 
> >> 
> >> std::vector> 
> >> function_pointers; 
> >> 
> >> function_pointers.resize(function_count); 
> >> 
> >> dealii::ParameterHandler prm; 
> >> 
> >> for (unsigned int f = 0; f < function_count; ++f) 
> >> { 
> >> prm.enter_subsection("parsed_function_"+std::to_string(f)); 
> >> { 
> >> function_pointers[f]->declare_parameters(prm); 
> >> } 
> >> prm.leave_subsection(); 
> >> } 
> >> 
> >> prm.read_input("default.prm", false, true); 
> >> 
> >> for (unsigned int f = 0; f < function_count; ++f) 
> >> { 
> >> prm.enter_subsection("parsed_function_"+std::to_string(f)); 
> >> { 
> >> function_pointers[f]->parse_parameters(prm); 
> >> } 
> >> prm.leave_subsection(); 
> >> } 
> >> 
> >> return 0; 
> >> } 
> >> 
> >> Thanks, 
> >> 
> >> Alex 
> >> 
> >> On Wednesday, March 22, 2017 at

Re: [deal.II] How can I resize a std::vector<Functions::ParsedFunction>?

2017-03-23 Thread Alex Zimmerman

Timo, thanks for the extra clarification. As I mentioned in my reply to 
Luca, I don't understand why declare_parameters appeared to succeed in my 
test, which failed at parse_parameters.

My issue is solved using the 
vector.push_back(new ParsedFunction) approach. 
My attempt at the resize approach isn't working; but push_back is fine. 

Related test code is here: 
https://github.com/alexanderzimmerman/nsb-pcm/tree/bugs/tests
ctest passes all tests except for issue_10_resize.cc 
.
 
I'm guessing that I have the syntax wrong somewhere.

Thanks for the help!


On Wednesday, March 22, 2017 at 10:01:21 PM UTC+1, Timo Heister wrote:
>
> Your call to resize() will create 4 shared_ptrs, but they are pointing 
> to nothing (NULL) because you are never allocating any objects and 
> assigning them. 
>
> Think of this example: 
> std::vector v; 
> v.resize(4); 
> now v[0] is a pointer to an int, but it is NULL unless you do something 
> like 
> v[0] = new int; 
>
>
>
> On Wed, Mar 22, 2017 at 3:09 PM, Alex Zimmerman 
>  wrote: 
> > Also for the record, I verified that the test passes if I use a 
> > std::vector without resizing: 
> > 
> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_alexanderzimmerman_nsb-2Dpcm_blob_bugs_tests_bug-5F10-5Fpass.cc=DwIFaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=bMNkhuzVlxuFjrw_sjeYQVtM1sMLS6wmKIkP9ILH7d4=v1-30ZnlLX9dscEMGIPAxToOfyXrgOcTjUUS54eVdJA=
>  
> > 
> > 
> > 
> > On Wednesday, March 22, 2017 at 7:17:25 PM UTC+1, Alex Zimmerman wrote: 
> >> 
> >> Luca, 
> >> 
> >> Thanks for the idea. Conceptually I don't understand how a shared_ptr 
> >> solves my problem. I can resize 
> >> std::vector>; 
> but 
> >> all of the shared pointers are empty. I imagine that they are supposed 
> to 
> >> point to ParsedFunction objects that I've instantiated somewhere. 
> >> 
> >> That being said, I think you are leading me in the correct direction. 
> In 
> >> the following test code, "default.prm" is successfully generated with 
> the 
> >> four parsed function sections. It seg faults when I call 
> parse_parameters; 
> >> so I'm still debugging: 
> >> 
> >> #include  
> >> 
> >> #include  
> >> #include  
> >> 
> >> int main(int /*argc*/, char** /*argv*/) 
> >> { 
> >> const unsigned int dim = 2; 
> >> const unsigned int function_count = 4; 
> >> 
> >> 
> std::vector> 
> >> function_pointers; 
> >> 
> >> function_pointers.resize(function_count); 
> >> 
> >> dealii::ParameterHandler prm; 
> >> 
> >> for (unsigned int f = 0; f < function_count; ++f) 
> >> { 
> >> prm.enter_subsection("parsed_function_"+std::to_string(f)); 
> >> { 
> >> function_pointers[f]->declare_parameters(prm); 
> >> } 
> >> prm.leave_subsection(); 
> >> } 
> >> 
> >> prm.read_input("default.prm", false, true); 
> >> 
> >> for (unsigned int f = 0; f < function_count; ++f) 
> >> { 
> >> prm.enter_subsection("parsed_function_"+std::to_string(f)); 
> >> { 
> >> function_pointers[f]->parse_parameters(prm); 
> >> } 
> >> prm.leave_subsection(); 
> >> } 
> >> 
> >> return 0; 
> >> } 
> >> 
> >> Thanks, 
> >> 
> >> Alex 
> >> 
> >> On Wednesday, March 22, 2017 at 4:45:03 PM UTC+1, Luca Heltai wrote: 
> >>> 
> >>> Hi Alex, 
> >>> 
> >>> I’d use a vector of 
> >>> 
> >>> std::vector > v; 
> >>> 
> >>> which are light objects, and can be resized and reshaped. Whenever you 
> do 
> >>> a push_back, you’d have to use 
> >>> 
> >>> v.push_back(std_cxx11::shared_pointer >>> ParsedFunction(…) ); 
> >>> 
> >>> then 
> >>> 
> >>> v[i]->value(…) 
> >>> 
> >>> would work. 
> >>> 
> >>> 
> >>> > On 22 Mar 2017, at 16:31, Alex Zimmerman  
> wrote: 
> >>> > 
> >>> > Fundamentally I am trying to allow for a variable number of 
> >>> > ParsedFunction objects to be specified in a parameter input file. 
> Maybe 
> >>> > there is a better approach which circumvents my issue below. I can 
> continue 
> >>> > my work for some time with this being a constant; but as soon as I 
> want to 
> >>> > extend to 3D or different geometries, this will be a roadblock. 
> >>> > 
> >>> > I've documented my issue on my GitHub repostiory: 
> >>> > 
> https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_alexanderzimmerman_nsb-2Dpcm_issues_10=DwIFaQ=Ngd-ta5yRYsqeUsEDgxhcqsYYY1Xs5ogLxWPA_2Wlc4=4k7iKXbjGC8LfYxVJJXiaYVu6FRWmEjX38S7JmlS9Vw=bMNkhuzVlxuFjrw_sjeYQVtM1sMLS6wmKIkP9ILH7d4=a-21yeSa2EGcnF7p79Efq_xBzJ8qavKiAtJHROzgbX0=
>  
>  . Here's a copy of 
> >>> > the text from my issue, which pretty much explains

[deal.II] making a package of PETSc vectors

[deal.II] Re: Vertices constraints inside the domain

Re: [deal.II] fixing one component of solution to the same value

Re: [deal.II] How can I resize a std::vector<Functions::ParsedFunction>?

Re: [deal.II] How can I resize a std::vector<Functions::ParsedFunction>?

5 matches

Site Navigation

Mail list logo

Footer information