Hi, Bruno:
Again, thank you so much! But I still could not figure out what is going
wrong, this is my code, could you please take a look. In order to save your
time,
I have abreviate it a bit.
template <int dim, int fe_degree>
class GravityDamUnderSeismic
{
public:
GravityDamUnderSeismic (const FiniteElement<dim> &finite_element);
void run ();
private:
void make_triangulation ();
void setup_system ();
void assemble_system ();
void solve_modes ();
void output_modes_results () const;
void compute_general_mass_and_stiffness ();
void compute_Taux_and_Tauy ();
void readin_ground_acceleration ();
void prepare_time_iteration ();
void start_time_iteration (double gamma, double beta);
void output_displacement_solutions (unsigned int current_time_step);
void output_stress_solutions (unsigned int current_time_step);
double Elasticity_Tensor(unsigned int i,unsigned int j,
unsigned int k,unsigned int l);
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
ConstraintMatrix constraints;
const unsigned int eigen_numbers;
const unsigned int time_steps;
PETScWrappers::SparseMatrix stiffness_matrix,
mass_matrix;
std::vector<PETScWrappers::MPI::Vector> eigenfunctions;
std::vector<double> eigenvalues;
std::vector<double> GeneralMass;
std::vector<double> GeneralStiffness;
std::vector<double> Tau_x, Tau_y;
std::vector<double> Ag_x, Ag_y;
Vector<double> displacement_solutions;
};
template <int dim>
class StressPostprocessor: public DataPostprocessor<dim>
{
public:
StressPostprocessor();
virtual void
compute_derived_quantities_vector (const std::vector<Vector<double> >
&uh,
const std::vector<std::vector<
Tensor<1,dim> > > &duh,
const std::vector<std::vector<
Tensor<2,dim> > > &dduh,
const std::vector<Point<dim> >
&normals,
const std::vector<Point<dim> >
&evaluation_points,
std::vector<Vector<double> >
&computed_quantities) const;
virtual std::vector<std::string> get_names () const;
virtual std::vector< DataComponentInterpretation::
DataComponentInterpretation>
get_data_component_interpretation () const;
virtual UpdateFlags get_needed_update_flags () const;
};
template <int dim>
std::vector<std::string> StressPostprocessor<dim>::get_names() const
{
std::vector<std::string> stress_solution_names;
stress_solution_names.push_back ("sigma_11");
stress_solution_names.push_back ("sigma_22");
stress_solution_names.push_back ("sigma_12");
stress_solution_names.push_back ("principle stress");
return stress_solution_names;
}
template <int dim>
std::vector<DataComponentInterpretation::DataComponentInterpretation>
StressPostprocessor<dim>::
get_data_component_interpretation () const
{
std::vector<DataComponentInterpretation::DataComponentInterpretation>
interpretation (3,
DataComponentInterpretation::component_is_part_of_vector
);
interpretation.push_back (DataComponentInterpretation::
component_is_scalar);
return interpretation;
}
template <int dim>
UpdateFlags
StressPostprocessor<dim>::get_needed_update_flags() const
{
return update_gradients | update_q_points;
}
template <int dim>
void StressPostprocessor<dim>::compute_derived_quantities_vector
(const std::vector<Vector<double> > &uh,
const std::vector<std::vector<Tensor<1,dim> > > &duh,
const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
const std::vector<Point<dim> > &/*normals*/,
const std::vector<Point<dim> > &
/*evaluation_points*/,
std::vector<Vector<double> > &
computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
Assert (duh.size() == n_quadrature_points,
ExcInternalError());
Assert (computed_quantities.size() == n_quadrature_points,
ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
for (unsigned int k = 0; k != dim; ++k)
{
computed_quantities[q][0] += lambda * duh[q][k][k] + 2 * mu *
duh[q][0][0];
computed_quantities[q][1] += lambda * duh[q][k][k] + 2 * mu *
duh[q][1][1];
}
computed_quantities[q][2] = mu * (duh[q][0][1] + duh[q][1][0]);
computed_quantities[q][3] = (computed_quantities[q][0] +
computed_quantities[q][1])/2
+ std::sqrt(std::pow((computed_quantities[q][0] -
computed_quantities[q][1])/2, 2)\
+ std::pow(computed_quantities[q][2], 2));
}
}
template <int dim, int fe_degree>
GravityDamUnderSeismic<dim, fe_degree>::GravityDamUnderSeismic(const
FiniteElement<dim> &finite_element)
:
fe (finite_element, dim),
dof_handler (triangulation),
eigen_numbers(4),
time_steps(2000)
{}
template <int dim, int fe_degree>
void GravityDamUnderSeismic<dim, fe_degree>::make_triangulation(){... ...}
template <int dim, int fe_degree>
void GravityDamUnderSeismic<dim, fe_degree>::setup_system(){... ...}
template <int dim, int fe_degree>
void GravityDamUnderSeismic<dim, fe_degree>::assemble_system(){... ...}
template <int dim, int fe_degree>
void GravityDamUnderSeismic<dim, fe_degree>::solve_modes(){... ...}
template <int dim, int fe_degree>
void GravityDamUnderSeismic<dim, fe_degree>::output_modes_results() const
{... ...}
... ... ...
... ... ...
template <int dim, int fe_degree>
void GravityDamUnderSeismic<dim, fe_degree>::start_time_iteration (double
gamma, double beta)
{
...
...
for (unsigned int i = 1; i != time_steps; ++i)
{
...
...
output_displacement_solutions(current_step);
output_stress_solutions(current_step);
}
template <int dim, int fe_degree>
void GravityDamUnderSeismic<dim, fe_degree>::output_displacement_solutions
(unsigned int current_time_step)
{
std::vector<DataComponentInterpretation::DataComponentInterpretation>
data_component_interpretation (dim, DataComponentInterpretation::
component_is_part_of_vector);
DataOut<dim> u_data;
std::string displacement_filename = "displacement_solution_"
+ Utilities::int_to_string(
current_time_step)
+ ".vtk";
u_data.attach_dof_handler (dof_handler);
u_data.add_data_vector (displacement_solutions,
std::string("displacement_solution_")
+ Utilities::int_to_string(current_time_step
),
DataOut<dim>::type_dof_data,
data_component_interpretation);
std::ofstream u_output_stream (displacement_filename.c_str());
u_data.write_vtk(u_output_stream);
u_output_stream.close();
}
template <int dim, int fe_degree>
void GravityDamUnderSeismic<dim, fe_degree>::output_stress_solutions(
unsigned int current_time_step)
{
StressPostprocessor<dim> stress_out;
DataOut<dim> stress_data;
std::string stress_filename = "stress_solution_"
+ Utilities::int_to_string(
current_time_step)
+ ".vtk";
stress_data.attach_dof_handler (dof_handler);
stress_data.add_data_vector(displacement_solutions, stress_out);
stress_data.build_patches ();
std::ofstream stess_output_stream (stress_filename.c_str());
stress_data.write_vtk (stess_output_stream);
stess_output_stream.close();
}
template <int dim, int fe_degree>
void GravityDamUnderSeismic<dim, fe_degree>::run()
{
double gamma = 0.5, beta = 1./6.;
make_triangulation ();
setup_system ();
assemble_system ();
solve_modes ();
output_modes_results ();
compute_general_mass_and_stiffness ();
compute_Taux_and_Tauy ();
readin_ground_acceleration ();
prepare_time_iteration ();
start_time_iteration (gamma, beta);
}
int main (int argc, char **argv)
{
...
...
gravity_dam_under_seismic.run ();
}
Great thanks!
Best regards,
David
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
