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
------------------------------------------------------------------------------
EMC VNX: the world's simplest storage, starting under $10K
The only unified storage solution that offers unified management
Up to 160% more powerful than alternatives and 25% more efficient.
Guaranteed. http://p.sf.net/sfu/emc-vnx-dev2dev
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
