Dear all,
In solving for linear elasticity, I would need to know the reaction forces.
This should be relatively easy - just use the stiffness matrix * solution
vector and we can get the reaction forces.
I used a new system to store the unpenalized stiffness matrix as well as
the force values:
// Initialize another system to store forces
LinearImplicitSystem & force_system =
equation_systems.add_system<LinearImplicitSystem> ("ForceSystem");
force_system.add_variable("fx");
force_system.add_variable("fy");
force_system.add_variable("fz");
I used `system` being the main equation system. During the assembly, I just
assembled two copies of the stiffness matrix, so I can leave one
unpenalized:
force_system.matrix->add_matrix (Ke, dof_indices);
dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
Lastly, I just used a vector_mult to get the actual reaction force:
force_system.matrix->vector_mult(*force_system.solution, *system.solution);
It worked quite conveniently, because I can also export the nodal forces
that are already associated with the mesh. However, I did see a warning:
Matrix A must be assembled before calling PetscVector::add_vector(v, A).
Please update your code, as this warning will become an error in a future
release.src/numerics/petsc_vector.C, line 235, compiled Sep 19 2018 at
23:21:53 ***
*** Warning, This code is deprecated, and likely to be removed in future
library versions! src/numerics/petsc_vector.C, line 236, compiled Sep 19
2018 at 23:21:53 ***
Any ideas in how I can get around this issue?
Best,
Shawn
--
Yuxiang "Shawn" Wang, PhD
yw...@virginia.edu
+1 (434) 284-0836
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
