From: Robert <libmesh@ro...>  20140606 08:36:46

Thanks for the hint. I have now the following: for (unsigned int t_step=0; t_step != n_timesteps; ++t_step) { system.time += deltat; equation_systems.reinit(); system.solve(); } But in my FEMSystem::element_time_derivative(...) I have to do the following to get the correct displacements for the new timestep: for (unsigned int i = 0; i < 3; ++i) { const std::vector<dof_id_type>& u_idx = context.get_dof_indices(i); DenseVector<Number> Uc(u_idx.size()); context.get_system().get_dof_map().extract_local_vector(*context.get_system().current_local_solution,u_idx,Uc); for (unsigned int j = 0; j < u_idx.size(); ++j) { U[j][i] = Uc.el(j); } } instead of using just context.get_elem_solution(i).el(j); Is there an easier way for doing this? Robert On Fri, May 30, 2014 at 08:47:28AM 0500, Roy Stogner wrote: > > On Fri, 30 May 2014, Robert wrote: > > > Are there any additional step necessary in each timestep to apply > > the boundary conditions or what else am I missing? > > One little step, one big one: > > To apply the boundary conditions, you'll want to run > EquationSystems::reinit() before each time step. > > Don't trust the accuracy without testing. If you're not updating > System::time to match the *end* of your time step rather than the > beginning when reinit gets called, then your functions will get called > with the wrong time and you'll have an O(deltat) error from that. >  > Roy 