From: Roy S. <roy...@ic...> - 2011-08-26 03:30:33
|
On Thu, 25 Aug 2011, rochan wrote: > Before trying the L2 projection I am trying to understand how get the > gradient of my solution at a particular point (where it is defined even for a > piecewise linear basis function of the original variable). The solution is > correctly obtained using the following: > > int main (int argc, char** argv) > { > LibMeshInit init (argc, argv); > > Mesh mesh(2); > > GmshIO(mesh).read("gmsh_ex3.1a.msh"); > mesh.prepare_for_use(); > mesh.print_info(); > > EquationSystems equation_systems(mesh); > equation_systems.add_system<LinearImplicitSystem> ("Elasticity"); > LinearImplicitSystem& system_elasticity = > equation_systems.get_system<LinearImplicitSystem>("Elasticity"); > > system_elasticity.add_variable("u1", FIRST); > system_elasticity.add_variable("u2", FIRST); > > system_elasticity.attach_assemble_function(assemble_elasticity); > system_elasticity.init(); > > equation_systems.parameters.set<unsigned int>("linear solver maximum > iterations") = 1500; > equation_systems.parameters.set<Real>("linear solver tolerance") = > 1.0e-10; > equation_systems.print_info(); > system_elasticity.solve(); > > printf("Solution_complete"); > > Now I want the global solution in a vector, // Use: NumericVector<Number> &distributed_solution_vector = *system_elasticity.solution; // or: std::vector<Number> serializing_your_solution_is_inefficient; system_elasticity.update_global_solution(serializing_your_solution_is_inefficient, 0); if (libMesh::processor_id == 0) Rochan::use_solution(); else Rochan::twiddle_thumbs(); > but the following is not working > > std::vector<Number> global_soln; > std::vector<std::string> names; > equation_systems.get_solution( global_soln, names ); > printf("%d",global_soln.size()); Wow, that's kind of a misleading API name, isn't it? Jason was writing a "get_discontinuous_constant_solution_components_only()", which if given such a name would have been more clearly identifiable as a super-specialized API that nobody would want to use directly. I wasn't watching the svn diffs closely, so I didn't notice that he'd called it "get_solution" instead... > Supppose I can get the solution vector, would the following that > uses Meshfunction work? > > AutoPtr<NumericVector<Number> > fine_soln = > NumericVector<Number>::build(); > AutoPtr<MeshFunction> gradient_values_u1; > fine_soln->init (global_soln.size(), true, SERIAL); > (*fine_soln) = global_soln; > gradient_values_u1 = AutoPtr<MeshFunction> > (new MeshFunction(equation_systems, > *fine_soln, > system_elasticity.get_dof_map(), > system_elasticity.variable_number("u1"))); > Gradient sol_gradient; > gradient_values_u1->init(); > > then later when looping over quadrature points > > sol_gradient = gradient_values_u1->gradient(gp_face_xyz[gp]); This is closer, but is unnecessary. It's even less efficient under the hood than the simplest API, the one I was trying to suggest on Wednesday: // outside your loops unsigned int u1_var = system_elasticity.variable_number("u1"); // when looping over quadrature points Gradient sol_gradient = system_elasticity.point_gradient(u1_var, your_xyz[qp]); --- Roy |