Dear Wolfgang, Thank you very much for your reply, and for clearing things up. Sorry if my original message wasn't clear.
I am solving a system of time-dependent equations, where time is discritised using the Euler method. One of the equations is hyperbolic, and I am considering whether to use the method of characteristics to solve this. I am not sure exactly which method I am going to use to solve the system as yet, as it depends on which method is more computationally efficient. I have two options: 1. At each timepoint solve the equations on the current mesh (using DGQ to stabilise the hyperbolic term), to give the updated "solution" at the current mesh points. As the solution does not move with the mesh, I would then have to move the mesh, and then interpolate the solution vector onto the "new" moved mesh. 2. At each timepoint solve one of the equations using the method of characteristics; in this case, part of the BlockVector solution (i.e. block 0) advects with the mesh, and some (i.e. block 1) doesn't. Therefore, after solving, I will have part of my solution defined on the "new" mesh, and some defined on the "old" mesh. I think I will go with method 1 to begin with, as then I don't have to worry about different bits of the BlockVector solution being defined on different meshes. My problem is that I need to interpolate the solution onto the new mesh - but the two meshes don't necessarily have the same boundary (consider a domain that is advecting with constant velocity, for example, in which case the two meshes would overlap) and then I don't think I can use FEFieldFunction? It may be that I am trying to do something that is beyond the current abilities of deal.II. I would be grateful of any light you could throw on the above problem! Many Thanks, Katie Katie Leonard DPhil student in Computational Biology, The University of Oxford. ________________________________________ From: [email protected] [[email protected]] on behalf of Wolfgang Bangerth [[email protected]] Sent: 06 June 2012 15:39 To: [email protected] Subject: Re: [deal.II] A moving boundary problem > I'm trying to write a code in deal.II to solve a moving boundary problem, > where > the boundary moves according to the velocity field $\mathbf{v}$. After solving > at each timestep, I want to interpolate the "old_solution" onto the new grid. > However, I am unsure about how to go about this. I can move the mesh as in > step-18, but this does not include interpolation of the solution? I was > thinking of using FEFieldFunction, but triangulation2 needs to be entirely > included within triangulation1, but if the boundary is moving, this may > not be the case. You are trying to implement something for which I gather you don't yet quite know what it is. You first need to figure out what it is that you really want, and then you will also know how to implement it. Specifically: - In the step-18 way, we assume that the quantities that we move from the old to the new mesh indeed "advect" with the mesh. For example, if you imagine a solid where in every point you have a quantity that's attached to the solid (e.g., crystal defect density; chemical composition; imprinted stress) then these are things that need to deform with the mesh if the mesh deforms with the body. - If you were to use the FEFieldFunction route, then you assume that the quantity that you had on the old mesh remains where it is while the body and associated mesh move. This would be the case if what you want to transfer is a quantity that does not move along with the material. Examples would be the intensity of an external X-ray beam going through tissue, etc. In such a case, you want to interpolate the value of this external field from the old mesh to the new mesh using exactly the same coordinate values (and not advected coordinates). In this case you will also have to figure out what the values should be at points that are now part of the domain but that previously weren't. Does this taxonomy help? > I tried using something like > > ------------------------------------------------------------------------------------------------------------------------------------------ > > BlockVector<double> x_sol = solution; > SolutionTransfer<dim, BlockVector<double> > solution_transfer(dof_handler); > > for (typename DoFHandler<dim>::active_cell_iterator > cell = dof_handler.begin_active (); > cell != dof_handler.end(); ++cell) > { > for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v) > { > if (vertex_touched[cell->vertex_index(v)] == false) > { > vertex_touched[cell->vertex_index(v)] = true; > Point<dim> vertex_displacement; > for (unsigned int d=0; d<dim; ++d) > { > vertex_displacement[d] = > time_step*solution(cell->vertex_dof_index(v,d)); > } > cell->vertex(v) += vertex_displacement; > } > } > } > > solution_transfer.interpolate(x_sol, solution); > > > -------------------------------------------------------------------------------------------------------------------- > > but even if I redistributed dofs etc, I still get the error > > void dealii::SolutionTransfer<dim, VECTOR, DH>::interpolate(const VECTOR&, > VECTOR&) const [with int dim = 2, VECTOR = dealii::BlockVector<double>, DH = > dealii::DoFHandler<2>] > The violated condition was: > in.size()==n_dofs_old > The name and call sequence of the exception was: > ExcDimensionMismatch(in.size(), n_dofs_old) > Additional Information: > Dimension 3202 not equal to 0 > > which I don't really understand, as my "n_dofs_old" should have size 3202 > as well, as the mesh is just moved! May be I have missed something quite simple. I think you are missing a call to one of the SolutionTransfer::prepare_* functions. Best W. ------------------------------------------------------------------------ Wolfgang Bangerth email: [email protected] www: http://www.math.tamu.edu/~bangerth/ _______________________________________________ dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii _______________________________________________ dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
