diff --git a/SAND_top_opt/CMakeLists.txt b/SAND_top_opt/CMakeLists.txt new file mode 100644 index 00000000..c416efee --- /dev/null +++ b/SAND_top_opt/CMakeLists.txt @@ -0,0 +1,40 @@ +# Set the name of the project and target: +SET(TARGET "SAND") + + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. + + +FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +FILE(GLOB_RECURSE TARGET_INC "include/*.h") +SET( TARGET_SRC ${TARGET_SRC} ${TARGET_INC} + ) + + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) +FIND_PACKAGE(deal.II 9.4.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/SAND_top_opt/doc/author b/SAND_top_opt/doc/author new file mode 100644 index 00000000..811986c0 --- /dev/null +++ b/SAND_top_opt/doc/author @@ -0,0 +1 @@ +Justin O'Connor diff --git a/SAND_top_opt/doc/builds-on b/SAND_top_opt/doc/builds-on new file mode 100644 index 00000000..dcafa033 --- /dev/null +++ b/SAND_top_opt/doc/builds-on @@ -0,0 +1 @@ +step-79 diff --git a/SAND_top_opt/doc/dependencies b/SAND_top_opt/doc/dependencies new file mode 100644 index 00000000..4df6ca25 --- /dev/null +++ b/SAND_top_opt/doc/dependencies @@ -0,0 +1 @@ +DEAL_II_WITH_CXX14 DEAL_II_WITH_TRILINOS \ No newline at end of file diff --git a/SAND_top_opt/doc/entry-name b/SAND_top_opt/doc/entry-name new file mode 100644 index 00000000..aa279b63 --- /dev/null +++ b/SAND_top_opt/doc/entry-name @@ -0,0 +1 @@ +Simultaneous Analysis and Design in Topology Optimization diff --git a/SAND_top_opt/doc/tooltip b/SAND_top_opt/doc/tooltip new file mode 100644 index 00000000..4f912101 --- /dev/null +++ b/SAND_top_opt/doc/tooltip @@ -0,0 +1 @@ +Simultaneous Analysis and Design in Topology Optimization using block preconditioning and matrix-free geometric multigrid. \ No newline at end of file diff --git a/SAND_top_opt/include/density_filter.h b/SAND_top_opt/include/density_filter.h new file mode 100644 index 00000000..51daba29 --- /dev/null +++ b/SAND_top_opt/include/density_filter.h @@ -0,0 +1,87 @@ +// +// Created by justin on 5/13/21. +// + +#ifndef SAND_DENSITY_FILTER_H +#define SAND_DENSITY_FILTER_H + + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +/* This object is designed to calculate and form the matrix corresponding to a convolution of the unfiltered density. + * Once formed, we have F*\sigma = \rho*/ +namespace SAND { + using namespace dealii; + namespace LA + { + using namespace dealii::LinearAlgebraTrilinos; + } + template + class DensityFilter { + public: + + MPI_Comm mpi_communicator; + std::vector owned_partitioning; + std::vector relevant_partitioning; + + DensityFilter(); + DynamicSparsityPattern filter_dsp; + LA::MPI::SparseMatrix filter_matrix; + LA::MPI::SparseMatrix filter_matrix_transpose; + SparsityPattern filter_sparsity_pattern; + void initialize(DoFHandler &dof_handler); + std::set find_relevant_neighbors(types::global_dof_index cell_index) const; + + private: + std::vector cell_m; + std::vector x_coord; + std::vector y_coord; + std::vector z_coord; + std::vector cell_m_part; + std::vector x_coord_part; + std::vector y_coord_part; + std::vector z_coord_part; + + ConditionalOStream pcout; + + }; +} +#endif //SAND_DENSITY_FILTER_H diff --git a/SAND_top_opt/include/input_information.h b/SAND_top_opt/include/input_information.h new file mode 100644 index 00000000..fc2ae7bd --- /dev/null +++ b/SAND_top_opt/include/input_information.h @@ -0,0 +1,61 @@ +// +// Created by justin on 3/11/21. +// + +#ifndef SAND_INPUT_INFORMATION_H +#define SAND_INPUT_INFORMATION_H +#include +#include "parameters_and_components.h" + +namespace SAND { + using namespace dealii; + + namespace Input + { + constexpr double volume_percentage = .5; + + //geometry options + constexpr unsigned int geometry_base = GeometryOptions::mbb; + + constexpr unsigned int dim = 3; + // constexpr unsigned int refinements = 3; + extern unsigned int refinements = 3; + + //nonlinear algorithm options + constexpr double initial_barrier_size = 25; + constexpr double min_barrier_size = 1e-12; + + constexpr double fraction_to_boundary = .7; + constexpr unsigned int max_steps=100; + + constexpr unsigned int barrier_reduction=BarrierOptions::loqo; + constexpr double required_norm = .0001; + + //density filter options + constexpr double filter_r = .25; + + //other options + constexpr double density_penalty_exponent = 3; + + //output options + constexpr bool output_full_preconditioned_matrix = false; + constexpr bool output_full_matrix = false; + constexpr bool output_parts_of_matrix = false; + + //Linear solver options + constexpr unsigned int solver_choice = SolverOptions::inexact_K_with_inexact_A_gmres; + constexpr bool use_eisenstat_walker = false; + constexpr double default_gmres_tolerance = 1e-6; + + extern unsigned int a_inv_iterations = 5; + extern unsigned int k_inv_iterations = 30; + + constexpr double a_rel_tol = 1e-12; + constexpr double k_rel_tol = 1e-12; + + //Material Options + constexpr double material_lambda = 1.; + constexpr double material_mu = 1.; + } +} +#endif //SAND_INPUT_INFORMATION_H \ No newline at end of file diff --git a/SAND_top_opt/include/kkt_system.h b/SAND_top_opt/include/kkt_system.h new file mode 100644 index 00000000..bb264a11 --- /dev/null +++ b/SAND_top_opt/include/kkt_system.h @@ -0,0 +1,210 @@ +// +// Created by justin on 2/17/21. +// + +#ifndef SAND_KKT_SYSTEM_H +#define SAND_KKT_SYSTEM_H +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +#include "../include/schur_preconditioner.h" +#include "../include/density_filter.h" +#include "matrix_free_elasticity.h" + +#include + +#include +#include +#include +namespace SAND { + + + + namespace LA + { + using namespace dealii::LinearAlgebraTrilinos; + } + using namespace dealii; + + template + class KktSystem { + + using SystemMFMatrixType = MF_Elasticity_Operator; + using LevelMFMatrixType = MF_Elasticity_Operator; + public: + MPI_Comm mpi_communicator; + std::vector owned_partitioning; + std::vector relevant_partitioning; + + KktSystem(); + + void + create_triangulation(); + + void + setup_boundary_values(); + + void + setup_filter_matrix(); + + void + setup_block_system(); + + void + assemble_block_system(const LA::MPI::BlockVector &state, const double barrier_size); + + LA::MPI::BlockVector + solve(const LA::MPI::BlockVector &state); + + LA::MPI::BlockVector + get_initial_state(); + + double + calculate_objective_value(const LA::MPI::BlockVector &state) const; + + double + calculate_barrier_distance(const LA::MPI::BlockVector &state) const; + + double + calculate_feasibility(const LA::MPI::BlockVector &state, const double barrier_size) const; + + double + calculate_convergence(const LA::MPI::BlockVector &state) const; + + void + output(const LA::MPI::BlockVector &state, const unsigned int j) const; + + void + calculate_initial_rhs_error(); + + double + calculate_rhs_norm(const LA::MPI::BlockVector &state, const double barrier_size) const; + + void + output_stl(const LA::MPI::BlockVector &state); + + private: + + LA::MPI::BlockVector + calculate_rhs(const LA::MPI::BlockVector &test_solution, const double barrier_size) const; + + BlockDynamicSparsityPattern dsp; + BlockSparsityPattern sparsity_pattern; + mutable LA::MPI::BlockSparseMatrix system_matrix; + mutable LA::MPI::BlockVector locally_relevant_solution; + mutable LA::MPI::BlockVector distributed_solution; + LA::MPI::BlockVector system_rhs; + parallel::distributed::Triangulation triangulation; + DoFHandler dof_handler; + DoFHandler dof_handler_displacement; + DoFHandler dof_handler_density; + + std::map displacement_to_system_dof_index_map; + MGLevelObject> level_displacement_to_system_dof_index_map; + + AffineConstraints constraints; + AffineConstraints displacement_constraints; + FESystem fe_nine; + FESystem fe_ten; + hp::FECollection fe_collection; + FESystem fe_displacement; + FE_DGQ fe_density; + const double density_ratio; + const double density_penalty_exponent; + + mutable DensityFilter density_filter; + + std::map boundary_values; + MGLevelObject> level_boundary_values; + ConditionalOStream pcout; + + double initial_rhs_error; + + MGConstrainedDoFs mg_constrained_dofs; + SystemMFMatrixType elasticity_matrix_mf; + MGLevelObject mg_matrices; + + OperatorCellData active_cell_data; + MGLevelObject> level_cell_data; + dealii::LinearAlgebra::distributed::Vector active_density_vector; + dealii::LinearAlgebra::distributed::Vector relevant_density_vector; + MGLevelObject> level_density_vector; + + MGTransferMatrixFree transfer; + MGTransferMatrixFree mg_transfer; + + using SmootherType = PreconditionChebyshev>; + mg::SmootherRelaxation> mg_smoother; + MGLevelObject smoother_data; + MGCoarseGridApplySmoother> mg_coarse; + MGLevelObject> level_dirichlet_boundary_dofs; + MGLevelObject> mg_level_constraints; + MGLevelObject> mg_interface_matrices; + + std::set dirichlet_boundary; + + LinearAlgebra::distributed::Vector distributed_displacement_sol; + LinearAlgebra::distributed::Vector distributed_displacement_rhs; + + }; +} + +#endif //SAND_KKT_SYSTEM_H diff --git a/SAND_top_opt/include/markov_filter.h b/SAND_top_opt/include/markov_filter.h new file mode 100644 index 00000000..8472fcee --- /dev/null +++ b/SAND_top_opt/include/markov_filter.h @@ -0,0 +1,30 @@ +// +// Created by justin on 2/17/21. +// + +#ifndef SAND_MARKOV_FILTER_H +#define SAND_MARKOV_FILTER_H +#include +#include + +using namespace dealii; + +class MarkovFilter +{ + +public: + void setup(const double objective_value, const double barrier_distance, const double feasibility, const double barrier_value); + void add_point(const double objective_value, const double barrier_distance, const double feasibility); + void update_barrier_value(const double barrier_value); + bool check_filter(const double objective_value, const double barrier_distance,const double feasibility) const; + +private: + double objective_value; + double barrier_distance; + double feasibility; + double barrier_value; + double filter_barrier_function_value; +}; + + +#endif //SAND_MARKOV_FILTER_H diff --git a/SAND_top_opt/include/matrix_free_elasticity.h b/SAND_top_opt/include/matrix_free_elasticity.h new file mode 100644 index 00000000..525136f9 --- /dev/null +++ b/SAND_top_opt/include/matrix_free_elasticity.h @@ -0,0 +1,108 @@ +#ifndef MATRIX_FREE_GMG_H +#define MATRIX_FREE_GMG_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../include/parameters_and_components.h" +#include "../include/input_information.h" + + + +namespace SAND +{ +using GMGNumberType = double; +using namespace dealii; + +template +struct OperatorCellData +{ + Table<2, VectorizedArray> density; + std::size_t + memory_consumption() const; +}; + +template +class MF_Elasticity_Operator + : public MatrixFreeOperators:: + Base> +{ +public: + + MPI_Comm mpi_communicator; + + using value_type = number; + + MF_Elasticity_Operator(); + + void set_cell_data (const OperatorCellData &data); + + void compute_diagonal() override; + + void clear() override; + +private: + + void apply_add (dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src) const override; + + + void local_apply (const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const; + + void local_apply_face (const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src, + const std::pair &face_range) const; + + + void local_apply_boundary_face (const dealii::MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src, + const std::pair &face_range) const; + + void local_compute_diagonal (const MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const unsigned int &, + const std::pair &cell_range) const; + + + const OperatorCellData *cell_data; + + ConditionalOStream pcout; + +}; + + +} + + + + +#endif // MATRIX_FREE_GMG_H diff --git a/SAND_top_opt/include/parameters_and_components.h b/SAND_top_opt/include/parameters_and_components.h new file mode 100644 index 00000000..978cd96c --- /dev/null +++ b/SAND_top_opt/include/parameters_and_components.h @@ -0,0 +1,68 @@ +// +// Created by justin on 3/11/21. +// + +#ifndef SAND_PARAMETERS_AND_COMPONENTS_H +#define SAND_PARAMETERS_AND_COMPONENTS_H +#include + + +namespace SAND { + using namespace dealii; + + namespace SolutionComponents { + template static constexpr unsigned int density_lower_slack_multiplier = 0; + template static constexpr unsigned int density_upper_slack_multiplier = 1; + template static constexpr unsigned int density_lower_slack = 2; + template static constexpr unsigned int density_upper_slack = 3; + template static constexpr unsigned int unfiltered_density = 4; + template static constexpr unsigned int displacement = 5; + template static constexpr unsigned int displacement_multiplier = 5 + dim; + template static constexpr unsigned int unfiltered_density_multiplier = 5 + 2 * dim; + template static constexpr unsigned int density = 6 + 2 * dim; + template static constexpr unsigned int total_volume_multiplier = 7 + 2 * dim; + } + + namespace SolutionBlocks { + static constexpr unsigned int density_lower_slack_multiplier = 0; + static constexpr unsigned int density_upper_slack_multiplier = 1; + static constexpr unsigned int density_lower_slack = 2; + static constexpr unsigned int density_upper_slack = 3; + static constexpr unsigned int unfiltered_density = 4; + static constexpr unsigned int displacement = 5; + static constexpr unsigned int displacement_multiplier = 6; + static constexpr unsigned int unfiltered_density_multiplier = 7; + static constexpr unsigned int density = 8; + static constexpr unsigned int total_volume_multiplier = 9; + } + + namespace BoundaryIds { + static constexpr types::boundary_id no_force = 101; + static constexpr types::boundary_id down_force = 102; + static constexpr types::boundary_id held_still = 103; + } + + namespace MaterialIds { + static constexpr types::material_id with_multiplier = 10; + static constexpr types::material_id without_multiplier = 9; + } + + namespace SolverOptions { + static constexpr unsigned int direct_solve = 1; + static constexpr unsigned int exact_preconditioner_with_gmres = 2; + static constexpr unsigned int inexact_K_with_exact_A_gmres = 3; + static constexpr unsigned int inexact_K_with_inexact_A_gmres = 4; + } + + namespace BarrierOptions { + static constexpr unsigned int loqo = 1; + static constexpr unsigned int monotone = 2; + static constexpr unsigned int mixed = 3; + } + namespace GeometryOptions { + static constexpr unsigned int mbb = 1; + static constexpr unsigned int l_shape = 2; + } + static constexpr unsigned int block_number = 10; +} +#endif //SAND_PARAMETERS_AND_COMPONENTS_H diff --git a/SAND_top_opt/include/sand_tools.h b/SAND_top_opt/include/sand_tools.h new file mode 100644 index 00000000..4d196647 --- /dev/null +++ b/SAND_top_opt/include/sand_tools.h @@ -0,0 +1,145 @@ +// +// Created by justin on 7/3/21. +// + +#ifndef SAND_MY_TOOLS_H +#define SAND_MY_TOOLS_H +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + + + +namespace SAND{ + namespace LA + { + using namespace dealii::LinearAlgebraTrilinos; + } + + using namespace dealii; + + + void build_matrix_element_by_element (const LinearOperator &op_X, + FullMatrix &X_matrix, + LA::MPI::Vector &exemplar_vector) + { + Threads::TaskGroup tasks; + for (unsigned int j=0; j &X_matrix) + { + const unsigned int n = X_matrix.n(); + const unsigned int m = X_matrix.m(); + std::ofstream Xmat(name); + for (unsigned int i = 0; i < m; i++) + { + Xmat << X_matrix(i, 0); + for (unsigned int j = 1; j < n; j++) + { + Xmat << "," << X_matrix(i, j); + } + Xmat << "\n"; + } + Xmat.close(); + } + + void print_matrix (std::string &name, SparseMatrix &X_matrix) + { + const unsigned int n = X_matrix.n(); + const unsigned int m = X_matrix.m(); + std::ofstream Xmat(name); + for (unsigned int i = 0; i < m; i++) + { + Xmat << X_matrix.el(i, 0); + for (unsigned int j = 1; j < n; j++) + { + Xmat << "," << X_matrix.el(i, j); + } + Xmat << "\n"; + } + Xmat.close(); + } + + void print_matrix (std::string &name, LA::MPI::SparseMatrix &X_matrix) + { + const unsigned int n = X_matrix.n(); + const unsigned int m = X_matrix.m(); + std::ofstream Xmat(name); + for (unsigned int i = 0; i < m; i++) + { + Xmat << X_matrix.el(i, 0); + for (unsigned int j = 1; j < n; j++) + { + Xmat << "," << X_matrix.el(i, j); + } + Xmat << "\n"; + } + Xmat.close(); + } + +} + + + + +#endif //SAND_MY_TOOLS_H diff --git a/SAND_top_opt/include/schur_preconditioner.h b/SAND_top_opt/include/schur_preconditioner.h new file mode 100644 index 00000000..8f06feb7 --- /dev/null +++ b/SAND_top_opt/include/schur_preconditioner.h @@ -0,0 +1,470 @@ +// +// Created by justin on 3/2/21. +// + +#ifndef SAND_SCHUR_PRECONDITIONER_H +#define SAND_SCHUR_PRECONDITIONER_H +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + + +#include +#include +#include + + +#include "../include/parameters_and_components.h" +#include "../include/poly_pre.h" +#include "matrix_free_elasticity.h" + +#include + +#include +#include + + +namespace SAND +{ + using MatrixType = dealii::TrilinosWrappers::SparseMatrix; + using VectorType = dealii::TrilinosWrappers::MPI::Vector; + using PayloadType = dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload; + using PayloadVectorType = typename PayloadType::VectorType; + namespace LA + { + using namespace dealii::LinearAlgebraTrilinos; + } + using namespace dealii; + + class VmultTrilinosSolverDirect : public TrilinosWrappers::SparseMatrix { + public: + VmultTrilinosSolverDirect(SolverControl &cn, + const TrilinosWrappers::SolverDirect::AdditionalData &data + ); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void vmult(LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const; + void initialize(LA::MPI::SparseMatrix &a_mat); + unsigned int m() const; + unsigned int n() const; + int get_size() + { + return size; + } + private: + mutable TrilinosWrappers::SolverDirect solver_direct; + int size; + }; + + class AMatWrapped : public TrilinosWrappers::SparseMatrix { + public: + AMatWrapped(LA::MPI::SparseMatrix &a_mat); + void vmult(LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const; + void Tvmult(LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const; + unsigned int m() const; + unsigned int n() const; + void set_exemplar_vector (const LA::MPI::Vector &exemplar_vector) + { + temp_dst = exemplar_vector; + temp_src = exemplar_vector; + } + + private: + const LA::MPI::SparseMatrix &a_mat; + mutable LA::MPI::Vector temp_src; + mutable LA::MPI::Vector temp_dst; + + }; + + template + class AInvMatMFGMG : public TrilinosWrappers::SparseMatrix { + public: + AInvMatMFGMG(MF_Elasticity_Operator &mf_elasticity_operator_in , PreconditionMG,MGTransferMatrixFree> + &mf_gmg_preconditioner_in,LA::MPI::SparseMatrix &a_mat, std::map &displacement_to_system_dof_index_map); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + unsigned int m() const; + unsigned int n() const; + void set_tol(double tolerance_in); + void set_iter(unsigned int iterations_in); + void set_exemplar_vector (const LA::MPI::Vector &exemplar_vector) + { + a_mat_wrapped.set_exemplar_vector(exemplar_vector); + } + double tolerance = Input::a_rel_tol; + unsigned int iterations = Input::a_inv_iterations; + private: + AMatWrapped a_mat_wrapped; + MF_Elasticity_Operator &mf_elasticity_operator; + PreconditionMG,MGTransferMatrixFree> &mf_gmg_preconditioner; + const std::map displacement_to_system_dof_index_map; + + mutable dealii::LinearAlgebra::distributed::Vector temp_src; + mutable dealii::LinearAlgebra::distributed::Vector temp_dst; + + }; + + + class GMatrix : public TrilinosWrappers::SparseMatrix { + public: + GMatrix(const LA::MPI::SparseMatrix &f_mat_in,const LA::MPI::SparseMatrix &f_t_mat_in, LA::MPI::SparseMatrix &d_8_mat_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void initialize(LA::MPI::Vector &exemplar_density_vector, LA::MPI::SparseMatrix &d_m_inv_mat_in); + unsigned int m() const; + unsigned int n() const; + private: + const LA::MPI::SparseMatrix &f_mat; + const LA::MPI::SparseMatrix &f_t_mat; + LA::MPI::SparseMatrix &d_8_mat; + LA::MPI::SparseMatrix d_m_inv_mat; + mutable LA::MPI::Vector temp_vect_1; + mutable LA::MPI::Vector temp_vect_2; + + + }; + + template + class HMatrix : public TrilinosWrappers::SparseMatrix { + public: + HMatrix(LA::MPI::SparseMatrix &a_mat_in, const LA::MPI::SparseMatrix &b_mat_in, const LA::MPI::SparseMatrix &c_mat_in, const LA::MPI::SparseMatrix &e_mat_in,TrilinosWrappers::PreconditionAMG &pre_amg_in, VmultTrilinosSolverDirect &a_inv_direct_in, AInvMatMFGMG &a_inv_mf_gmg_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void initialize(LA::MPI::Vector &exemplar_density_vector, LA::MPI::Vector &exemplar_displacement_vector, LA::MPI::SparseMatrix &d_m_inv_mat_in); + unsigned int m() const; + unsigned int n() const; + private: + LA::MPI::SparseMatrix &a_mat; + const LA::MPI::SparseMatrix &b_mat; + const LA::MPI::SparseMatrix &c_mat; + const LA::MPI::SparseMatrix &e_mat; + LA::MPI::SparseMatrix d_m_inv_mat; + TrilinosWrappers::PreconditionAMG &pre_amg; + VmultTrilinosSolverDirect &a_inv_direct; + mutable LA::MPI::Vector temp_vect_1; + mutable LA::MPI::Vector temp_vect_2; + mutable LA::MPI::Vector temp_vect_3; + mutable LA::MPI::Vector temp_vect_4; + mutable LA::MPI::Vector temp_vect_5; + mutable LA::MPI::Vector temp_vect_6; + mutable LA::MPI::Vector temp_vect_7; + AInvMatMFGMG &a_inv_mf_gmg; + + }; + + template + class HMatrixDirect : public TrilinosWrappers::SparseMatrix { + public: + HMatrixDirect(LA::MPI::SparseMatrix &a_mat_in, const LA::MPI::SparseMatrix &b_mat_in, const LA::MPI::SparseMatrix &c_mat_in, const LA::MPI::SparseMatrix &e_mat_in,TrilinosWrappers::PreconditionAMG &pre_amg_in, VmultTrilinosSolverDirect &a_inv_direct_in, AInvMatMFGMG &a_inv_mf_gmg_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void initialize(LA::MPI::Vector &exemplar_density_vector, LA::MPI::Vector &exemplar_displacement_vector); + unsigned int m() const; + unsigned int n() const; + private: + LA::MPI::SparseMatrix &a_mat; + const LA::MPI::SparseMatrix &b_mat; + const LA::MPI::SparseMatrix &c_mat; + const LA::MPI::SparseMatrix &e_mat; + TrilinosWrappers::PreconditionAMG &pre_amg; + VmultTrilinosSolverDirect &a_inv_direct; + mutable LA::MPI::Vector temp_vect_1; + mutable LA::MPI::Vector temp_vect_2; + mutable LA::MPI::Vector temp_vect_3; + mutable LA::MPI::Vector temp_vect_4; + mutable LA::MPI::Vector temp_vect_5; + mutable LA::MPI::Vector temp_vect_6; + mutable LA::MPI::Vector temp_vect_7; + AInvMatMFGMG &a_inv_mf_gmg; + + }; + + template + class KinvMatrix : public TrilinosWrappers::SparseMatrix { + public: + KinvMatrix(HMatrix &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void initialize(LA::MPI::Vector &exemplar_density_vector, LA::MPI::SparseMatrix &d_m_inv_mat_in); + unsigned int m() const; + unsigned int n() const; + private: + HMatrix &h_mat; + GMatrix &g_mat; + const LA::MPI::SparseMatrix &d_m_mat; + LA::MPI::SparseMatrix d_m_inv_mat; + mutable LA::MPI::Vector temp_vect_1; + mutable LA::MPI::Vector temp_vect_2; + mutable LA::MPI::Vector temp_vect_3; + mutable LA::MPI::Vector temp_vect_4; + }; + + template + class KinvMatrixPart : public Subscriptor { + public: + KinvMatrixPart(HMatrix &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void initialize(LA::MPI::Vector &exemplar_density_vector); + unsigned int m() const; + unsigned int n() const; + private: + HMatrix &h_mat; + GMatrix &g_mat; + const LA::MPI::SparseMatrix &d_m_mat; + LA::MPI::SparseMatrix d_m_inv_mat; + mutable LA::MPI::Vector temp_vect_1; + mutable LA::MPI::Vector temp_vect_2; + mutable LA::MPI::Vector temp_vect_3; + mutable LA::MPI::Vector temp_vect_4; + }; + + template + class KinvMatrixDirect : public TrilinosWrappers::SparseMatrix { + public: + KinvMatrixDirect(HMatrixDirect &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void initialize(LA::MPI::Vector &exemplar_density_vector); + unsigned int m() const; + unsigned int n() const; + private: + HMatrixDirect &h_mat; + GMatrix &g_mat; + const LA::MPI::SparseMatrix &d_m_mat; + LA::MPI::SparseMatrix d_m_inv_mat; + mutable LA::MPI::Vector temp_vect_1; + mutable LA::MPI::Vector temp_vect_2; + mutable LA::MPI::Vector temp_vect_3; + mutable LA::MPI::Vector temp_vect_4; + }; + + template + class JinvMatrix : public TrilinosWrappers::SparseMatrix { + public: + JinvMatrix(HMatrix &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void initialize(LA::MPI::Vector &exemplar_density_vector, LA::MPI::SparseMatrix &d_m_inv_mat_in); + unsigned int m() const; + unsigned int n() const; + private: + HMatrix &h_mat; + GMatrix &g_mat; + const LA::MPI::SparseMatrix &d_m_mat; + LA::MPI::SparseMatrix d_m_inv_mat; + mutable LA::MPI::Vector temp_vect_1; + mutable LA::MPI::Vector temp_vect_2; + mutable LA::MPI::Vector temp_vect_3; + mutable LA::MPI::Vector temp_vect_4; + }; + + template + class JinvMatrixPart : public Subscriptor { + public: + JinvMatrixPart(HMatrix &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void initialize(LA::MPI::Vector &exemplar_density_vector); + unsigned int m() const; + unsigned int n() const; + private: + HMatrix &h_mat; + GMatrix &g_mat; + const LA::MPI::SparseMatrix &d_m_mat; + LA::MPI::SparseMatrix d_m_inv_mat; + mutable LA::MPI::Vector temp_vect_1; + mutable LA::MPI::Vector temp_vect_2; + mutable LA::MPI::Vector temp_vect_3; + mutable LA::MPI::Vector temp_vect_4; + }; + + template + class JinvMatrixDirect : public TrilinosWrappers::SparseMatrix { + public: + JinvMatrixDirect(HMatrixDirect &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + void initialize(LA::MPI::Vector &exemplar_density_vector); + unsigned int m() const; + unsigned int n() const; + private: + HMatrixDirect &h_mat; + GMatrix &g_mat; + const LA::MPI::SparseMatrix &d_m_mat; + LA::MPI::SparseMatrix d_m_inv_mat; + mutable LA::MPI::Vector temp_vect_1; + mutable LA::MPI::Vector temp_vect_2; + mutable LA::MPI::Vector temp_vect_3; + mutable LA::MPI::Vector temp_vect_4; + }; + + + + template + class TopOptSchurPreconditioner: public Subscriptor { + public: + TopOptSchurPreconditioner(LA::MPI::BlockSparseMatrix &matrix_in, DoFHandler &big_dof_handler_in, MF_Elasticity_Operator &mf_elasticity_operator_in , PreconditionMG,MGTransferMatrixFree> + &mf_gmg_preconditioner_in, std::map &displacement_to_system_dof_index_map); + void initialize (LA::MPI::BlockSparseMatrix &matrix, const std::map &boundary_values, const DoFHandler &dof_handler, const LA::MPI::BlockVector &distributed_state); + void vmult(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const; + void Tvmult(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const; + void vmult_add(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const; + void Tvmult_add(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const; + void clear(); + unsigned int m() const; + unsigned int n() const; + void get_sparsity_pattern(BlockDynamicSparsityPattern &bdsp); + + void assemble_mass_matrix(const LA::MPI::BlockVector &state, const hp::FECollection &fe_system, const DoFHandler &dof_handler, const AffineConstraints &constraints, const BlockSparsityPattern &bsp); + + void print_stuff(); + + LA::MPI::BlockSparseMatrix &system_matrix; + + private: + MPI_Comm mpi_communicator; + unsigned int n_rows; + unsigned int n_columns; + unsigned int n_block_rows; + unsigned int n_block_columns; + void vmult_step_1(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const; + void vmult_step_2(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const; + void vmult_step_3(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const; + void vmult_step_4(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const; + void vmult_step_5(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const; + + BlockSparsityPattern mass_sparsity; + LA::MPI::BlockSparseMatrix approx_h_mat; + + SolverControl other_solver_control; + mutable SolverBicgstab other_bicgstab; + mutable SolverGMRES other_gmres; + mutable SolverCG other_cg; + + LA::MPI::SparseMatrix &a_mat; + const LA::MPI::SparseMatrix &b_mat; + const LA::MPI::SparseMatrix &c_mat; + const LA::MPI::SparseMatrix &e_mat; + const LA::MPI::SparseMatrix &f_mat; + const LA::MPI::SparseMatrix &f_t_mat; + const LA::MPI::SparseMatrix &d_m_mat; + const LA::MPI::SparseMatrix &d_1_mat; + const LA::MPI::SparseMatrix &d_2_mat; + const LA::MPI::SparseMatrix &m_vect; + const DoFHandler &big_dof_handler; + + LA::MPI::SparseMatrix d_3_mat; + LA::MPI::SparseMatrix d_4_mat; + LA::MPI::SparseMatrix d_5_mat; + LA::MPI::SparseMatrix d_6_mat; + LA::MPI::SparseMatrix d_7_mat; + LA::MPI::SparseMatrix d_8_mat; + LA::MPI::SparseMatrix d_m_inv_mat; + + mutable LA::MPI::Vector pre_j; + mutable LA::MPI::Vector pre_k; + mutable LA::MPI::Vector g_d_m_inv_density; + mutable LA::MPI::Vector k_g_d_m_inv_density; + + std::string solver_type; + TrilinosWrappers::SolverDirect::AdditionalData additional_data; + SolverControl direct_solver_control; + mutable VmultTrilinosSolverDirect a_inv_direct; + + mutable AInvMatMFGMG a_inv_mf_gmg; + ConditionalOStream pcout; + mutable TimerOutput timer; + + mutable TrilinosWrappers::PreconditionAMG pre_amg; + + mutable GMatrix g_mat; + HMatrix h_mat; + + JinvMatrix j_inv_mat; + KinvMatrix k_inv_mat; + + JinvMatrixPart j_inv_part; + KinvMatrixPart k_inv_part; + + mutable int num_mults; + + }; + + template + class PolyPreJ { + + public: + PolyPreJ(const JinvMatrixPart &inner_matrix_in, const int degree_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + // void vmult(LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const; + + + private: + const JinvMatrixPart &inner_matrix; + const int degree; + }; + + template + class PolyPreK { + + public: + PolyPreK(const KinvMatrixPart &inner_matrix_in, const int degree_in); + void vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const; + // void vmult(LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const; + + + private: + const KinvMatrixPart &inner_matrix; + const int degree; + }; + + +} +#endif //SAND_SCHUR_PRECONDITIONER_H diff --git a/SAND_top_opt/include/watchdog.h b/SAND_top_opt/include/watchdog.h new file mode 100644 index 00000000..3ce4dceb --- /dev/null +++ b/SAND_top_opt/include/watchdog.h @@ -0,0 +1,69 @@ +#include +#include +#include +#include +#include +#include "../include/markov_filter.h" +#include "../include/kkt_system.h" +#include "../include/input_information.h" +#include +#include +#include +#include +#include + +///Above are fairly normal files to include. I also use the sparse direct package, which requiresBLAS/LAPACK +/// to perform a direct solve while I work on a fast iterative solver for this problem. + +namespace SAND { + namespace LA + { + using namespace dealii::LinearAlgebraTrilinos; + } + + + using namespace dealii; + + /// Below is the main class for solving this problem. It handles the nonlinear solver portion of the problem, + /// taking information from the KKTSystem class for step directions, and calculating step lengths. This class + /// not only takes those steps, but handles the barrier parameter for the log barrier used. + template + class NonlinearWatchdog { + public: + NonlinearWatchdog(); + + void + run(); + + private: + MPI_Comm mpi_communicator; + std::pair + calculate_max_step_size(const LA::MPI::BlockVector &state, const LA::MPI::BlockVector &step) const; + + const LA::MPI::BlockVector + find_max_step(const LA::MPI::BlockVector &state); + + LA::MPI::BlockVector + take_scaled_step(const LA::MPI::BlockVector &state,const LA::MPI::BlockVector &max_step) const; + + bool + check_convergence(const LA::MPI::BlockVector &state) const; + + void + update_barrier(LA::MPI::BlockVector ¤t_state); + + void + perform_initial_setup(); + + void + nonlinear_step(LA::MPI::BlockVector ¤t_state, LA::MPI::BlockVector ¤t_step, const unsigned int max_uphill_steps, unsigned int &iteration_number); + + KktSystem kkt_system; + MarkovFilter markov_filter; + double barrier_size; + bool mixed_barrier_monotone_mode; + ConditionalOStream pcout; + TimerOutput overall_timer; + }; + +} // namespace SAND diff --git a/SAND_top_opt/readme.md b/SAND_top_opt/readme.md new file mode 100644 index 00000000..0ab4f67c --- /dev/null +++ b/SAND_top_opt/readme.md @@ -0,0 +1,7 @@ +Topology optimization is a class of algorithms designed to optimize a design or structure to accomplish some goal. It is part of a process of computer generated design that allows engineers to design better products faster. + +Several topology optimization algorithms have been developed to generate these designs, and having several methods is itself a good thing. The minimization problems resulting from topology optimization have many local minima, and different algorithms often hit different local minima. Topology optimization is a tool used by engineers, and having more generated designs allows engineers to better make decisions with these designs. + +One such algorithm that has piqued the imagination of developers is called Simultaneous Analysis and Design (SAND), especially in the context of Interior Point Optimization (IPO). This method is known to generate extremely optimal designs, and is good at avoiding local minima. However, this method is not used in practice, due to its computational cost. + +This code is an attempt at an effective algorithm to generate a design using it. The most expensive part is a linear solve of a block system. A multi-tier preconditioning approach is used to solve this system in a reasonable amount of time. \ No newline at end of file diff --git a/SAND_top_opt/source/density_filter.cc b/SAND_top_opt/source/density_filter.cc new file mode 100644 index 00000000..4e64a56f --- /dev/null +++ b/SAND_top_opt/source/density_filter.cc @@ -0,0 +1,268 @@ +// +// Created by justin on 5/13/21. +// +#include "../include/density_filter.h" +#include "../include/input_information.h" +#include +#include +#include +#include +#include +#include +#include + +namespace SAND +{ + using namespace dealii; + + template + DensityFilter::DensityFilter() : + mpi_communicator(MPI_COMM_WORLD), + pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + { + } + + ///When initialized, this function takes the current triangulation and creates a matrix corresponding to a + /// convolution being applied to a piecewise constant function on that triangulation + /// + template + void + DensityFilter::initialize(DoFHandler &dof_handler) { + ///Start by making vectors to fill with information about the x,y,z coords of centers of cells + /// + std::vector block_component(10, 2); + block_component[SolutionBlocks::density] = 0; + block_component[SolutionBlocks::displacement] = 1; + const std::vector dofs_per_block = + DoFTools::count_dofs_per_fe_block(dof_handler, block_component); + const unsigned int n_p = dofs_per_block[0]; + IndexSet local_owned = dof_handler.locally_owned_dofs().get_view(0, n_p); + x_coord.resize(n_p); + y_coord.resize(n_p); + z_coord.resize(n_p); + auto row_sum = z_coord; + auto row_sum_full = z_coord; + cell_m.resize(n_p); + x_coord_part.resize(n_p); + y_coord_part.resize(n_p); + z_coord_part.resize(n_p); + cell_m_part.resize(n_p); + + filter_dsp.reinit(dofs_per_block[0], + dofs_per_block[0]); + filter_sparsity_pattern.copy_from(filter_dsp); + + // const auto owned_dofs = dof_handler.locally_owned_dofs().get_view(0, dofs_per_block[0]); + + // filter_matrix.reinit(owned_dofs, filter_sparsity_pattern, MPI_COMM_WORLD); + + + // //identity filter + // for (const auto &cell : dof_handler.active_cell_iterators()) + // { + // if(cell->is_locally_owned()) + // { + // std::vector i(cell->get_fe().n_dofs_per_cell()); + // cell->get_dof_indices(i); + // filter_matrix.add(i[cell->get_fe().component_to_system_index(0, 0)], i[cell->get_fe().component_to_system_index(0, 0)], 1.0); + // } + // } + + std::set neighbor_ids; + std::set::cell_iterator> cells_to_check; + std::set::cell_iterator> cells_to_check_temp; + ///finds neighbors whose values would be relevant, and adds them to the sparsity pattern of the matrix + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if(cell->is_locally_owned()) + { + std::vector i(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(i); + const unsigned int i_val = i[cell->get_fe().component_to_system_index(0, 0)]; + x_coord_part[i_val] = cell->center()[0] ; + y_coord_part[i_val] = cell->center()[1] ; + cell_m_part[i_val] = cell->measure(); + if (dim==3) + { + z_coord_part[i_val] = cell->center()[2] ; + } + } + } + MPI_Allreduce(x_coord_part.data(), x_coord.data(), n_p, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(y_coord_part.data(), y_coord.data(), n_p, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(z_coord_part.data(), z_coord.data(), n_p, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(cell_m_part.data(), cell_m.data(), n_p, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + for (const auto &cell : dof_handler.active_cell_iterators()) { + if(cell->is_locally_owned()) + { + std::vector i(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(i); + for (const auto &neighbor_cell_index : find_relevant_neighbors(i[cell->get_fe().component_to_system_index(0, 0)])) + { + filter_dsp.add(i[cell->get_fe().component_to_system_index(0, 0)], neighbor_cell_index); + } + } + } + filter_sparsity_pattern.copy_from(filter_dsp); + const auto owned_dofs = dof_handler.locally_owned_dofs().get_view(4 * n_p, 5 * n_p); + + filter_matrix.reinit(owned_dofs, filter_sparsity_pattern, MPI_COMM_WORLD); + filter_matrix_transpose.reinit(owned_dofs, filter_sparsity_pattern, MPI_COMM_WORLD); + + /// adds values to the matrix corresponding to the max radius - distance + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if(cell->is_locally_owned()) + { + std::vector i(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(i); + auto cell_index = i[cell->get_fe().component_to_system_index(0, 0)]; + double value_total = 0; + for (const auto &neighbor_cell_index : find_relevant_neighbors(cell_index)) + { + double d_x = std::abs(x_coord[cell_index]-x_coord[neighbor_cell_index]); + double d_y = std::abs(y_coord[cell_index]-y_coord[neighbor_cell_index]); + double d; + if (dim==3) + { + double d_z = std::abs(z_coord[cell_index]-z_coord[neighbor_cell_index]); + d = std::pow(d_x*d_x + d_y*d_y + d_z*d_z , .5); + } + else + { + d = std::pow(d_x*d_x + d_y*d_y , .5); + } + ///value should be (max radius - distance between cells)*cell measure + double value = (Input::filter_r - d)*cell_m[neighbor_cell_index]; + value_total += value; + } + row_sum[cell_index] = value_total; + + } + } + MPI_Allreduce(row_sum.data(), row_sum_full.data(), n_p, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if(cell->is_locally_owned()) + { + std::vector i(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(i); + auto cell_index = i[cell->get_fe().component_to_system_index(0, 0)]; + for (const auto &neighbor_cell_index : find_relevant_neighbors(cell_index)) + { + double d_x = std::abs(x_coord[cell_index]-x_coord[neighbor_cell_index]); + double d_y = std::abs(y_coord[cell_index]-y_coord[neighbor_cell_index]); + double d; + if (dim==3) + { + double d_z = std::abs(z_coord[cell_index]-z_coord[neighbor_cell_index]); + d = std::pow(d_x*d_x + d_y*d_y + d_z*d_z , .5); + } + else + { + d = std::pow(d_x*d_x + d_y*d_y , .5); + } + ///value should be (max radius - distance between cells)*cell measure + double value = (Input::filter_r - d)*cell_m[neighbor_cell_index]; + filter_matrix.set(cell_index, neighbor_cell_index, value/row_sum_full[cell_index]); + filter_matrix_transpose.set(cell_index, neighbor_cell_index, value/row_sum_full[neighbor_cell_index]); + } + } + } + double sum_sum = 0; + double x_sum = 0; + double y_sum = 0; + double z_sum = 0; + for (int i=0; i< row_sum_full.size(); i++) + { + sum_sum += std::abs(row_sum_full[i]); + x_sum += std::abs(x_coord[i]); + y_sum += std::abs(y_coord[i]); + z_sum += std::abs(z_coord[i]); + } + filter_matrix.compress(VectorOperation::insert); + filter_matrix_transpose.compress(VectorOperation::insert); + + ///here we normalize the filter so it computes an average. Sum of values in a row should be 1 + // for (const auto &cell : dof_handler.active_cell_iterators()) + // { + // if(cell->is_locally_owned()) + // { + // std::vector i(cell->get_fe().n_dofs_per_cell()); + // cell->get_dof_indices(i); + // const int i_ind = cell->get_fe().component_to_system_index(0, 0); + // double denominator = 0; + // typename LA::MPI::SparseMatrix::iterator iter = filter_matrix.begin( + // i[i_ind]); + // for (; iter != filter_matrix.end(i[i_ind]); iter++) + // { + // denominator = denominator + iter->value(); + // } + // iter = filter_matrix.begin(i[i_ind]); + // for (; iter != filter_matrix.end(i[i_ind]); iter++) + // { + // iter->value() = iter->value() / denominator; + // } + // } + // } + // pcout << "FILTER TO 5" << std::endl; + // LA::MPI::Vector test_density_start; + // test_density_start.reinit(local_owned, mpi_communicator); + // LA::MPI::Vector test_density_end; + // test_density_end.reinit(local_owned, mpi_communicator); + // test_density_end = 0.; + // for (int i = 0; i + std::set + DensityFilter::find_relevant_neighbors(types::global_dof_index cell_index) const + { + double d_x,d_y,d_z; + std::set relevant_cells; + for (unsigned int i=0; i < x_coord.size(); i++) + { + d_x = std::abs(x_coord[cell_index]-x_coord[i]); + + if (d_x < Input::filter_r) + { + d_y = std::abs(y_coord[cell_index]-y_coord[i]); + + if ((d_x*d_x + d_y*d_y) < (Input::filter_r*Input::filter_r)) + { + + if (dim == 3) + { + d_z = std::abs(z_coord[cell_index]-z_coord[i]); + + if ((d_x*d_x + d_y*d_y + d_z*d_z) < (Input::filter_r*Input::filter_r)) + { + relevant_cells.insert(i); + } + } + else + { + relevant_cells.insert(i); + } + + } + } + } + return relevant_cells; + + } + +}//SAND namespace + template class SAND::DensityFilter<2>; + template class SAND::DensityFilter<3>; diff --git a/SAND_top_opt/source/kkt_system.cc b/SAND_top_opt/source/kkt_system.cc new file mode 100644 index 00000000..12b5c8b9 --- /dev/null +++ b/SAND_top_opt/source/kkt_system.cc @@ -0,0 +1,3220 @@ +// +// Created by justin on 2/17/21. +// +#include "../include/kkt_system.h" +#include +#include +#include +#include + +#include +#include +#include +#include +#include + + +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + + +#include + +#include "../include/input_information.h" +#include "../include/matrix_free_elasticity.h" +#include "../include/poly_pre.h" + +#include +#include + +/// This problem initializes with a FESystem composed of 2×dim FE_Q(1) elements, and 8 FE_DGQ(0) elements. +/// The piecewise constant functions are for density-related variables,and displacement-related variables are assigned to the FE_Q(1) elements. +namespace SAND { + +///Necessary functions for going between Trilinos vectors and multigrid-compatible distributed vectors. +namespace ChangeVectorTypes +{ +template +void copy_from_displacement_to_system_vector(LA::MPI::Vector &out, + const dealii::LinearAlgebra::distributed::Vector &in, + std::map & displacement_to_system_dof_index_map) +{ +// dealii::LinearAlgebra::ReadWriteVector rwv( +// out.locally_owned_elements()); +// rwv.import(in, VectorOperation::insert); + for (const auto &index_pair : displacement_to_system_dof_index_map) + { + out[index_pair.second] = in[index_pair.first]; + } +// out.import(rwv, VectorOperation::insert); +} + +template +void copy_from_system_to_displacement_vector(dealii::LinearAlgebra::distributed::Vector &out, + const LA::MPI::Vector &in, + std::map & displacement_to_system_dof_index_map) +{ +// dealii::LinearAlgebra::ReadWriteVector rwv( +// out.locally_owned_elements()); +// rwv.import(in, VectorOperation::insert); + ConditionalOStream pcout (std::cout,(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 1)); + for (const auto &index_pair : displacement_to_system_dof_index_map) + { + out[index_pair.first] = in[index_pair.second]; + } +// out.import(rwv, VectorOperation::insert); +} + +template +void copy(LA::MPI::Vector & out, + const dealii::LinearAlgebra::distributed::Vector &in) +{ + dealii::LinearAlgebra::ReadWriteVector rwv( + out.locally_owned_elements()); + rwv.import(in, VectorOperation::insert); + out.import(rwv, VectorOperation::insert); +} +template +void copy(dealii::LinearAlgebra::distributed::Vector &out, + const LA::MPI::Vector & in) +{ + dealii::LinearAlgebra::ReadWriteVector rwv; + rwv.reinit(in); + out.import(rwv, VectorOperation::insert); +} + + +} // namespace ChangeVectorTypes + +///The KKTSystem class calculates the Hessian and Gradient of the Lagrangian of the system, and solves the resulting system to be used +/// as a step direction for the overarching solver. +template +KktSystem::KktSystem() + : + mpi_communicator(MPI_COMM_WORLD), + triangulation(mpi_communicator, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy), + dof_handler(triangulation), + dof_handler_displacement(triangulation), + dof_handler_density(triangulation), + /*fe should have 1 FE_DGQ(0) element for density, dim FE_Q finite elements for displacement, + * another dim FE_Q elements for the lagrange multiplier on the FE constraint, and 2 more FE_DGQ(0) + * elements for the upper and lower bound constraints */ + fe_nine(FE_DGQ(0) ^ 5, + (FESystem(FE_Q(1) ^ dim)) ^ 2, + FE_DGQ(0) ^ 2, + FE_Nothing() ^ 1), + fe_ten(FE_DGQ(0) ^ 5, + (FESystem(FE_Q(1) ^ dim)) ^ 2, + FE_DGQ(0) ^ 2, + FE_DGQ(0) ^ 1), + fe_displacement(FE_Q(1) ^ dim), + fe_density(0), + density_ratio(Input::volume_percentage), + density_penalty_exponent(Input::density_penalty_exponent), + density_filter(), + pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) +{ + fe_collection.push_back(fe_nine); + fe_collection.push_back(fe_ten); + +} + + +///A function used once at the beginning of the program, this creates a matrix H so that H* unfiltered density = filtered density + +template +void +KktSystem::setup_filter_matrix() { + pcout << "IN KKT FILTER SETUP FUNCTION" << std::endl; + density_filter.initialize(dof_handler); +} + +///This triangulation matches the problem description + +template +void +KktSystem::create_triangulation() { + + ///Start by defining the sub-blocks of the DoFHandler + + std::vector sub_blocks(2*dim+8, 0); + + sub_blocks[0]=0; + sub_blocks[1]=1; + sub_blocks[2]=2; + sub_blocks[3]=3; + sub_blocks[4]=4; + for(int i=0; i(0, 0), + Point(width, height)); + + triangulation.refine_global(Input::refinements); + + /*Set BCIDs */ + for (const auto &cell: dof_handler.active_cell_iterators()) + { + if(cell->is_locally_owned()) + { + cell->set_active_fe_index(0); + cell->set_material_id(MaterialIds::without_multiplier); + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) + { + if (cell->face(face_number)->at_boundary()) + { + const auto center = cell->face(face_number)->center(); + + if (std::fabs(center(1) - downforce_y) < 1e-12) + { + if (std::fabs(center(0) - downforce_x) < downforce_size) + { + cell->face(face_number)->set_boundary_id(BoundaryIds::down_force); + } + else + { + cell->face(face_number)->set_boundary_id(BoundaryIds::no_force); + } + } + } + } + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) { + if (std::abs(cell->vertex(vertex_number)(0)) + std::abs(cell->vertex(vertex_number)(1)) < + 1e-10) { + cell->set_active_fe_index(1); + cell->set_material_id(MaterialIds::with_multiplier); + } + } + } + } + + } + else if (dim == 3) + { + GridGenerator::subdivided_hyper_rectangle(triangulation, + {width_refine, height_refine, depth_refine}, + Point(0, 0, 0), + Point(width, height, depth)); + + triangulation.refine_global(Input::refinements); + + for (const auto &cell: dof_handler.active_cell_iterators()) { + if (cell->is_locally_owned()) + { + cell->set_active_fe_index(0); + cell->set_material_id(MaterialIds::without_multiplier); + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) { + if (cell->face(face_number)->at_boundary()) { + const auto center = cell->face(face_number)->center(); + + if (std::fabs(center(1) - downforce_y) < 1e-12) { + if (std::fabs(center(0) - downforce_x) < downforce_size) { + cell->face(face_number)->set_boundary_id(BoundaryIds::down_force); + } else { + cell->face(face_number)->set_boundary_id(BoundaryIds::no_force); + } + } + } + } + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) { + if (std::abs(cell->vertex(vertex_number)(0)) + std::abs(cell->vertex(vertex_number)(1)) + + std::abs(cell->vertex(vertex_number)(2)) < 1e-10) { + cell->set_active_fe_index(1); + cell->set_material_id(MaterialIds::with_multiplier); + } + } + } + } + + } else { + throw; + } + ///L-shaped cantilever with re-entrant corner + } else if (Input::geometry_base == GeometryOptions::l_shape) { + const double width = 2; + const unsigned int width_refine = 2; + const double height = 2; + const unsigned int height_refine = 2; + const double depth = 1; + const unsigned int depth_refine = 1; + const double downforce_x = 2; + const double downforce_y = 1; + const double downforce_z = .5; + const double downforce_size = .3; + + if (dim == 2) { + GridGenerator::subdivided_hyper_L(triangulation, + {width_refine, height_refine}, + Point(0, 0), + Point(width, height), + {-1, -1}); + + triangulation.refine_global(Input::refinements); + + /*Set BCIDs */ + for (const auto &cell: dof_handler.active_cell_iterators()) { + if (cell->is_locally_owned()) + { + cell->set_active_fe_index(0); + cell->set_material_id(MaterialIds::without_multiplier); + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) { + if (cell->face(face_number)->at_boundary()) { + const auto center = cell->face(face_number)->center(); + + if (std::fabs(center(0) - downforce_x) < 1e-12) { + if (std::fabs(center(1) - downforce_y) < downforce_size) { + cell->face(face_number)->set_boundary_id(BoundaryIds::down_force); + } else { + cell->face(face_number)->set_boundary_id(BoundaryIds::no_force); + } + } + } + } + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) { + if (std::abs(cell->vertex(vertex_number)(0)) + std::abs(cell->vertex(vertex_number)(1)) < + 1e-10) { + cell->set_active_fe_index(1); + cell->set_material_id(MaterialIds::with_multiplier); + } + } + } + } + + } else if (dim == 3) { + GridGenerator::subdivided_hyper_L(triangulation, + {width_refine, height_refine, depth_refine}, + Point(0, 0, 0), + Point(width, height, depth), + {-1, -1, depth_refine}); + + triangulation.refine_global(Input::refinements); + + /*Set BCIDs */ + for (const auto &cell: dof_handler.active_cell_iterators()) { + if(cell->is_locally_owned()) + { + cell->set_active_fe_index(0); + cell->set_material_id(MaterialIds::without_multiplier); + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) { + if (cell->face(face_number)->at_boundary()) { + const auto center = cell->face(face_number)->center(); + + if (std::fabs(center(0) - downforce_x) < 1e-12) { + if (std::fabs(center(1) - downforce_y) < downforce_size) { + cell->face(face_number)->set_boundary_id(BoundaryIds::down_force); + if (std::fabs(center(2) - downforce_z) < downforce_size) { + cell->face(face_number)->set_boundary_id(BoundaryIds::down_force); + } else { + cell->face(face_number)->set_boundary_id(BoundaryIds::no_force); + } + } else { + cell->face(face_number)->set_boundary_id(BoundaryIds::no_force); + } + } + } + } + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) { + if (std::abs(cell->vertex(vertex_number)(0)) + std::abs(cell->vertex(vertex_number)(1)) < + 1e-10) { + cell->set_active_fe_index(1); + cell->set_material_id(MaterialIds::with_multiplier); + } + } + } + } + } else { + throw; + } + } else { + throw; + } + + dof_handler.distribute_dofs(fe_collection); + DoFRenumbering::component_wise(dof_handler, sub_blocks); + + dof_handler_displacement.distribute_dofs(fe_displacement); + dof_handler_displacement.distribute_mg_dofs(); + + displacement_to_system_dof_index_map.clear(); + +} + +/// Only individual points are given Dirichlet Boundary Conditions. +/// For example, in the MBB caes, The bottom corners are kept in place in the y direction +/// and the bottom left also in the x direction. +/// Because deal.ii is formulated to enforce boundary conditions along regions of the boundary, +/// we do this to ensure these BCs are only enforced at points. +template +void +KktSystem::setup_boundary_values() +{ + if (Input::geometry_base == GeometryOptions::mbb) + { + if (dim == 2) + { + for (const auto &cell: dof_handler.active_cell_iterators()) + { + if(cell->is_locally_owned()) + { + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) + { + if (cell->face(face_number)->at_boundary()) + { + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) + { + const auto vert = cell->vertex(vertex_number); + /*Find bottom left corner*/ + if (std::fabs(vert(0) - 0) < 1e-12 && std::fabs( + vert(1) - 0) < 1e-12) + { + + const unsigned int x_displacement = + cell->vertex_dof_index(vertex_number, 0, cell->active_fe_index()); + const unsigned int y_displacement = + cell->vertex_dof_index(vertex_number, 1, cell->active_fe_index()); + const unsigned int x_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 2, cell->active_fe_index()); + const unsigned int y_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 3, cell->active_fe_index()); + /*set bottom left BC*/ + boundary_values[x_displacement] = 0; + boundary_values[y_displacement] = 0; + boundary_values[x_displacement_multiplier] = 0; + boundary_values[y_displacement_multiplier] = 0; + } + /*Find bottom right corner*/ + if (std::fabs(vert(0) - 6) < 1e-12 && std::fabs( + vert(1) - 0) < 1e-12) + { +// const unsigned int x_displacement = +// cell->vertex_dof_index(vertex_number, 0, cell->active_fe_index()); + const unsigned int y_displacement = + cell->vertex_dof_index(vertex_number, 1, cell->active_fe_index()); +// const unsigned int x_displacement_multiplier = +// cell->vertex_dof_index(vertex_number, 2, cell->active_fe_index()); + const unsigned int y_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 3, cell->active_fe_index()); +// boundary_values[x_displacement] = 0; + boundary_values[y_displacement] = 0; +// boundary_values[x_displacement_multiplier] = 0; + boundary_values[y_displacement_multiplier] = 0; + } + } + } + } + } + + } + const unsigned int n_levels = triangulation.n_global_levels(); + level_dirichlet_boundary_dofs.resize(0,n_levels-1); + level_boundary_values.resize(0,n_levels-1); + mg_level_constraints.resize(0,n_levels-1); + + for(unsigned int level = 0; level < n_levels; ++level) + { + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler_displacement, + level, + relevant_dofs); + mg_level_constraints[level].reinit(relevant_dofs); + } + + + for (auto cell=dof_handler_displacement.begin_active(n_levels-1); + cell!=dof_handler_displacement.end_active(n_levels-1); + ++cell) + { + if(cell->is_locally_owned()) + { + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) + { + if (cell->face(face_number)->at_boundary()) + { + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) + { + const auto vert = cell->vertex(vertex_number); + /*Find bottom left corner*/ + if (std::fabs(vert(0) - 0) < 1e-12 && std::fabs( + vert(1) - 0) < 1e-12) + { + for (unsigned int level = 0; level < n_levels; ++level) + { + const unsigned int x_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 0, cell->active_fe_index()); + const unsigned int y_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 1, cell->active_fe_index()); + + /*set bottom left BC*/ + + level_dirichlet_boundary_dofs[level].insert(x_displacement); + level_dirichlet_boundary_dofs[level].insert(y_displacement); + + level_boundary_values[level][x_displacement] = 0; + level_boundary_values[level][y_displacement] = 0; + } + + } + /*Find bottom right corner*/ + if (std::fabs(vert(0) - 6) < 1e-12 && std::fabs( + vert(1) - 0) < 1e-12) + { + for (unsigned int level = 0; level < n_levels; ++level) + { + +// const unsigned int x_displacement = +// cell->mg_vertex_dof_index(level, vertex_number, 0,cell->active_fe_index()); + const unsigned int y_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 1, cell->active_fe_index()); + +// level_dirichlet_boundary_dofs[level].insert(x_displacement); + level_dirichlet_boundary_dofs[level].insert(y_displacement); + +// level_boundary_values[level][x_displacement] = 0; + level_boundary_values[level][y_displacement] = 0; + } + } + } + } + } + } + } + for (unsigned int level = 0; level < n_levels; ++level) + { + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler_displacement,level,relevant_dofs); + mg_level_constraints[level].add_lines(level_dirichlet_boundary_dofs[level]); + mg_level_constraints[level].make_consistent_in_parallel( + dof_handler_displacement.locally_owned_mg_dofs(level), + relevant_dofs, + mpi_communicator + ); + mg_level_constraints[level].close(); + } + + } else if (dim == 3) + { + pcout << "setting up BVs" << std::endl; + for (const auto &cell: dof_handler.active_cell_iterators()) + { + if(cell->is_locally_owned()) + { + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) + { + if (cell->face(face_number)->at_boundary()) + { + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) { + const auto vert = cell->vertex(vertex_number); + /*Find bottom left corner*/ + if (std::fabs(vert(0) - 0) < 1e-12 && + std::fabs(vert(1) - 0) < 1e-12 && + ((std::fabs(vert(2) - 0) < 1e-12) || (std::fabs(vert(2) - 1) < 1e-12))) + { + + + const unsigned int x_displacement = + cell->vertex_dof_index(vertex_number, 0, cell->active_fe_index()); + const unsigned int y_displacement = + cell->vertex_dof_index(vertex_number, 1, cell->active_fe_index()); + const unsigned int z_displacement = + cell->vertex_dof_index(vertex_number, 2, cell->active_fe_index()); + const unsigned int x_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 3, cell->active_fe_index()); + const unsigned int y_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 4, cell->active_fe_index()); + const unsigned int z_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 5, cell->active_fe_index()); + + boundary_values[x_displacement] = 0; + boundary_values[y_displacement] = 0; + boundary_values[z_displacement] = 0; + boundary_values[x_displacement_multiplier] = 0; + boundary_values[y_displacement_multiplier] = 0; + boundary_values[z_displacement_multiplier] = 0; + } + /*Find bottom right corner*/ + if (std::fabs(vert(0) - 6) < 1e-12 && + std::fabs(vert(1) - 0) < 1e-12 && + ((std::fabs(vert(2) - 0) < 1e-12) || (std::fabs(vert(2) - 1) < 1e-12))) + { + // const unsigned int x_displacement = + // cell->vertex_dof_index(vertex_number, 0, cell->active_fe_index()); + const unsigned int y_displacement = + cell->vertex_dof_index(vertex_number, 1, cell->active_fe_index()); + const unsigned int z_displacement = + cell->vertex_dof_index(vertex_number, 2, cell->active_fe_index()); + // const unsigned int x_displacement_multiplier = + // cell->vertex_dof_index(vertex_number, 3, cell->active_fe_index()); + const unsigned int y_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 4, cell->active_fe_index()); + const unsigned int z_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 5, cell->active_fe_index()); + // boundary_values[x_displacement] = 0; + boundary_values[y_displacement] = 0; + boundary_values[z_displacement] = 0; + // boundary_values[x_displacement_multiplier] = 0; + boundary_values[y_displacement_multiplier] = 0; + boundary_values[z_displacement_multiplier] = 0; + } + } + } + } + } + + } + + const unsigned int n_levels = triangulation.n_global_levels(); + level_dirichlet_boundary_dofs.resize(0,n_levels-1); + level_boundary_values.resize(0,n_levels-1); + mg_level_constraints.resize(0,n_levels-1); + + for(unsigned int level = 0; level < n_levels; ++level) + { + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler_displacement, + level, + relevant_dofs); + mg_level_constraints[level].reinit(relevant_dofs); + } + + for (auto cell=dof_handler_displacement.begin_active(n_levels-1); + cell!=dof_handler_displacement.end_active(n_levels-1); + ++cell) + { + if(cell->is_locally_owned()) + { + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) + { + if (cell->face(face_number)->at_boundary()) + { + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) + { + const auto vert = cell->vertex(vertex_number); + /*Find bottom left corner*/ + if (std::fabs(vert(0) - 0) < 1e-12 && + std::fabs(vert(1) - 0) < 1e-12 && + ((std::fabs(vert(2) - 0) < 1e-12) || (std::fabs(vert(2) - 1) < 1e-12))) + { + for (unsigned int level = 0; level < n_levels; ++level) + { + const unsigned int x_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 0, cell->active_fe_index()); + const unsigned int y_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 1, cell->active_fe_index()); + const unsigned int z_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 2, cell->active_fe_index()); + /*set bottom left BC*/ + level_boundary_values[level][x_displacement] = 0; + level_boundary_values[level][y_displacement] = 0; + level_boundary_values[level][z_displacement] = 0; + level_dirichlet_boundary_dofs[level].insert(x_displacement); + level_dirichlet_boundary_dofs[level].insert(y_displacement); + level_dirichlet_boundary_dofs[level].insert(z_displacement); + } + } + /*Find bottom right corner*/ + if (std::fabs(vert(0) - 6) < 1e-12 && + std::fabs(vert(1) - 0) < 1e-12 && + ((std::fabs(vert(2) - 0) < 1e-12) || + (std::fabs(vert(2) - 1) < 1e-12))) + { + for (unsigned int level = 0; level < n_levels; ++level) + { + const unsigned int y_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 1, cell->active_fe_index()); + const unsigned int z_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 2, cell->active_fe_index()); + level_boundary_values[level][y_displacement] = 0; + level_boundary_values[level][z_displacement] = 0; + + level_dirichlet_boundary_dofs[level].insert(y_displacement); + level_dirichlet_boundary_dofs[level].insert(z_displacement); + } + } + } + } + } + } + } + for (unsigned int level = 0; level < n_levels; ++level) + { + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler_displacement,level,relevant_dofs); + mg_level_constraints[level].add_lines(level_dirichlet_boundary_dofs[level]); + mg_level_constraints[level].make_consistent_in_parallel( + dof_handler_displacement.locally_owned_mg_dofs(level), + relevant_dofs, + mpi_communicator + ); + mg_level_constraints[level].close(); + } + + } + else + { + throw; + } + } else if (Input::geometry_base == GeometryOptions::l_shape) { + if (dim == 2) + { + for (const auto &cell: dof_handler.active_cell_iterators()) + { + if(cell->is_locally_owned()) + { + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) { + if (cell->face(face_number)->at_boundary()) { + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) { + const auto vert = cell->vertex(vertex_number); + /*Find top left corner*/ + if (std::fabs(vert(0) - 0) < 1e-12 && std::fabs( + vert(1) - 2) < 1e-12) { + + const unsigned int x_displacement = + cell->vertex_dof_index(vertex_number, 0, cell->active_fe_index()); + const unsigned int y_displacement = + cell->vertex_dof_index(vertex_number, 1, cell->active_fe_index()); + const unsigned int x_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 2, cell->active_fe_index()); + const unsigned int y_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 3, cell->active_fe_index()); + /*set bottom left BC*/ + boundary_values[x_displacement] = 0; + boundary_values[y_displacement] = 0; + boundary_values[x_displacement_multiplier] = 0; + boundary_values[y_displacement_multiplier] = 0; + } + /*Find top right corner*/ + if (std::fabs(vert(0) - 1) < 1e-12 && std::fabs( + vert(1) - 2) < 1e-12) { + const unsigned int x_displacement = + cell->vertex_dof_index(vertex_number, 0, cell->active_fe_index()); + const unsigned int y_displacement = + cell->vertex_dof_index(vertex_number, 1, cell->active_fe_index()); + const unsigned int x_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 2, cell->active_fe_index()); + const unsigned int y_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 3, cell->active_fe_index()); + boundary_values[x_displacement] = 0; + boundary_values[y_displacement] = 0; + boundary_values[x_displacement_multiplier] = 0; + boundary_values[y_displacement_multiplier] = 0; + } + } + } + } + } + } + const unsigned int n_levels = triangulation.n_global_levels(); + for (unsigned int level = 0; level < n_levels; ++level) + { + for (auto cell=dof_handler_displacement.begin_active(level); + cell!=dof_handler.end_active(level); + ++cell) + { + if(cell->is_locally_owned()) + { + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) + { + if (cell->face(face_number)->at_boundary()) + { + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) + { + const auto vert = cell->vertex(vertex_number); + /*Find bottom left corner*/ + if (std::fabs(vert(0) - 0) < 1e-12 && std::fabs( + vert(1) - 2) < 1e-12) + { + + const unsigned int x_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 0, 0); + const unsigned int y_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 1, 0); + /*set bottom left BC*/ + level_boundary_values[level][x_displacement] = 0; + level_boundary_values[level][y_displacement] = 0; + + level_dirichlet_boundary_dofs[level].insert(x_displacement); + level_dirichlet_boundary_dofs[level].insert(y_displacement); + + } + /*Find bottom right corner*/ + if (std::fabs(vert(0) - 1) < 1e-12 && std::fabs( + vert(1) - 2) < 1e-12) + { + const unsigned int x_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 0, 0); + const unsigned int y_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 1, 0); + level_boundary_values[level][x_displacement] = 0; + level_boundary_values[level][y_displacement] = 0; + + level_dirichlet_boundary_dofs[level].insert(x_displacement); + level_dirichlet_boundary_dofs[level].insert(y_displacement); + } + } + } + } + } + } + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler_displacement,level,relevant_dofs); + mg_level_constraints[level].add_lines(level_dirichlet_boundary_dofs[level]); + mg_level_constraints[level].make_consistent_in_parallel( + dof_handler_displacement.locally_owned_mg_dofs(level), + relevant_dofs, + mpi_communicator + ); + mg_level_constraints[level].close(); + } + + } + else if (dim == 3) + { + for (const auto &cell: dof_handler.active_cell_iterators()) { + if(cell->is_locally_owned()) + { + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) { + if (cell->face(face_number)->at_boundary()) { + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) { + const auto vert = cell->vertex(vertex_number); + /*Find bottom left corner*/ + if (std::fabs(vert(0) - 0) < 1e-12 && std::fabs( + vert(1) - 2) < 1e-12 && ((std::fabs( + vert(2) - 0) < 1e-12) || (std::fabs( + vert(2) - 1) < 1e-12))) { + + + const unsigned int x_displacement = + cell->vertex_dof_index(vertex_number, 0, cell->active_fe_index()); + const unsigned int y_displacement = + cell->vertex_dof_index(vertex_number, 1, cell->active_fe_index()); + const unsigned int z_displacement = + cell->vertex_dof_index(vertex_number, 2, cell->active_fe_index()); + const unsigned int x_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 3, cell->active_fe_index()); + const unsigned int y_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 4, cell->active_fe_index()); + const unsigned int z_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 5, cell->active_fe_index()); + + boundary_values[x_displacement] = 0; + boundary_values[y_displacement] = 0; + boundary_values[z_displacement] = 0; + boundary_values[x_displacement_multiplier] = 0; + boundary_values[y_displacement_multiplier] = 0; + boundary_values[z_displacement_multiplier] = 0; + } + /*Find bottom right corner*/ + if (std::fabs(vert(0) - 1) < 1e-12 && std::fabs( + vert(1) - 2) < 1e-12 && ((std::fabs( + vert(2) - 0) < 1e-12) || (std::fabs( + vert(2) - 1) < 1e-12))) { + const unsigned int x_displacement = + cell->vertex_dof_index(vertex_number, 0, cell->active_fe_index()); + const unsigned int y_displacement = + cell->vertex_dof_index(vertex_number, 1, cell->active_fe_index()); + const unsigned int z_displacement = + cell->vertex_dof_index(vertex_number, 2, cell->active_fe_index()); + const unsigned int x_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 3, cell->active_fe_index()); + const unsigned int y_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 4, cell->active_fe_index()); + const unsigned int z_displacement_multiplier = + cell->vertex_dof_index(vertex_number, 5, cell->active_fe_index()); + boundary_values[x_displacement] = 0; + boundary_values[y_displacement] = 0; + boundary_values[z_displacement] = 0; + boundary_values[x_displacement_multiplier] = 0; + boundary_values[y_displacement_multiplier] = 0; + boundary_values[z_displacement_multiplier] = 0; + } + } + } + } + } + } + const unsigned int n_levels = triangulation.n_global_levels(); + for (unsigned int level = 0; level < n_levels; ++level) + { + for (auto cell=dof_handler_displacement.begin_active(level); + cell!=dof_handler.end_active(level); + ++cell) + { + if(cell->is_locally_owned()) + { + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) + { + if (cell->face(face_number)->at_boundary()) + { + for (unsigned int vertex_number = 0; + vertex_number < GeometryInfo::vertices_per_cell; + ++vertex_number) + { + const auto vert = cell->vertex(vertex_number); + /*Find bottom left corner*/ + if (std::fabs(vert(0) - 0) < 1e-12 && std::fabs( + vert(1) - 2) < 1e-12 && ((std::fabs( + vert(2) - 0) < 1e-12) || (std::fabs( + vert(2) - 1) < 1e-12))) + { + + const unsigned int x_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 0, 0); + const unsigned int y_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 1, 0); + const unsigned int z_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 2, 0); + /*set bottom left BC*/ + level_boundary_values[level][x_displacement] = 0; + level_boundary_values[level][y_displacement] = 0; + level_boundary_values[level][z_displacement] = 0; + + level_dirichlet_boundary_dofs[level].insert(x_displacement); + level_dirichlet_boundary_dofs[level].insert(y_displacement); + level_dirichlet_boundary_dofs[level].insert(z_displacement); + } + /*Find bottom right corner*/ + if (std::fabs(vert(0) - 1) < 1e-12 && std::fabs( + vert(1) - 2) < 1e-12 && ((std::fabs( + vert(2) - 0) < 1e-12) || (std::fabs( + vert(2) - 1) < 1e-12))) + { + const unsigned int x_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 0, 0); + const unsigned int y_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 1, 0); + const unsigned int z_displacement = + cell->mg_vertex_dof_index(level, vertex_number, 2, 0); + level_boundary_values[level][x_displacement] = 0; + level_boundary_values[level][y_displacement] = 0; + level_boundary_values[level][z_displacement] = 0; + + level_dirichlet_boundary_dofs[level].insert(x_displacement); + level_dirichlet_boundary_dofs[level].insert(y_displacement); + level_dirichlet_boundary_dofs[level].insert(z_displacement); + } + } + } + } + } + } + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler_displacement,level,relevant_dofs); + mg_level_constraints[level].add_lines(level_dirichlet_boundary_dofs[level]); + mg_level_constraints[level].make_consistent_in_parallel( + dof_handler_displacement.locally_owned_mg_dofs(level), + relevant_dofs, + mpi_communicator + ); + mg_level_constraints[level].close(); + } + + } else { + throw; + } + } + + +} + + +///This makes a giant 10-by-10 block matrix that when assembled will represents the 10 KKT equations that +/// come from this problem, and also sets up the necessary block vectors. The +/// sparsity pattern for this matrix includes the sparsity pattern for the filter matrix. It also initializes +/// any block vectors we will use. +template +void +KktSystem::setup_block_system() { + + //MAKE n_u and n_P + + /*Setup 10 by 10 block matrix*/ + std::vector block_component(10, 2); + + block_component[0] = 0; + block_component[5] = 1; + + const std::vector dofs_per_block = + DoFTools::count_dofs_per_fe_block(dof_handler, block_component); + const unsigned int n_p = dofs_per_block[0]; + const unsigned int n_u = dofs_per_block[1]; + + pcout << "n_p: " << n_p << " n_u: " << n_u << std::endl; + + IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs(); + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + dsp.reinit(10, 10); + owned_partitioning.resize(10); + owned_partitioning[0] = dof_handler.locally_owned_dofs().get_view(0, n_p); + owned_partitioning[1] = dof_handler.locally_owned_dofs().get_view(n_p, 2 * n_p); + owned_partitioning[2] = dof_handler.locally_owned_dofs().get_view(2 * n_p, 3 * n_p); + owned_partitioning[3] = dof_handler.locally_owned_dofs().get_view(3 * n_p, 4 * n_p); + owned_partitioning[4] = dof_handler.locally_owned_dofs().get_view(4 * n_p, 5 * n_p); + owned_partitioning[5] = dof_handler.locally_owned_dofs().get_view(5 * n_p, 5 * n_p + n_u); + owned_partitioning[6] = dof_handler.locally_owned_dofs().get_view(5 * n_p + n_u, 5 * n_p + 2 * n_u); + owned_partitioning[7] = dof_handler.locally_owned_dofs().get_view(5 * n_p + 2 * n_u, 6 * n_p + 2 * n_u); + owned_partitioning[8] = dof_handler.locally_owned_dofs().get_view(6 * n_p + 2 * n_u, 7 * n_p + 2 * n_u); + owned_partitioning[9] = dof_handler.locally_owned_dofs().get_view(7 * n_p + 2 * n_u, 7 * n_p + 2 * n_u + 1); + relevant_partitioning.resize(10); + relevant_partitioning[0] = locally_relevant_dofs.get_view(0, n_p); + relevant_partitioning[1] = locally_relevant_dofs.get_view(n_p, 2 * n_p); + relevant_partitioning[2] = locally_relevant_dofs.get_view(2 * n_p, 3 * n_p); + relevant_partitioning[3] = locally_relevant_dofs.get_view(3 * n_p, 4 * n_p); + relevant_partitioning[4] = locally_relevant_dofs.get_view(4 * n_p, 5 * n_p); + relevant_partitioning[5] = locally_relevant_dofs.get_view(5 * n_p, 5 * n_p + n_u); + relevant_partitioning[6] = locally_relevant_dofs.get_view(5 * n_p + n_u, 5 * n_p + 2 * n_u); + relevant_partitioning[7] = locally_relevant_dofs.get_view(5 * n_p + 2 * n_u, 6 * n_p + 2 * n_u); + relevant_partitioning[8] = locally_relevant_dofs.get_view(6 * n_p + 2 * n_u, 7 * n_p + 2 * n_u); + relevant_partitioning[9] = locally_relevant_dofs.get_view(7 * n_p + 2 * n_u, 7 * n_p + 2 * n_u + 1); + + const std::vector block_sizes = {n_p, n_p, n_p, n_p, n_p, n_u, n_u, n_p, n_p, 1}; + + for (unsigned int k = 0; k < 10; k++) { + for (unsigned int j = 0; j < 10; j++) { + dsp.block(j, k).reinit(block_sizes[j], block_sizes[k]); + } + } + dsp.collect_sizes(); + Table<2, DoFTools::Coupling> coupling(2 * dim + 8, 2 * dim + 8); + //Coupling for density + coupling[SolutionComponents::density][SolutionComponents::density] = DoFTools::always; + + for (unsigned int i = 0; i < dim; i++) { + coupling[SolutionComponents::density][SolutionComponents::displacement + + i] = DoFTools::always; + coupling[SolutionComponents::displacement + + i][SolutionComponents::density] = DoFTools::always; + } + + coupling[SolutionComponents::density][SolutionComponents::unfiltered_density_multiplier] = DoFTools::always; + coupling[SolutionComponents::unfiltered_density_multiplier][SolutionComponents::density] = DoFTools::always; + + for (unsigned int i = 0; i < dim; i++) { + coupling[SolutionComponents::density][SolutionComponents::displacement_multiplier + + i] = DoFTools::always; + coupling[SolutionComponents::displacement_multiplier + + i][SolutionComponents::density] = DoFTools::always; + } + + //Coupling for displacement + for (unsigned int i = 0; i < dim; i++) { + + for (unsigned int k = 0; k < dim; k++) { + coupling[SolutionComponents::displacement + i][ + SolutionComponents::displacement_multiplier + + k] = DoFTools::always; + coupling[SolutionComponents::displacement_multiplier + k][ + SolutionComponents::displacement + + i] = DoFTools::always; + } + } + + // coupling for unfiltered density + coupling[SolutionComponents::unfiltered_density][SolutionComponents::density_lower_slack_multiplier] = DoFTools::always; + coupling[SolutionComponents::density_lower_slack_multiplier][SolutionComponents::unfiltered_density] = DoFTools::always; + + coupling[SolutionComponents::unfiltered_density][SolutionComponents::density_upper_slack_multiplier] = DoFTools::always; + coupling[SolutionComponents::density_upper_slack_multiplier][SolutionComponents::unfiltered_density] = DoFTools::always; + + coupling[SolutionComponents::unfiltered_density][SolutionComponents::unfiltered_density_multiplier] = DoFTools::always; + coupling[SolutionComponents::unfiltered_density_multiplier][SolutionComponents::unfiltered_density] = DoFTools::always; + + // Coupling for lower slack + coupling[SolutionComponents::density_lower_slack][SolutionComponents::density_lower_slack] = DoFTools::always; + + coupling[SolutionComponents::density_lower_slack][SolutionComponents::density_lower_slack_multiplier] = DoFTools::always; + coupling[SolutionComponents::density_lower_slack_multiplier][SolutionComponents::density_lower_slack] = DoFTools::always; + + // + coupling[SolutionComponents::density_upper_slack][SolutionComponents::density_upper_slack] = DoFTools::always; + coupling[SolutionComponents::density_upper_slack][SolutionComponents::density_upper_slack_multiplier] = DoFTools::always; + coupling[SolutionComponents::density_upper_slack_multiplier][SolutionComponents::density_upper_slack] = DoFTools::always; + + coupling[SolutionComponents::density_upper_slack_multiplier][SolutionComponents::density_upper_slack_multiplier] = DoFTools::always; + constraints.reinit(locally_relevant_dofs); + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler,constraints); + constraints.close(); + + system_matrix.clear(); + + // DoFTools::make_sparsity_pattern(dof_handler, coupling, dsp, constraints, false); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, dof_handler.locally_owned_dofs(), mpi_communicator, + locally_relevant_dofs); + //adds the row into the sparsity pattern for the total volume constraint + // for (const auto &cell: dof_handler.active_cell_iterators()) + // { + // if (cell->is_locally_owned()) + // { + // std::vector i(cell->get_fe().n_dofs_per_cell()); + // cell->get_dof_indices(i); + // dsp.block(SolutionBlocks::density, SolutionBlocks::total_volume_multiplier).add(i[cell->get_fe().component_to_system_index(0, 0)], 0); + // dsp.block(SolutionBlocks::total_volume_multiplier, SolutionBlocks::density).add(0, i[cell->get_fe().component_to_system_index(0, 0)]); + // } + // } + // Because of the single volume multiplier element only being on one processor, this works, and the above does not. + for (unsigned int i = 0; iis_locally_owned()) + { + std::vector i(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(i); + const unsigned int cell_index = i[cell->get_fe().component_to_system_index(0, 0)]; + for (const auto &neighbor_cell_index : density_filter.find_relevant_neighbors(cell_index)) + { + dsp.block(SolutionBlocks::unfiltered_density_multiplier, + SolutionBlocks::unfiltered_density).add(cell_index, neighbor_cell_index); + dsp.block(SolutionBlocks::unfiltered_density, + SolutionBlocks::unfiltered_density_multiplier).add(cell_index, neighbor_cell_index); + } + } + } + + SparsityTools::distribute_sparsity_pattern( + dsp, + Utilities::MPI::all_gather(mpi_communicator, + dof_handler.locally_owned_dofs()), + mpi_communicator, + locally_relevant_dofs); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + system_matrix.reinit(owned_partitioning, dsp, mpi_communicator); + + locally_relevant_solution.reinit(owned_partitioning, relevant_partitioning, mpi_communicator); + distributed_solution.reinit(owned_partitioning, mpi_communicator); + system_rhs.reinit(owned_partitioning, mpi_communicator); + + locally_relevant_solution.collect_sizes(); + distributed_solution.collect_sizes(); + system_rhs.collect_sizes(); + system_matrix.collect_sizes(); + IndexSet locally_owned_displacement_dofs = dof_handler_displacement.locally_owned_dofs(); + std::vector displacement_dof_indices; + std::vector system_dof_indices; + for (const auto &displacement_cell : dof_handler_displacement.active_cell_iterators()) + if (displacement_cell->is_locally_owned()) + { + typename DoFHandler::active_cell_iterator system_cell (&displacement_cell->get_triangulation(), + displacement_cell->level(), + displacement_cell->index(), + &dof_handler); + + displacement_dof_indices.resize (displacement_cell->get_fe().dofs_per_cell); + system_dof_indices.resize (system_cell->get_fe().dofs_per_cell); + + displacement_cell->get_dof_indices (displacement_dof_indices); + system_cell->get_dof_indices (system_dof_indices); + + for (unsigned int i=0; iget_fe().component_to_system_index( + displacement_cell->get_fe().system_to_component_index(i).first+SolutionComponents::displacement, + displacement_cell->get_fe().system_to_component_index(i).second + )]; + } + + } + } + const types::global_dof_index disp_start_index = system_matrix.get_row_indices().block_start( + SolutionBlocks::displacement); + for (auto &index_pair : displacement_to_system_dof_index_map) + index_pair.second -=disp_start_index; + for (auto &index_pair : displacement_to_system_dof_index_map) + { + if(index_pair.first != index_pair.second) + { + std::cout << "inexact matching for index: " << index_pair.first << " and " << index_pair.second << std::endl; + } + } + + +} + +///The equations describing the newtons method for finding 0s in the KKT conditions are implemented here. +template +void +KktSystem::assemble_block_system(const LA::MPI::BlockVector &distributed_state, const double barrier_size) { + /*Remove any values from old iterations*/ + + LA::MPI::BlockVector relevant_state(owned_partitioning, relevant_partitioning, mpi_communicator); + relevant_state = distributed_state; + + system_matrix.reinit(owned_partitioning, dsp, mpi_communicator); + locally_relevant_solution = 0; + system_rhs = 0; + + QGauss nine_quadrature(fe_nine.degree + 1); + QGauss ten_quadrature(fe_ten.degree + 1); + + hp::QCollection q_collection; + q_collection.push_back(nine_quadrature); + q_collection.push_back(ten_quadrature); + + hp::FEValues hp_fe_values(fe_collection, + q_collection, + update_values | update_quadrature_points | + update_JxW_values | update_gradients); + + QGauss common_face_quadrature(fe_ten.degree + 1); + + FEFaceValues fe_nine_face_values(fe_nine, + common_face_quadrature, + update_JxW_values | + update_gradients | update_values); + FEFaceValues fe_ten_face_values(fe_ten, + common_face_quadrature, + update_normal_vectors | + update_values); + + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + + const FEValuesExtractors::Scalar densities(SolutionComponents::density); + const FEValuesExtractors::Vector displacements(SolutionComponents::displacement); + const FEValuesExtractors::Scalar unfiltered_densities(SolutionComponents::unfiltered_density); + const FEValuesExtractors::Vector displacement_multipliers(SolutionComponents::displacement_multiplier); + const FEValuesExtractors::Scalar unfiltered_density_multipliers( + SolutionComponents::unfiltered_density_multiplier); + const FEValuesExtractors::Scalar density_lower_slacks(SolutionComponents::density_lower_slack); + const FEValuesExtractors::Scalar density_lower_slack_multipliers( + SolutionComponents::density_lower_slack_multiplier); + const FEValuesExtractors::Scalar density_upper_slacks(SolutionComponents::density_upper_slack); + const FEValuesExtractors::Scalar density_upper_slack_multipliers( + SolutionComponents::density_upper_slack_multiplier); + const FEValuesExtractors::Scalar total_volume_multiplier( + SolutionComponents::total_volume_multiplier); + + const Functions::ConstantFunction lambda(Input::material_lambda), mu(Input::material_mu); + + distributed_solution = distributed_state; + LA::MPI::BlockVector filtered_unfiltered_density_solution = distributed_solution; + LA::MPI::BlockVector filter_adjoint_unfiltered_density_multiplier_solution = distributed_solution; + filtered_unfiltered_density_solution.block(SolutionBlocks::unfiltered_density) = 0; + filter_adjoint_unfiltered_density_multiplier_solution.block(SolutionBlocks::unfiltered_density_multiplier) = 0; + density_filter.filter_matrix.vmult(filtered_unfiltered_density_solution.block(SolutionBlocks::unfiltered_density),distributed_solution.block(SolutionBlocks::unfiltered_density)); + density_filter.filter_matrix_transpose.vmult(filter_adjoint_unfiltered_density_multiplier_solution.block(SolutionBlocks::unfiltered_density_multiplier),distributed_solution.block(SolutionBlocks::unfiltered_density_multiplier)); + + LA::MPI::BlockVector relevant_filtered_unfiltered_density_solution = locally_relevant_solution; + LA::MPI::BlockVector relevant_filter_adjoint_unfiltered_density_multiplier_solution = locally_relevant_solution; + relevant_filtered_unfiltered_density_solution =filtered_unfiltered_density_solution; + relevant_filter_adjoint_unfiltered_density_multiplier_solution = filter_adjoint_unfiltered_density_multiplier_solution; + for (const auto &cell: dof_handler.active_cell_iterators()) { + if(cell->is_locally_owned()) + { + hp_fe_values.reinit(cell); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + cell_matrix.reinit(cell->get_fe().n_dofs_per_cell(), + cell->get_fe().n_dofs_per_cell()); + cell_rhs.reinit(cell->get_fe().n_dofs_per_cell()); + + const unsigned int n_q_points = fe_values.n_quadrature_points; + + std::vector old_density_values(n_q_points); + std::vector> old_displacement_values(n_q_points); + std::vector old_displacement_divs(n_q_points); + std::vector> old_displacement_symmgrads( + n_q_points); + std::vector> old_displacement_multiplier_values( + n_q_points); + std::vector old_displacement_multiplier_divs(n_q_points); + std::vector> old_displacement_multiplier_symmgrads( + n_q_points); + std::vector old_lower_slack_multiplier_values(n_q_points); + std::vector old_upper_slack_multiplier_values(n_q_points); + std::vector old_lower_slack_values(n_q_points); + std::vector old_upper_slack_values(n_q_points); + std::vector old_unfiltered_density_values(n_q_points); + std::vector old_unfiltered_density_multiplier_values(n_q_points); + std::vector filtered_unfiltered_density_values(n_q_points); + std::vector filter_adjoint_unfiltered_density_multiplier_values(n_q_points); + std::vector lambda_values(n_q_points); + std::vector mu_values(n_q_points); + + const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); + + cell_matrix = 0; + cell_rhs = 0; + local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(local_dof_indices); + + lambda.value_list(fe_values.get_quadrature_points(), lambda_values); + mu.value_list(fe_values.get_quadrature_points(), mu_values); + + fe_values[densities].get_function_values(relevant_state, + old_density_values); + fe_values[displacements].get_function_values(relevant_state, + old_displacement_values); + fe_values[displacements].get_function_divergences(relevant_state, + old_displacement_divs); + fe_values[displacements].get_function_symmetric_gradients( + relevant_state, old_displacement_symmgrads); + fe_values[displacement_multipliers].get_function_values( + relevant_state, old_displacement_multiplier_values); + fe_values[displacement_multipliers].get_function_divergences( + relevant_state, old_displacement_multiplier_divs); + fe_values[displacement_multipliers].get_function_symmetric_gradients( + relevant_state, old_displacement_multiplier_symmgrads); + fe_values[density_lower_slacks].get_function_values( + relevant_state, old_lower_slack_values); + fe_values[density_lower_slack_multipliers].get_function_values( + relevant_state, old_lower_slack_multiplier_values); + fe_values[density_upper_slacks].get_function_values( + relevant_state, old_upper_slack_values); + fe_values[density_upper_slack_multipliers].get_function_values( + relevant_state, old_upper_slack_multiplier_values); + fe_values[unfiltered_densities].get_function_values( + relevant_state, old_unfiltered_density_values); + fe_values[unfiltered_density_multipliers].get_function_values( + relevant_state, old_unfiltered_density_multiplier_values); + fe_values[unfiltered_densities].get_function_values( + relevant_filtered_unfiltered_density_solution, filtered_unfiltered_density_values); + fe_values[unfiltered_density_multipliers].get_function_values( + relevant_filter_adjoint_unfiltered_density_multiplier_solution, + filter_adjoint_unfiltered_density_multiplier_values); + + Tensor<1, dim> traction; + traction[1] = -1; + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { + + for (unsigned int i = 0; i < dofs_per_cell; ++i) { + const SymmetricTensor<2, dim> displacement_phi_i_symmgrad = + fe_values[displacements].symmetric_gradient(i, q_point); + const double displacement_phi_i_div = + fe_values[displacements].divergence(i, q_point); + + const SymmetricTensor<2, dim> displacement_multiplier_phi_i_symmgrad = + fe_values[displacement_multipliers].symmetric_gradient(i, + q_point); + const double displacement_multiplier_phi_i_div = + fe_values[displacement_multipliers].divergence(i, + q_point); + + + const double density_phi_i = fe_values[densities].value(i, + q_point); + const double unfiltered_density_phi_i = fe_values[unfiltered_densities].value(i, + q_point); + const double unfiltered_density_multiplier_phi_i = fe_values[unfiltered_density_multipliers].value( + i, q_point); + + const double lower_slack_multiplier_phi_i = + fe_values[density_lower_slack_multipliers].value(i, + q_point); + + const double lower_slack_phi_i = + fe_values[density_lower_slacks].value(i, q_point); + + const double upper_slack_phi_i = + fe_values[density_upper_slacks].value(i, q_point); + + const double upper_slack_multiplier_phi_i = + fe_values[density_upper_slack_multipliers].value(i, + q_point); + + + for (unsigned int j = 0; j < dofs_per_cell; ++j) { + const SymmetricTensor<2, dim> displacement_phi_j_symmgrad = + fe_values[displacements].symmetric_gradient(j, + q_point); + const double displacement_phi_j_div = + fe_values[displacements].divergence(j, q_point); + + const SymmetricTensor<2, dim> displacement_multiplier_phi_j_symmgrad = + fe_values[displacement_multipliers].symmetric_gradient( + j, q_point); + const double displacement_multiplier_phi_j_div = + fe_values[displacement_multipliers].divergence(j, + q_point); + + const double density_phi_j = fe_values[densities].value( + j, q_point); + + const double unfiltered_density_phi_j = fe_values[unfiltered_densities].value(j, + q_point); + const double unfiltered_density_multiplier_phi_j = fe_values[unfiltered_density_multipliers].value( + j, q_point); + + + const double lower_slack_phi_j = + fe_values[density_lower_slacks].value(j, q_point); + + const double upper_slack_phi_j = + fe_values[density_upper_slacks].value(j, q_point); + + const double lower_slack_multiplier_phi_j = + fe_values[density_lower_slack_multipliers].value(j, + q_point); + + const double upper_slack_multiplier_phi_j = + fe_values[density_upper_slack_multipliers].value(j, + q_point); + + //Equation 0 + cell_matrix(i, j) += + fe_values.JxW(q_point) * + ( + -density_phi_i * unfiltered_density_multiplier_phi_j + + - density_penalty_exponent * (density_penalty_exponent - 1) + * std::pow( + old_density_values[q_point], + density_penalty_exponent - 2) + * density_phi_i + * density_phi_j + * (old_displacement_multiplier_divs[q_point] * old_displacement_divs[q_point] + * lambda_values[q_point] + + 2 * mu_values[q_point] + * (old_displacement_symmgrads[q_point] * + old_displacement_multiplier_symmgrads[q_point])) + + - density_penalty_exponent * std::pow( + old_density_values[q_point], + density_penalty_exponent - 1) + * density_phi_i + * (displacement_multiplier_phi_j_div * old_displacement_divs[q_point] + * lambda_values[q_point] + + 2 * mu_values[q_point] + * + (old_displacement_symmgrads[q_point] * + displacement_multiplier_phi_j_symmgrad)) + + - density_penalty_exponent * std::pow( + old_density_values[q_point], + density_penalty_exponent - 1) + * density_phi_i + * (displacement_phi_j_div * old_displacement_multiplier_divs[q_point] + * lambda_values[q_point] + + 2 * mu_values[q_point] + * (old_displacement_multiplier_symmgrads[q_point] * + displacement_phi_j_symmgrad))); + //Equation 1 + + cell_matrix(i, j) += + fe_values.JxW(q_point) * ( + -density_penalty_exponent * std::pow( + old_density_values[q_point], + density_penalty_exponent - 1) + * density_phi_j + * (old_displacement_multiplier_divs[q_point] * displacement_phi_i_div + * lambda_values[q_point] + + 2 * mu_values[q_point] + * (old_displacement_multiplier_symmgrads[q_point] * + displacement_phi_i_symmgrad)) + + - std::pow(old_density_values[q_point], + density_penalty_exponent) + * (displacement_multiplier_phi_j_div * displacement_phi_i_div + * lambda_values[q_point] + + 2 * mu_values[q_point] + * (displacement_multiplier_phi_j_symmgrad * displacement_phi_i_symmgrad)) + + ); + + //Equation 2 has to do with the filter, which is calculated elsewhere. + cell_matrix(i, j) += + fe_values.JxW(q_point) * ( + -1 * unfiltered_density_phi_i * lower_slack_multiplier_phi_j + + unfiltered_density_phi_i * upper_slack_multiplier_phi_j); + + //Equation 3 - Primal Feasibility + + cell_matrix(i, j) += + fe_values.JxW(q_point) * ( + + -1 * density_penalty_exponent * std::pow( + old_density_values[q_point], + density_penalty_exponent - 1) + * density_phi_j + * (old_displacement_divs[q_point] * displacement_multiplier_phi_i_div + * lambda_values[q_point] + + 2 * mu_values[q_point] + * (old_displacement_symmgrads[q_point] * + displacement_multiplier_phi_i_symmgrad)) + + + -1 * std::pow(old_density_values[q_point], + density_penalty_exponent) + * (displacement_phi_j_div * displacement_multiplier_phi_i_div + * lambda_values[q_point] + + 2 * mu_values[q_point] + * + (displacement_phi_j_symmgrad * displacement_multiplier_phi_i_symmgrad))); + + //Equation 4 - more primal feasibility + cell_matrix(i, j) += + -1 * fe_values.JxW(q_point) * lower_slack_multiplier_phi_i * + (unfiltered_density_phi_j - lower_slack_phi_j); + + //Equation 5 - more primal feasibility + cell_matrix(i, j) += + -1 * fe_values.JxW(q_point) * upper_slack_multiplier_phi_i * ( + -1 * unfiltered_density_phi_j - upper_slack_phi_j); + + //Equation 6 - more primal feasibility - part with filter added later + cell_matrix(i, j) += + -1 * fe_values.JxW(q_point) * unfiltered_density_multiplier_phi_i * ( + density_phi_j); + + //Equation 7 - complementary slackness + cell_matrix(i, j) += fe_values.JxW(q_point) * + (lower_slack_phi_i * lower_slack_multiplier_phi_j + + lower_slack_phi_i * lower_slack_phi_j * + old_lower_slack_multiplier_values[q_point] / + old_lower_slack_values[q_point]); + //Equation 8 - complementary slackness + cell_matrix(i, j) += fe_values.JxW(q_point) * + (upper_slack_phi_i * upper_slack_multiplier_phi_j + + upper_slack_phi_i * upper_slack_phi_j * + old_upper_slack_multiplier_values[q_point] / + old_upper_slack_values[q_point]); + } + + } + } + + + MatrixTools::local_apply_boundary_values(boundary_values, local_dof_indices, + cell_matrix, cell_rhs, true); + + + constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); + + } + + } + // MPI_BARRIER(MPI_COMM_WORLD); + system_matrix.compress(VectorOperation::add); + system_rhs = calculate_rhs(distributed_state, barrier_size); + double cell_measure; + for (const auto &cell: dof_handler.active_cell_iterators()) + { + if(cell->is_locally_owned()) + { + std::vector i(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(i); + const unsigned int cell_i = i[cell->get_fe().component_to_system_index(0, 0)]; + + // typename LA::MPI::SparseMatrix::iterator iter = density_filter.filter_matrix.begin(cell_i); + for (const unsigned int j : density_filter.find_relevant_neighbors(cell_i)) + { + // unsigned int j = iter->column(); + double value = density_filter.filter_matrix(cell_i,j) * cell->measure(); + double value_transpose = density_filter.filter_matrix_transpose(cell_i,j) * cell->measure(); + + system_matrix.block(SolutionBlocks::unfiltered_density_multiplier, + SolutionBlocks::unfiltered_density).set(cell_i, j, value); + system_matrix.block(SolutionBlocks::unfiltered_density, + SolutionBlocks::unfiltered_density_multiplier).set(cell_i, j, value_transpose); + } + + cell_measure = cell->measure(); + + system_matrix.block(SolutionBlocks::density, SolutionBlocks::total_volume_multiplier).set(cell_i, 0, + cell->measure()); + system_matrix.block(SolutionBlocks::total_volume_multiplier,SolutionBlocks::density).set(0,cell_i, + cell->measure()); + } + } + system_matrix.compress(VectorOperation::insert); + // for (const auto &cell: dof_handler.active_cell_iterators()) + // { + // if(cell->is_locally_owned()) + // { + // std::vector i(cell->get_fe().n_dofs_per_cell()); + // cell->get_dof_indices(i); + // const unsigned int cell_i = i[cell->get_fe().component_to_system_index(0, 0)]; + + // for (const unsigned int j : density_filter.find_relevant_neighbors(cell_i)) { + // double value = system_matrix.block(SolutionBlocks::unfiltered_density_multiplier, + // SolutionBlocks::unfiltered_density).el(j,cell_i); + + // system_matrix.block(SolutionBlocks::unfiltered_density, + // SolutionBlocks::unfiltered_density_multiplier).set(j, cell_i, value); + // } + // } + // } + // system_matrix.compress(VectorOperation::insert); + + + pcout << "assembled " << std::endl; + +} + +///For use in the filter, this calculates the objective value we are working to minimize. +template +double +KktSystem::calculate_objective_value(const LA::MPI::BlockVector &distributed_state) const { + /*Remove any values from old iterations*/ + + locally_relevant_solution = distributed_state; + + + QGauss nine_quadrature(fe_nine.degree + 1); + QGauss ten_quadrature(fe_ten.degree + 1); + + hp::QCollection q_collection; + q_collection.push_back(nine_quadrature); + q_collection.push_back(ten_quadrature); + + hp::FEValues hp_fe_values(fe_collection, + q_collection, + update_values | update_quadrature_points | + update_JxW_values | update_gradients); + + QGauss common_face_quadrature(fe_ten.degree + 1); + + FEFaceValues fe_nine_face_values(fe_nine, + common_face_quadrature, + update_JxW_values | + update_gradients | update_values); + FEFaceValues fe_ten_face_values(fe_ten, + common_face_quadrature, + update_normal_vectors | + update_values); + + FullMatrix cell_matrix; + Vector cell_rhs; + + const FEValuesExtractors::Vector displacements(SolutionComponents::displacement); + + Tensor<1, dim> traction; + traction[1] = -1; + distributed_solution = distributed_state; + double objective_value = 0; + for (const auto &cell: dof_handler.active_cell_iterators()) + { + if(cell->is_locally_owned()) + { + hp_fe_values.reinit(cell); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); + const unsigned int n_q_points = fe_values.n_quadrature_points; + const unsigned int n_face_q_points = common_face_quadrature.size(); + + std::vector> old_displacement_values(n_q_points); + fe_values[displacements].get_function_values( + locally_relevant_solution, old_displacement_values); + + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) + { + if (cell->face(face_number)->at_boundary() && cell->face(face_number)->boundary_id() + == BoundaryIds::down_force) + { + for (unsigned int face_q_point = 0; + face_q_point < n_face_q_points; ++face_q_point) { + for (unsigned int i = 0; i < dofs_per_cell; ++i) { + if (cell->material_id() == MaterialIds::without_multiplier) { + fe_nine_face_values.reinit(cell, face_number); + objective_value += traction + * fe_nine_face_values[displacements].value(i, + face_q_point) + * fe_nine_face_values.JxW(face_q_point); + } else { + fe_ten_face_values.reinit(cell, face_number); + objective_value += traction + * fe_ten_face_values[displacements].value(i, + face_q_point) + * fe_ten_face_values.JxW(face_q_point); + } + } + } + } + } + } + } + double objective_value_out; + MPI_Allreduce(&objective_value, &objective_value_out, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + return objective_value; +} + + +///As the KKT System knows which vectors correspond to the slack variables, the sum of the logs of the slacks is computed here for use in the filter. +template +double +KktSystem::calculate_barrier_distance(const LA::MPI::BlockVector &state) const { + double barrier_distance_log_sum = 0; + unsigned int vect_size = state.block(SolutionBlocks::density_lower_slack).size(); + distributed_solution = state; + for (unsigned int k = 0; k < vect_size; k++) { + if (distributed_solution.block(SolutionBlocks::density_lower_slack).in_local_range(k)) + barrier_distance_log_sum += std::log(state.block(SolutionBlocks::density_lower_slack)[k]); + } + for (unsigned int k = 0; k < vect_size; k++) { + if (distributed_solution.block(SolutionBlocks::density_upper_slack).in_local_range(k)) + barrier_distance_log_sum += std::log(state.block(SolutionBlocks::density_upper_slack)[k]); + } + double out_barrier_distance_log_sum; + MPI_Allreduce(&barrier_distance_log_sum, &out_barrier_distance_log_sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + return out_barrier_distance_log_sum; +} + +///Calculates the norm of the RHS. While not the KKT norm, we also expect this to be 0 at a minimum. +template +double +KktSystem::calculate_rhs_norm(const LA::MPI::BlockVector &state, const double barrier_size) const { + return calculate_rhs(state, barrier_size).l2_norm(); +} + + +///Feasibility conditions appear on the RHS of the linear system, so I compute the RHS to find it. Could probably be combined with the objective value finding part to make it faster. +template +double +KktSystem::calculate_feasibility(const LA::MPI::BlockVector &state, const double barrier_size) const { + LA::MPI::BlockVector test_rhs = calculate_rhs(state, barrier_size); + + double norm = 0; + + distributed_solution = state; + double full_prod1 =1; + double full_prod2 = 1; + + for (unsigned int k = 0; k < state.block(SolutionBlocks::density_upper_slack).size(); k++) { + double prod1 = 1; + double prod2 = 1; + if(state.block(SolutionBlocks::density_upper_slack).in_local_range(k)) + { + prod1 = prod1 * state.block(SolutionBlocks::density_upper_slack)[k] + * state.block(SolutionBlocks::density_upper_slack)[k]; + } + if(state.block(SolutionBlocks::density_lower_slack).in_local_range(k)) + { + prod2 = prod2 * state.block(SolutionBlocks::density_lower_slack)[k] + * state.block(SolutionBlocks::density_lower_slack)[k]; + } + if(state.block(SolutionBlocks::density_upper_slack_multiplier).in_local_range(k)) + { + prod1 = prod1 * state.block(SolutionBlocks::density_upper_slack_multiplier)[k] + * state.block(SolutionBlocks::density_upper_slack_multiplier)[k]; + } + if(state.block(SolutionBlocks::density_lower_slack_multiplier).in_local_range(k)) + { + prod2 = prod2 * state.block(SolutionBlocks::density_lower_slack_multiplier)[k] + * state.block(SolutionBlocks::density_lower_slack_multiplier)[k]; + } + MPI_Allreduce(&prod1, &full_prod1, 1, MPI_DOUBLE, MPI_PROD, MPI_COMM_WORLD); + MPI_Allreduce(&prod2, &full_prod2, 1, MPI_DOUBLE, MPI_PROD, MPI_COMM_WORLD); + norm = norm + full_prod1 + full_prod2; + } + + norm += std::pow(test_rhs.block(SolutionBlocks::displacement).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::density).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::unfiltered_density).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::displacement_multiplier).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::unfiltered_density_multiplier).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::total_volume_multiplier).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::density_upper_slack_multiplier).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::density_lower_slack_multiplier).l2_norm(), 2); + + return norm; +} + +///calculates the KKT norm of the system, representing how close the program is to convergence. +template +double +KktSystem::calculate_convergence(const LA::MPI::BlockVector &state) const { + LA::MPI::BlockVector test_rhs = calculate_rhs(state, Input::min_barrier_size); + double norm = 0; + + distributed_solution = state; + double full_prod1 =1; + double full_prod2 = 1; + + for (unsigned int k = 0; k < state.block(SolutionBlocks::density_upper_slack).size(); k++) { + double prod1 = 1; + double prod2 = 1; + if(state.block(SolutionBlocks::density_upper_slack).in_local_range(k)) + { + prod1 = prod1 * state.block(SolutionBlocks::density_upper_slack)[k] + * state.block(SolutionBlocks::density_upper_slack)[k]; + } + if(state.block(SolutionBlocks::density_lower_slack).in_local_range(k)) + { + prod2 = prod2 * state.block(SolutionBlocks::density_lower_slack)[k] + * state.block(SolutionBlocks::density_lower_slack)[k]; + } + if(state.block(SolutionBlocks::density_upper_slack_multiplier).in_local_range(k)) + { + prod1 = prod1 * state.block(SolutionBlocks::density_upper_slack_multiplier)[k] + * state.block(SolutionBlocks::density_upper_slack_multiplier)[k]; + } + if(state.block(SolutionBlocks::density_lower_slack_multiplier).in_local_range(k)) + { + prod2 = prod2 * state.block(SolutionBlocks::density_lower_slack_multiplier)[k] + * state.block(SolutionBlocks::density_lower_slack_multiplier)[k]; + } + MPI_Allreduce(&prod1, &full_prod1, 1, MPI_DOUBLE, MPI_PROD, MPI_COMM_WORLD); + MPI_Allreduce(&prod2, &full_prod2, 1, MPI_DOUBLE, MPI_PROD, MPI_COMM_WORLD); + norm = norm + full_prod1 + full_prod2; + } + + + norm += std::pow(test_rhs.block(SolutionBlocks::displacement).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::density).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::unfiltered_density).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::displacement_multiplier).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::unfiltered_density_multiplier).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::total_volume_multiplier).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::density_upper_slack_multiplier).l2_norm(), 2); + norm += std::pow(test_rhs.block(SolutionBlocks::density_lower_slack_multiplier).l2_norm(), 2); + + norm = std::pow(norm, .5); + + pcout << "KKT norm: " << norm << std::endl; + return norm; +} + +/// Makes the RHS of the KKT equations +template +LA::MPI::BlockVector +KktSystem::calculate_rhs(const LA::MPI::BlockVector &distributed_state, const double barrier_size) const { + LA::MPI::BlockVector test_rhs (system_rhs); + LA::MPI::BlockVector state (locally_relevant_solution); + state = distributed_state; + test_rhs = 0.; + + QGauss nine_quadrature(fe_nine.degree + 1); + QGauss ten_quadrature(fe_ten.degree + 1); + + hp::QCollection q_collection; + q_collection.push_back(nine_quadrature); + q_collection.push_back(ten_quadrature); + + hp::FEValues hp_fe_values(fe_collection, + q_collection, + update_values | update_quadrature_points | + update_JxW_values | update_gradients); + + QGauss common_face_quadrature(fe_ten.degree + 1); + + FEFaceValues fe_nine_face_values(fe_nine, + common_face_quadrature, + update_JxW_values | + update_gradients | update_values); + FEFaceValues fe_ten_face_values(fe_ten, + common_face_quadrature, + update_normal_vectors | + update_values); + + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + + const FEValuesExtractors::Scalar densities(SolutionComponents::density); + const FEValuesExtractors::Vector displacements(SolutionComponents::displacement); + const FEValuesExtractors::Scalar unfiltered_densities(SolutionComponents::unfiltered_density); + const FEValuesExtractors::Vector displacement_multipliers(SolutionComponents::displacement_multiplier); + const FEValuesExtractors::Scalar unfiltered_density_multipliers( + SolutionComponents::unfiltered_density_multiplier); + const FEValuesExtractors::Scalar density_lower_slacks(SolutionComponents::density_lower_slack); + const FEValuesExtractors::Scalar density_lower_slack_multipliers( + SolutionComponents::density_lower_slack_multiplier); + const FEValuesExtractors::Scalar density_upper_slacks(SolutionComponents::density_upper_slack); + const FEValuesExtractors::Scalar density_upper_slack_multipliers( + SolutionComponents::density_upper_slack_multiplier); + const FEValuesExtractors::Scalar total_volume_multiplier( + SolutionComponents::total_volume_multiplier); + + + const unsigned int n_face_q_points = common_face_quadrature.size(); + + const Functions::ConstantFunction lambda(1.), mu(1.); + + locally_relevant_solution = state; + distributed_solution = state; + LA::MPI::BlockVector filtered_unfiltered_density_solution (distributed_solution); + LA::MPI::BlockVector filter_adjoint_unfiltered_density_multiplier_solution (distributed_solution); + filtered_unfiltered_density_solution.block(SolutionBlocks::unfiltered_density) = 0; + filter_adjoint_unfiltered_density_multiplier_solution.block(SolutionBlocks::unfiltered_density_multiplier) = 0; + + density_filter.filter_matrix.vmult(filtered_unfiltered_density_solution.block(SolutionBlocks::unfiltered_density),distributed_solution.block(SolutionBlocks::unfiltered_density)); + density_filter.filter_matrix_transpose.vmult(filter_adjoint_unfiltered_density_multiplier_solution.block(SolutionBlocks::unfiltered_density_multiplier),distributed_solution.block(SolutionBlocks::unfiltered_density_multiplier)); + + LA::MPI::BlockVector relevant_filtered_unfiltered_density_solution (locally_relevant_solution); + LA::MPI::BlockVector relevant_filter_adjoint_unfiltered_density_multiplier_solution (locally_relevant_solution); + relevant_filtered_unfiltered_density_solution = filtered_unfiltered_density_solution; + relevant_filter_adjoint_unfiltered_density_multiplier_solution = filter_adjoint_unfiltered_density_multiplier_solution; + + double old_volume_multiplier_temp = 0; + double old_volume_multiplier; + if(distributed_state.block(SolutionBlocks::total_volume_multiplier).in_local_range(0)) + { + old_volume_multiplier_temp = state.block(SolutionBlocks::total_volume_multiplier)[0]; + } + MPI_Allreduce(&old_volume_multiplier_temp, &old_volume_multiplier, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + for (const auto &cell: dof_handler.active_cell_iterators()) { + if(cell->is_locally_owned()) + { + hp_fe_values.reinit(cell); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + cell_matrix.reinit(cell->get_fe().n_dofs_per_cell(), + cell->get_fe().n_dofs_per_cell()); + cell_rhs.reinit(cell->get_fe().n_dofs_per_cell()); + + const unsigned int n_q_points = fe_values.n_quadrature_points; + + std::vector old_density_values(n_q_points); + std::vector> old_displacement_values(n_q_points); + std::vector old_displacement_divs(n_q_points); + std::vector> old_displacement_symmgrads( + n_q_points); + std::vector> old_displacement_multiplier_values( + n_q_points); + std::vector old_displacement_multiplier_divs(n_q_points); + std::vector> old_displacement_multiplier_symmgrads( + n_q_points); + std::vector old_lower_slack_multiplier_values(n_q_points); + std::vector old_upper_slack_multiplier_values(n_q_points); + std::vector old_lower_slack_values(n_q_points); + std::vector old_upper_slack_values(n_q_points); + std::vector old_unfiltered_density_values(n_q_points); + std::vector old_unfiltered_density_multiplier_values(n_q_points); + std::vector filtered_unfiltered_density_values(n_q_points); + std::vector filter_adjoint_unfiltered_density_multiplier_values(n_q_points); + std::vector lambda_values(n_q_points); + std::vector mu_values(n_q_points); + + const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); + + cell_matrix = 0; + cell_rhs = 0; + local_dof_indices.resize(cell->get_fe().n_dofs_per_cell()); + + cell->get_dof_indices(local_dof_indices); + + lambda.value_list(fe_values.get_quadrature_points(), lambda_values); + mu.value_list(fe_values.get_quadrature_points(), mu_values); + + fe_values[densities].get_function_values(state, + old_density_values); + fe_values[displacements].get_function_values(state, + old_displacement_values); + fe_values[displacements].get_function_divergences(state, + old_displacement_divs); + fe_values[displacements].get_function_symmetric_gradients( + state, old_displacement_symmgrads); + fe_values[displacement_multipliers].get_function_values( + state, old_displacement_multiplier_values); + fe_values[displacement_multipliers].get_function_divergences( + state, old_displacement_multiplier_divs); + fe_values[displacement_multipliers].get_function_symmetric_gradients( + state, old_displacement_multiplier_symmgrads); + fe_values[density_lower_slacks].get_function_values( + state, old_lower_slack_values); + fe_values[density_lower_slack_multipliers].get_function_values( + state, old_lower_slack_multiplier_values); + fe_values[density_upper_slacks].get_function_values( + state, old_upper_slack_values); + fe_values[density_upper_slack_multipliers].get_function_values( + state, old_upper_slack_multiplier_values); + fe_values[unfiltered_densities].get_function_values( + state, old_unfiltered_density_values); + fe_values[unfiltered_density_multipliers].get_function_values( + state, old_unfiltered_density_multiplier_values); + fe_values[unfiltered_densities].get_function_values( + relevant_filtered_unfiltered_density_solution, filtered_unfiltered_density_values); + fe_values[unfiltered_density_multipliers].get_function_values( + relevant_filter_adjoint_unfiltered_density_multiplier_solution, + filter_adjoint_unfiltered_density_multiplier_values); + + + Tensor<1, dim> traction; + traction[1] = -1; + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { + + for (unsigned int i = 0; i < dofs_per_cell; ++i) { + const SymmetricTensor<2, dim> displacement_phi_i_symmgrad = + fe_values[displacements].symmetric_gradient(i, q_point); + const double displacement_phi_i_div = + fe_values[displacements].divergence(i, q_point); + + const SymmetricTensor<2, dim> displacement_multiplier_phi_i_symmgrad = + fe_values[displacement_multipliers].symmetric_gradient(i, + q_point); + const double displacement_multiplier_phi_i_div = + fe_values[displacement_multipliers].divergence(i, + q_point); + + + const double density_phi_i = fe_values[densities].value(i, + q_point); + const double unfiltered_density_phi_i = fe_values[unfiltered_densities].value(i, + q_point); + const double unfiltered_density_multiplier_phi_i = fe_values[unfiltered_density_multipliers].value( + i, q_point); + + const double lower_slack_multiplier_phi_i = + fe_values[density_lower_slack_multipliers].value(i, + q_point); + + const double lower_slack_phi_i = + fe_values[density_lower_slacks].value(i, q_point); + + const double upper_slack_phi_i = + fe_values[density_upper_slacks].value(i, q_point); + + const double upper_slack_multiplier_phi_i = + fe_values[density_upper_slack_multipliers].value(i, + q_point); + + //rhs eqn 0 + cell_rhs(i) += + -1 * fe_values.JxW(q_point) * ( + -1 * density_penalty_exponent * + std::pow(old_density_values[q_point], density_penalty_exponent - 1) * density_phi_i + * (old_displacement_multiplier_divs[q_point] * old_displacement_divs[q_point] + * lambda_values[q_point] + + 2 * mu_values[q_point] * (old_displacement_symmgrads[q_point] + * old_displacement_multiplier_symmgrads[q_point])) + - density_phi_i * old_unfiltered_density_multiplier_values[q_point] + + old_volume_multiplier * density_phi_i + ); + + //rhs eqn 1 - boundary terms counted later + cell_rhs(i) += + -1 * fe_values.JxW(q_point) + * ( + -1 * std::pow(old_density_values[q_point], density_penalty_exponent) + * ( + old_displacement_multiplier_divs[q_point] * displacement_phi_i_div + * lambda_values[q_point] + + 2 * mu_values[q_point] * (old_displacement_multiplier_symmgrads[q_point] + * displacement_phi_i_symmgrad) + ) + ); + + //rhs eqn 2 + cell_rhs(i) += + -1 * fe_values.JxW(q_point) * ( + unfiltered_density_phi_i * + filter_adjoint_unfiltered_density_multiplier_values[q_point] + + unfiltered_density_phi_i * old_upper_slack_multiplier_values[q_point] + + -1 * unfiltered_density_phi_i * old_lower_slack_multiplier_values[q_point] + ); + + //rhs eqn 3 - boundary terms counted later + cell_rhs(i) += + -1 * fe_values.JxW(q_point) * ( + -1 * std::pow(old_density_values[q_point], density_penalty_exponent) + * (old_displacement_divs[q_point] * displacement_multiplier_phi_i_div + * lambda_values[q_point] + + 2 * mu_values[q_point] * (displacement_multiplier_phi_i_symmgrad + * old_displacement_symmgrads[q_point])) + ); + + //rhs eqn 4 + cell_rhs(i) += + -1 * fe_values.JxW(q_point) * + (-1 * lower_slack_multiplier_phi_i + * (old_unfiltered_density_values[q_point] - old_lower_slack_values[q_point]) + ); + + //rhs eqn 5 + cell_rhs(i) += + -1 * fe_values.JxW(q_point) * ( + -1 * upper_slack_multiplier_phi_i + * (1 - old_unfiltered_density_values[q_point] + - old_upper_slack_values[q_point])); + + //rhs eqn 6 + if (std::abs(old_density_values[q_point] - filtered_unfiltered_density_values[q_point])>1e-12) + { + cell_rhs(i) += + -1 * fe_values.JxW(q_point) * ( + -1 * unfiltered_density_multiplier_phi_i + * (old_density_values[q_point] - filtered_unfiltered_density_values[q_point]) + ); + } + + + //rhs eqn 7 + cell_rhs(i) += + -1 * fe_values.JxW(q_point) * + (lower_slack_phi_i * + (old_lower_slack_multiplier_values[q_point] - + barrier_size / old_lower_slack_values[q_point])); + + //rhs eqn 8 + cell_rhs(i) += + -1 * fe_values.JxW(q_point) * + (upper_slack_phi_i * + (old_upper_slack_multiplier_values[q_point] - + barrier_size / old_upper_slack_values[q_point])); + + } + + } + + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) { + if (cell->face(face_number)->at_boundary() && cell->face( + face_number)->boundary_id() == BoundaryIds::down_force) { + for (unsigned int face_q_point = 0; + face_q_point < n_face_q_points; ++face_q_point) { + for (unsigned int i = 0; i < dofs_per_cell; ++i) { + if (cell->material_id() == MaterialIds::without_multiplier) { + fe_nine_face_values.reinit(cell, face_number); + cell_rhs(i) += -1 + * traction + * fe_nine_face_values[displacements].value(i, + face_q_point) + * fe_nine_face_values.JxW(face_q_point); + + cell_rhs(i) += -1 * traction + * fe_nine_face_values[displacement_multipliers].value( + i, face_q_point) + * fe_nine_face_values.JxW(face_q_point); + } else { + fe_ten_face_values.reinit(cell, face_number); + cell_rhs(i) += -1 + * traction + * fe_ten_face_values[displacements].value(i, + face_q_point) + * fe_ten_face_values.JxW(face_q_point); + + cell_rhs(i) += -1 * traction + * fe_ten_face_values[displacement_multipliers].value( + i, face_q_point) + * fe_ten_face_values.JxW(face_q_point); + } + } + } + } + } + + + MatrixTools::local_apply_boundary_values(boundary_values, local_dof_indices, + cell_matrix, cell_rhs, true); + constraints.distribute_local_to_global( + cell_rhs, local_dof_indices, test_rhs); + } + } + + test_rhs.compress(VectorOperation::add); + double total_volume_temp = 0; + double goal_volume_temp = 0; + double total_volume, goal_volume; + + distributed_solution = state; + for (const auto &cell: dof_handler.active_cell_iterators()) { + if(cell->is_locally_owned()) + { + std::vector i(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(i); + if (distributed_solution.block(SolutionBlocks::density).in_local_range(i[cell->get_fe().component_to_system_index(0, 0)])) + { + total_volume_temp += cell->measure() * state.block(SolutionBlocks::density)[i[cell->get_fe().component_to_system_index(0, 0)]]; + goal_volume_temp += cell->measure() * Input::volume_percentage; + } + } + } + + MPI_Allreduce(&total_volume_temp, &total_volume, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&goal_volume_temp, &goal_volume, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + + if (test_rhs.block(SolutionBlocks::total_volume_multiplier).in_local_range(0)) + { + test_rhs.block(SolutionBlocks::total_volume_multiplier)[0] = goal_volume - total_volume; + } + test_rhs.compress(VectorOperation::insert); + + return test_rhs; + +} + +///Solves the big system to get the newton step +template +LA::MPI::BlockVector +KktSystem::solve(const LA::MPI::BlockVector &state) { + double gmres_tolerance; + if (Input::use_eisenstat_walker) { + gmres_tolerance = std::max( + std::min( + .1 * system_rhs.l2_norm() / (initial_rhs_error), + .001 + ), + Input::default_gmres_tolerance); + } + else { + gmres_tolerance = Input::default_gmres_tolerance*system_rhs.l2_norm(); + } + + + locally_relevant_solution=state; + distributed_solution = state; + + SolverControl solver_control(10000, gmres_tolerance); + + // ************ BEGIN MAKING MF GMG ELASTICITY PRECONDITIONER *************************** + using SystemMFMatrixType = MF_Elasticity_Operator; + using LevelMFMatrixType = MF_Elasticity_Operator; + + elasticity_matrix_mf.clear(); + mg_matrices.clear_elements(); + + std::map< types::global_dof_index, Point< dim > > support_points; + std::map< types::global_dof_index, Point< dim > > support_points_displacement; + + MappingQGeneric generic_map_displacement(1); + MappingQGeneric generic_map_1(1); + MappingQGeneric generic_map_2(1); + + hp::MappingCollection< dim, dim > hp_generic_map; + + hp_generic_map.push_back(generic_map_1); + hp_generic_map.push_back(generic_map_2); + + DoFTools::map_dofs_to_support_points(generic_map_displacement, dof_handler_displacement, support_points_displacement); + DoFTools::map_dofs_to_support_points(hp_generic_map, dof_handler, support_points); + + + const types::global_dof_index disp_mult_start_index = system_matrix.get_row_indices().block_start(SolutionBlocks::displacement_multiplier); + + for (const auto &support_points_displacement_pair : support_points_displacement) + { + if (support_points_displacement_pair.second != support_points[support_points_displacement_pair.first+disp_mult_start_index]) + pcout << "d = " << support_points_displacement_pair.first << ", points are " << support_points_displacement_pair.second << " and " << support_points[support_points_displacement_pair.first+disp_mult_start_index] << std::endl; + } + + MPI_Barrier(MPI_COMM_WORLD); + + std::vector locally_owned_dofs = Utilities::MPI::all_gather(mpi_communicator, dof_handler_displacement.locally_owned_dofs()); + IndexSet locally_active_dofs; + DoFTools::extract_locally_active_dofs(dof_handler_displacement, locally_active_dofs); + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof_handler_displacement, locally_relevant_dofs); + // AffineConstraints temp_displacement_constraints; + + + + displacement_constraints.clear(); + displacement_constraints.reinit(locally_relevant_dofs); //FIXME SHOULD THIS BE RELEVANT??? + displacement_constraints.copy_from(mg_level_constraints[triangulation.n_global_levels()-1]); + std::cout << "displacement constraint number: " << displacement_constraints.n_constraints() <::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points); + std::shared_ptr> system_mf_storage( + new MatrixFree()); + system_mf_storage->reinit(generic_map_1, + dof_handler_displacement, + displacement_constraints, + QGauss<1>(fe_displacement.degree + 1), + additional_data); + elasticity_matrix_mf.initialize(system_mf_storage); + } + + + LinearAlgebra::distributed::Vector distributed_displacement_sol; + LinearAlgebra::distributed::Vector distributed_displacement_rhs; + + elasticity_matrix_mf.initialize_dof_vector(distributed_displacement_sol); + elasticity_matrix_mf.initialize_dof_vector(distributed_displacement_rhs); + + ChangeVectorTypes::copy_from_system_to_displacement_vector(distributed_displacement_sol,distributed_solution.block(SolutionBlocks::displacement),displacement_to_system_dof_index_map); + + ChangeVectorTypes::copy_from_system_to_displacement_vector(distributed_displacement_rhs,system_rhs.block(SolutionBlocks::displacement),displacement_to_system_dof_index_map); + + const unsigned int n_levels = triangulation.n_global_levels(); + mg_matrices.resize(0, n_levels - 1); + + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(dof_handler_displacement); + const std::set empty_boundary_set; + // mg_constrained_dofs.make_zero_boundary_constraints(dof_handler_displacement,empty_boundary_set); + + + for (unsigned int level = 0; level < n_levels; ++level) + { + mg_constrained_dofs.add_user_constraints(level, mg_level_constraints[level]); + } + + // for (unsigned int level = 0; level < n_levels; ++level) + // { + // mg_level_constraints[level].print(std::cout); + // } + + for (unsigned int level = 0; level < n_levels; ++level) + { + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler_displacement, + level, + relevant_dofs); + + typename MatrixFree::AdditionalData additional_data; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + additional_data.mapping_update_flags = + (update_gradients | update_JxW_values | update_quadrature_points); + additional_data.mg_level = level; + + std::shared_ptr> mg_mf_storage_level( + new MatrixFree()); + mg_mf_storage_level->reinit(generic_map_1, + dof_handler_displacement, + mg_level_constraints[level], + QGauss<1>(fe_displacement.degree + 1), + additional_data); + + mg_matrices[level].clear(); + mg_matrices[level].initialize(mg_mf_storage_level, + mg_constrained_dofs, + level); + } + + + //+++++++++++++++++++++++++EVALUATE MATRIX LEVEL DENSITIES HERE +++++++++++++++++++++++++++++++++++++ + + + dof_handler_density.distribute_dofs(fe_density); + + DoFRenumbering::component_wise(dof_handler_density); + DoFRenumbering::hierarchical(dof_handler_density); + + dof_handler_density.distribute_mg_dofs(); + + active_density_vector.reinit(dof_handler_density.locally_owned_dofs(),triangulation.get_communicator()); + + ChangeVectorTypes::copy(active_density_vector,distributed_solution.block(SolutionBlocks::density)); + + + const unsigned int n_cells = elasticity_matrix_mf.get_matrix_free()->n_cell_batches(); + // { + + + // QGauss nine_quadrature(2); + // QGauss ten_quadrature(2); + + // hp::QCollection q_collection; + // q_collection.push_back(nine_quadrature); + // q_collection.push_back(ten_quadrature); + + // hp::FEValues hp_fe_values(fe_collection, + // q_collection, + // update_values | update_quadrature_points | + // update_JxW_values | update_gradients); + + + + // for (const auto &cell : dof_handler.active_cell_iterators()) + // if (cell->is_locally_owned()) + // { + + // hp_fe_values.reinit(cell); + // const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + + // const unsigned int dofs_per_cell = fe_values.dofs_per_cell; + // const unsigned int n_q_points = fe_values.n_quadrature_points; + + // std::vector local_dof_indices (dofs_per_cell); + // cell->get_dof_indices (local_dof_indices); + // Vector cell_vector (dofs_per_cell); + // Vector local_projection (dofs_per_cell); + // FullMatrix local_mass_matrix (dofs_per_cell, dofs_per_cell); + + // std::vector rhs_values(n_q_points); + // std::vector old_density_values(n_q_points); + + // const FEValuesExtractors::Scalar densities(SolutionComponents::density); + // fe_values[densities].get_function_values(locally_relevant_solution, old_density_values); + // double cell_density = old_density_values[0]; + + // for (unsigned int i=0; i i(cell->get_fe().n_dofs_per_cell()); + // cell->get_dof_indices(i); + // const unsigned int i_val = i[cell->get_fe().component_to_system_index(0, 0)]; + // active_density_vector[i_val] = cell_density; + + + // } + + // // active_density_vector.compress(VectorOperation::insert); + // } + // MAKE ACTIVE_CELL_DATA + std::vector local_dof_indices(fe_density.dofs_per_cell); + active_cell_data.density.reinit(TableIndices<2>(n_cells, 1)); + for (unsigned int cell=0; celln_active_entries_per_cell_batch(cell); + + for (unsigned int i=0; i::active_cell_iterator FEQ_cell =elasticity_matrix_mf.get_matrix_free()->get_cell_iterator(cell,i); + typename DoFHandler::active_cell_iterator DG_cell(&(triangulation), + FEQ_cell->level(), + FEQ_cell->index(), + &dof_handler_density); + + DG_cell->get_active_or_mg_dof_indices(local_dof_indices); + + active_cell_data.density(cell, 0)[i] = active_density_vector(local_dof_indices[0]); + } + } + + elasticity_matrix_mf.set_cell_data(active_cell_data); + + //MAKE LEVEL DENSITY VECTOR + + level_cell_data.resize(0,n_levels-1); + level_density_vector = 0.; + level_density_vector.resize(0,n_levels-1); + + transfer.build(dof_handler_density); + + transfer.interpolate_to_mg(dof_handler_density, + level_density_vector, + active_density_vector); + + // MAKE LEVEL_CELL_DATA + for (unsigned int level=0; leveln_cell_batches(); + + level_cell_data[level].density.reinit(TableIndices<2>(n_cells, 1)); + for (unsigned int cell=0; celln_active_entries_per_cell_batch(cell); + for (unsigned int i=0; i::level_cell_iterator FEQ_cell = mg_matrices[level].get_matrix_free()->get_cell_iterator(cell,i); + typename DoFHandler::level_cell_iterator DG_cell(&(triangulation), + FEQ_cell->level(), + FEQ_cell->index(), + &dof_handler_density); + DG_cell->get_active_or_mg_dof_indices(local_dof_indices); + + + level_cell_data[level].density(cell, 0)[i] = level_density_vector[level](local_dof_indices[0]); + } + + } + + // Store density tables and other data into the multigrid level matrix-free objects. + + mg_matrices[level].set_cell_data (level_cell_data[level]); + + } + + + + + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + mg_transfer.initialize_constraints(mg_constrained_dofs); + mg_transfer.build(dof_handler_displacement); + + smoother_data.resize(0, triangulation.n_global_levels() - 1); + + for (unsigned int level = 0; level < triangulation.n_global_levels(); + ++level) + { + if (level > 0) + { + smoother_data[level].smoothing_range = 15.; + smoother_data[level].degree = 10; + smoother_data[level].eig_cg_n_iterations = 10; + } + else + { + smoother_data[0].smoothing_range = 1e-3; + smoother_data[0].degree = numbers::invalid_unsigned_int; + smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m(); + } + mg_matrices[level].compute_diagonal(); + smoother_data[level].preconditioner = + mg_matrices[level].get_matrix_diagonal_inverse(); + } + + mg_smoother.initialize(mg_matrices, smoother_data); + + mg_coarse.initialize(mg_smoother); + mg::Matrix> mg_matrix(mg_matrices); + + mg_interface_matrices.resize(0, triangulation.n_global_levels() - 1); + for (unsigned int level = 0; level < triangulation.n_global_levels(); + ++level) + { + mg_interface_matrices[level].initialize(mg_matrices[level]); + } + mg::Matrix> mg_interface(mg_interface_matrices); + + Multigrid> mg( + mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother); + mg.set_edge_matrices(mg_interface, mg_interface); + mg.set_cycle(Multigrid>::v_cycle); + mg.set_minlevel(0); + PreconditionMG,MGTransferMatrixFree> + mf_gmg_preconditioner(dof_handler_displacement, mg, mg_transfer); + + +// *************TEST SOLVE************************* +// time the solve + +// output(distributed_solution, 66); + +// TimerOutput t(pcout, TimerOutput::never, TimerOutput::wall_times); +// { +// TimerOutput::Scope t_scope(t, "Solve_mfgmg"); +// elasticity_matrix_mf.initialize_dof_vector(distributed_displacement_sol); +// elasticity_matrix_mf.initialize_dof_vector(distributed_displacement_rhs); + +// locally_relevant_solution = system_rhs; +// ChangeVectorTypes::copy_from_system_to_displacement_vector(distributed_displacement_rhs,locally_relevant_solution.block(SolutionBlocks::displacement),displacement_to_system_dof_index_map); + +// pcout << "real rhs norm: " << system_rhs.block(SolutionBlocks::displacement).l2_norm() << std::endl; + +// locally_relevant_solution = distributed_solution; +// ChangeVectorTypes::copy_from_system_to_displacement_vector(distributed_displacement_sol,locally_relevant_solution.block(SolutionBlocks::displacement),displacement_to_system_dof_index_map); + +// SolverControl test_solver_control_1(20, 1e-6); +// SolverCG> CG_Solve_1(test_solver_control_1); + +// pcout << "pre norm: " << distributed_displacement_rhs.l2_norm() << std::endl; + +// try +// { +// // mf_gmg_preconditioner.vmult(distributed_displacement_sol, distributed_displacement_rhs); +// // CG_Solve.solve(elasticity_matrix_mf, distributed_displacement_sol, -1* distributed_displacement_rhs, ); +// CG_Solve_1.solve(elasticity_matrix_mf, distributed_displacement_sol, -1* distributed_displacement_rhs, mf_gmg_preconditioner); +// } +// catch(std::exception &exc) +// { +// pcout << "mfgmg diff: " << solver_control.initial_value()/solver_control.last_value() << std::endl; +// } + +// pcout << "mfgmg solved in " << test_solver_control_1.last_step() << " steps" << std::endl; +// // try +// // { +// // // CG_Solve.solve(elasticity_matrix_mf, distributed_displacement_sol, -1* distributed_displacement_rhs, PreconditionIdentity()); +// // CG_Solve_2.solve(system_matrix.block(SolutionBlocks::displacement,SolutionBlocks::displacement_multiplier), distributed_solution.block(SolutionBlocks::displacement_multiplier), system_rhs.block(SolutionBlocks::displacement), PreconditionIdentity() ); +// // } +// // catch(std::exception &exc) +// // { +// // std::cout << "solve failed in " << test_solver_control_2.last_step() << " steps" << std::endl; +// // throw; +// // } + +// // std::cout << "solved in " << test_solver_control_2.last_step() << " steps" << std::endl; + + +// } +// ChangeVectorTypes::copy_from_displacement_to_system_vector(distributed_solution.block(SolutionBlocks::displacement), distributed_displacement_sol,displacement_to_system_dof_index_map); +// displacement_constraints.distribute(distributed_solution.block(SolutionBlocks::displacement)); + +// pcout << distributed_displacement_sol.linfty_norm() << "+++++++++++++" << std::endl; + +// int a = Utilities::MPI::n_mpi_processes(mpi_communicator); + +// ChangeVectorTypes::copy(distributed_solution.block(SolutionBlocks::density),active_density_vector); + + + +// TrilinosWrappers::PreconditionAMG amg_pre; +// amg_pre.initialize(system_matrix.block(SolutionBlocks::displacement,SolutionBlocks::displacement_multiplier)); +// { +// TimerOutput::Scope t_scope(t, "Solve_AMG"); + +// SolverControl test_solver_control_2(50000, 1e-6); +// SolverCG CG_Solve_2(test_solver_control_2); + +// distributed_solution.block(SolutionBlocks::displacement_multiplier) = 0.; + +// try +// { +// CG_Solve_2.solve(system_matrix.block(SolutionBlocks::displacement,SolutionBlocks::displacement_multiplier), distributed_solution.block(SolutionBlocks::displacement_multiplier), system_rhs.block(SolutionBlocks::displacement), amg_pre); +// } +// catch(std::exception &exc) +// { +// std::cout << "solve failed in " << test_solver_control_2.last_step() << " steps" << std::endl; +// throw; +// } + +// std::cout << "amg solved in " << test_solver_control_2.last_step() << " steps" << std::endl; + +// } +// distributed_solution.block(SolutionBlocks::displacement_multiplier)=0; +// { +// TimerOutput::Scope t_scope(t, "Solve_CG"); +// SolverControl test_solver_control_3(50000, 1e-6); +// SolverCG CG_Solve_3(test_solver_control_3); +// try +// { +// CG_Solve_3.solve(system_matrix.block(SolutionBlocks::displacement,SolutionBlocks::displacement_multiplier), distributed_solution.block(SolutionBlocks::displacement), system_rhs.block(SolutionBlocks::displacement), PreconditionIdentity()); +// } +// catch(std::exception &exc) +// { +// std::cout << "solve failed in " << test_solver_control_3.last_step() << " steps" << std::endl; +// throw; +// } + +// std::cout << "CG solved in " << test_solver_control_3.last_step() << " steps" << std::endl; + +// } + + +// t.print_summary(); +// MPI_Abort(mpi_communicator, 1); + +// ***************END TEST SOLVE************************* + + + + TopOptSchurPreconditioner preconditioner(system_matrix, dof_handler, elasticity_matrix_mf, mf_gmg_preconditioner, displacement_to_system_dof_index_map); + // pcout << "about to solve" << std::endl; + // preconditioner.initialize(system_matrix, boundary_values, dof_handler, distributed_solution); + // FullMatrix out; + // out.reinit(system_matrix.m(),system_matrix.n()); + // LA::MPI::BlockVector e_j (system_rhs); + // LA::MPI::BlockVector r_j (system_rhs); + // LA::MPI::BlockVector r2_j (system_rhs); + // for (unsigned int j=0; j B_fgmres(solver_control); + // B_fgmres.solve(system_matrix, distributed_solution, system_rhs, preconditioner); + // pcout << solver_control.last_step() << " steps to solve with GMRES" << std::endl; + // break; + // } + case SolverOptions::inexact_K_with_inexact_A_gmres: { + pcout << "size of rhs block 0 : " << system_rhs.block(0).l1_norm()<< std::endl; + pcout << "size of rhs block 1 : " << system_rhs.block(1).l1_norm()<< std::endl; + pcout << "size of rhs block 2 : " << system_rhs.block(2).l1_norm()<< std::endl; + pcout << "size of rhs block 3 : " << system_rhs.block(3).l1_norm()<< std::endl; + pcout << "size of rhs block 4 : " << system_rhs.block(4).l1_norm()<< std::endl; + pcout << "size of rhs block 5 : " << system_rhs.block(5).l1_norm()<< std::endl; + pcout << "size of rhs block 6 : " << system_rhs.block(6).l1_norm()<< std::endl; + pcout << "size of rhs block 7 : " << system_rhs.block(7).l1_norm()<< std::endl; + pcout << "size of rhs block 8 : " << system_rhs.block(8).l1_norm()<< std::endl; + pcout << "size of rhs block 9 : " << system_rhs.block(9).l1_norm()<< std::endl; + + preconditioner.initialize(system_matrix, boundary_values, dof_handler, distributed_solution); + pcout << "preconditioner initialized" << std::endl; + distributed_solution = 0.; + SolverFGMRES C_fgmres(solver_control); + C_fgmres.solve(system_matrix, distributed_solution, system_rhs, preconditioner); + pcout << solver_control.last_step() << " steps to solve with FGMRES" << std::endl; + break; + } + default: + throw; + + } + + constraints.distribute(distributed_solution); + pcout << "size of distributed solution block 0 : " << distributed_solution.block(0).l1_norm()<< std::endl; + pcout << "size of distributed solution block 1 : " << distributed_solution.block(1).l1_norm()<< std::endl; + pcout << "size of distributed solution block 2 : " << distributed_solution.block(2).l1_norm()<< std::endl; + pcout << "size of distributed solution block 3 : " << distributed_solution.block(3).l1_norm()<< std::endl; + pcout << "size of distributed solution block 4 : " << distributed_solution.block(4).l1_norm()<< std::endl; + pcout << "size of distributed solution block 5 : " << distributed_solution.block(5).l1_norm()<< std::endl; + pcout << "size of distributed solution block 6 : " << distributed_solution.block(6).l1_norm()<< std::endl; + pcout << "size of distributed solution block 7 : " << distributed_solution.block(7).l1_norm()<< std::endl; + pcout << "size of distributed solution block 8 : " << distributed_solution.block(8).l1_norm()<< std::endl; + pcout << "size of distributed solution block 9 : " << distributed_solution.block(9).l1_norm()<< std::endl; + output(distributed_solution,100); + return distributed_solution; +} + +///Calculates and stores the first RHS norm for comparison with future RHS norm values +template +void +KktSystem::calculate_initial_rhs_error() { + initial_rhs_error = system_rhs.l2_norm(); +} + +///Creates an initial state vector used as an initial guess for the nonlinear solver. +template +LA::MPI::BlockVector +KktSystem::get_initial_state() { + + std::vector block_component(10, 2); + block_component[SolutionBlocks::density] = 0; + block_component[SolutionBlocks::displacement] = 1; + const std::vector dofs_per_block = + DoFTools::count_dofs_per_fe_block(dof_handler, block_component); + const unsigned int n_p = dofs_per_block[0]; + const unsigned int n_u = dofs_per_block[1]; + const std::vector block_sizes = {n_p, n_p, n_p, n_p, n_p, n_u, n_u, n_p, n_p, 1}; + + LA::MPI::BlockVector state(owned_partitioning, mpi_communicator); + { + using namespace SolutionBlocks; + state.block(density).add(density_ratio); + state.block(unfiltered_density).add(density_ratio); + state.block(unfiltered_density_multiplier) + .add(density_ratio); + state.block(density_lower_slack).add(density_ratio); + state.block(density_lower_slack_multiplier).add(50); + state.block(density_upper_slack).add(1 - density_ratio); + state.block(density_upper_slack_multiplier).add(50); + state.block(total_volume_multiplier).add(1); + state.block(displacement).add(0); + state.block(displacement_multiplier).add(0); + // state.compress(VectorOperation::add); + + // RANDOM PART HERE + // for(unsigned int k = 0; k +void +KktSystem::output(const LA::MPI::BlockVector &state, const unsigned int j) const { + locally_relevant_solution = state; + std::vector solution_names(1, "low_slack_multiplier"); + std::vector data_component_interpretation( + 1, DataComponentInterpretation::component_is_scalar); + solution_names.emplace_back("upper_slack_multiplier"); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + solution_names.emplace_back("low_slack"); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + solution_names.emplace_back("upper_slack"); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + solution_names.emplace_back("unfiltered_density"); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + for (unsigned int i = 0; i < dim; i++) { + solution_names.emplace_back("displacement"); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_part_of_vector); + } + for (unsigned int i = 0; i < dim; i++) { + solution_names.emplace_back("displacement_multiplier"); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_part_of_vector); + } + solution_names.emplace_back("density_multiplier"); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + solution_names.emplace_back("density"); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + solution_names.emplace_back("volume_multiplier"); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(locally_relevant_solution, + solution_names, + DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches(); + std::string output("solution" + std::to_string(j) + ".vtu"); + data_out.write_vtu_in_parallel(output, MPI_COMM_WORLD); + +} + + + +///Outputs to a 3d-printable file. Not yet usable in parallel. +template<> +void +KktSystem<2>::output_stl(const LA::MPI::BlockVector &state) { + double height = .25; + const int dim = 2; + std::ofstream stlfile; + stlfile.open("bridge.stl"); + stlfile << "solid bridge\n" << std::scientific; + + for (const auto &cell: dof_handler.active_cell_iterators()) { + if (state.block( + SolutionBlocks::density)[cell->get_fe().component_to_system_index(0, 0)] > 0.5) { + const Tensor<1, dim> edge_directions[2] = {cell->vertex(1) - + cell->vertex(0), + cell->vertex(2) - + cell->vertex(0)}; + const Tensor<2, dim> edge_tensor( + {{edge_directions[0][0], edge_directions[0][1]}, + {edge_directions[1][0], edge_directions[1][1]}}); + const bool is_right_handed_cell = (determinant(edge_tensor) > 0); + if (is_right_handed_cell) { + /* Write one side at z = 0. */ + stlfile << " facet normal " << 0.000000e+00 << " " + << 0.000000e+00 << " " << -1.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << cell->vertex(0)[0] << " " + << cell->vertex(0)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " vertex " << cell->vertex(2)[0] << " " + << cell->vertex(2)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " vertex " << cell->vertex(1)[0] << " " + << cell->vertex(1)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + stlfile << " facet normal " << 0.000000e+00 << " " + << 0.000000e+00 << " " << -1.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << cell->vertex(1)[0] << " " + << cell->vertex(1)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " vertex " << cell->vertex(2)[0] << " " + << cell->vertex(2)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " vertex " << cell->vertex(3)[0] << " " + << cell->vertex(3)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + /* Write one side at z = height. */ + stlfile << " facet normal " << 0.000000e+00 << " " + << 0.000000e+00 << " " << 1.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << cell->vertex(0)[0] << " " + << cell->vertex(0)[1] << " " << height << "\n"; + stlfile << " vertex " << cell->vertex(1)[0] << " " + << cell->vertex(1)[1] << " " << height << "\n"; + stlfile << " vertex " << cell->vertex(2)[0] << " " + << cell->vertex(2)[1] << " " << height << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + stlfile << " facet normal " << 0.000000e+00 << " " + << 0.000000e+00 << " " << 1.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << cell->vertex(1)[0] << " " + << cell->vertex(1)[1] << " " << height << "\n"; + stlfile << " vertex " << cell->vertex(3)[0] << " " + << cell->vertex(3)[1] << " " << height << "\n"; + stlfile << " vertex " << cell->vertex(2)[0] << " " + << cell->vertex(2)[1] << " " << height << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + } else /* The cell has a left-handed set up */ + { + /* Write one side at z = 0. */ + stlfile << " facet normal " << 0.000000e+00 << " " + << 0.000000e+00 << " " << -1.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << cell->vertex(0)[0] << " " + << cell->vertex(0)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " vertex " << cell->vertex(1)[0] << " " + << cell->vertex(1)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " vertex " << cell->vertex(2)[0] << " " + << cell->vertex(2)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + stlfile << " facet normal " << 0.000000e+00 << " " + << 0.000000e+00 << " " << -1.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << cell->vertex(1)[0] << " " + << cell->vertex(1)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " vertex " << cell->vertex(3)[0] << " " + << cell->vertex(3)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " vertex " << cell->vertex(2)[0] << " " + << cell->vertex(2)[1] << " " << 0.000000e+00 << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + /* Write one side at z = height. */ + stlfile << " facet normal " << 0.000000e+00 << " " + << 0.000000e+00 << " " << 1.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << cell->vertex(0)[0] << " " + << cell->vertex(0)[1] << " " << height << "\n"; + stlfile << " vertex " << cell->vertex(2)[0] << " " + << cell->vertex(2)[1] << " " << height << "\n"; + stlfile << " vertex " << cell->vertex(1)[0] << " " + << cell->vertex(1)[1] << " " << height << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + stlfile << " facet normal " << 0.000000e+00 << " " + << 0.000000e+00 << " " << 1.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << cell->vertex(1)[0] << " " + << cell->vertex(1)[1] << " " << height << "\n"; + stlfile << " vertex " << cell->vertex(2)[0] << " " + << cell->vertex(2)[1] << " " << height << "\n"; + stlfile << " vertex " << cell->vertex(3)[0] << " " + << cell->vertex(3)[1] << " " << height << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + } + for (unsigned int face_number = 0; + face_number < GeometryInfo::faces_per_cell; + ++face_number) { + const typename DoFHandler::face_iterator face = + cell->face(face_number); + if ((face->at_boundary()) || + (!face->at_boundary() && + (state.block( + SolutionBlocks::density)[cell->neighbor(face_number)->get_fe().component_to_system_index(0, 0)] < + 0.5))) { + const Tensor<1, dim> normal_vector = + (face->center() - cell->center()); + const double normal_norm = normal_vector.norm(); + if ((face->vertex(0)[0] - face->vertex(0)[0]) * + (face->vertex(1)[1] - face->vertex(0)[1]) * + 0.000000e+00 + + (face->vertex(0)[1] - face->vertex(0)[1]) * (0 - 0) * + normal_vector[0] + + (height - 0) * + (face->vertex(1)[0] - face->vertex(0)[0]) * + normal_vector[1] - + (face->vertex(0)[0] - face->vertex(0)[0]) * (0 - 0) * + normal_vector[1] - + (face->vertex(0)[1] - face->vertex(0)[1]) * + (face->vertex(1)[0] - face->vertex(0)[0]) * + normal_vector[0] - + (height - 0) * + (face->vertex(1)[1] - face->vertex(0)[1]) * 0 > + 0) { + stlfile << " facet normal " + << normal_vector[0] / normal_norm << " " + << normal_vector[1] / normal_norm << " " + << 0.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << face->vertex(0)[0] + << " " << face->vertex(0)[1] << " " + << 0.000000e+00 << "\n"; + stlfile << " vertex " << face->vertex(0)[0] + << " " << face->vertex(0)[1] << " " << height + << "\n"; + stlfile << " vertex " << face->vertex(1)[0] + << " " << face->vertex(1)[1] << " " + << 0.000000e+00 << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + stlfile << " facet normal " + << normal_vector[0] / normal_norm << " " + << normal_vector[1] / normal_norm << " " + << 0.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << face->vertex(0)[0] + << " " << face->vertex(0)[1] << " " << height + << "\n"; + stlfile << " vertex " << face->vertex(1)[0] + << " " << face->vertex(1)[1] << " " << height + << "\n"; + stlfile << " vertex " << face->vertex(1)[0] + << " " << face->vertex(1)[1] << " " + << 0.000000e+00 << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + } else { + stlfile << " facet normal " + << normal_vector[0] / normal_norm << " " + << normal_vector[1] / normal_norm << " " + << 0.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << face->vertex(0)[0] + << " " << face->vertex(0)[1] << " " + << 0.000000e+00 << "\n"; + stlfile << " vertex " << face->vertex(1)[0] + << " " << face->vertex(1)[1] << " " + << 0.000000e+00 << "\n"; + stlfile << " vertex " << face->vertex(0)[0] + << " " << face->vertex(0)[1] << " " << height + << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + stlfile << " facet normal " + << normal_vector[0] / normal_norm << " " + << normal_vector[1] / normal_norm << " " + << 0.000000e+00 << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << face->vertex(0)[0] + << " " << face->vertex(0)[1] << " " << height + << "\n"; + stlfile << " vertex " << face->vertex(1)[0] + << " " << face->vertex(1)[1] << " " + << 0.000000e+00 << "\n"; + stlfile << " vertex " << face->vertex(1)[0] + << " " << face->vertex(1)[1] << " " << height + << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + } + } + } + } + } + stlfile << "endsolid bridge"; +} + +///Outputs to a 3d-printable file. Not yet usable in parallel. +template<> +void +KktSystem<3>::output_stl(const LA::MPI::BlockVector &state) +{ + std::ofstream stlfile; + stlfile.open("bridge.stl"); + stlfile << "solid bridge\n" << std::scientific; + const int dim = 3; + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (state.block( + SolutionBlocks::unfiltered_density)[cell->get_fe().component_to_system_index(0, 0)] > 0.5) + { + for (const auto n : cell->face_indices()) + { + bool create_boundary = false; + if (cell->at_boundary(n)) + { + create_boundary = true; + } + else if (state.block( + SolutionBlocks::unfiltered_density)[cell->neighbor(n)->get_fe().component_to_system_index(0, 0)] <= 0.5) + { + create_boundary = true; + } + + if (create_boundary) + { + const auto face = cell->face(n); + const Tensor<1,dim> normal_vector = face->center() - + cell->center(); + double normal_norm = normal_vector.norm(); + const Tensor<1,dim> edge_vectors_1 = face->vertex(1) - face->vertex(0); + const Tensor<1,dim> edge_vectors_2 = face->vertex(2) - face->vertex(0); + + const Tensor<2, dim> edge_tensor ( + {{edge_vectors_1[0], edge_vectors_1[1],edge_vectors_1[2]}, + {edge_vectors_2[0], edge_vectors_2[1],edge_vectors_2[2]}, + {normal_vector[0], normal_vector[1], normal_vector[2]}}); + const bool is_right_handed_cell = (determinant(edge_tensor) > 0); + + if (is_right_handed_cell) + { + stlfile << " facet normal " + << normal_vector[0] / normal_norm << " " + << normal_vector[1] / normal_norm << " " + << normal_vector[2] / normal_norm << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << face->vertex(0)[0] + << " " << face->vertex(0)[1] << " " + << face->vertex(0)[2] << "\n"; + stlfile << " vertex " << face->vertex(1)[0] + << " " << face->vertex(1)[1] << " " + << face->vertex(1)[2] << "\n"; + stlfile << " vertex " << face->vertex(2)[0] + << " " << face->vertex(2)[1] << " " + << face->vertex(2)[2] << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + stlfile << " facet normal " + << normal_vector[0] / normal_norm << " " + << normal_vector[1] / normal_norm << " " + << normal_vector[2] / normal_norm << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << face->vertex(1)[0] + << " " << face->vertex(1)[1] << " " << face->vertex(1)[2] + << "\n"; + stlfile << " vertex " << face->vertex(3)[0] + << " " << face->vertex(3)[1] << " " << face->vertex(3)[2] + << "\n"; + stlfile << " vertex " << face->vertex(2)[0] + << " " << face->vertex(2)[1] << " " + << face->vertex(2)[2] << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + } + else + { + stlfile << " facet normal " + << normal_vector[0] / normal_norm << " " + << normal_vector[1] / normal_norm << " " + << normal_vector[2] / normal_norm << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << face->vertex(0)[0] + << " " << face->vertex(0)[1] << " " + << face->vertex(0)[2] << "\n"; + stlfile << " vertex " << face->vertex(2)[0] + << " " << face->vertex(2)[1] << " " + << face->vertex(2)[2] << "\n"; + stlfile << " vertex " << face->vertex(1)[0] + << " " << face->vertex(1)[1] << " " + << face->vertex(1)[2] << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + stlfile << " facet normal " + << normal_vector[0] / normal_norm << " " + << normal_vector[1] / normal_norm << " " + << normal_vector[2] / normal_norm << "\n"; + stlfile << " outer loop\n"; + stlfile << " vertex " << face->vertex(1)[0] + << " " << face->vertex(1)[1] << " " << face->vertex(1)[2] + << "\n"; + stlfile << " vertex " << face->vertex(2)[0] + << " " << face->vertex(2)[1] << " " << face->vertex(2)[2] + << "\n"; + stlfile << " vertex " << face->vertex(3)[0] + << " " << face->vertex(3)[1] << " " + << face->vertex(3)[2] << "\n"; + stlfile << " endloop\n"; + stlfile << " endfacet\n"; + } + + } + + } + } + } +} +} + +template class SAND::KktSystem<2>; +template class SAND::KktSystem<3>; diff --git a/SAND_top_opt/source/main.cc b/SAND_top_opt/source/main.cc new file mode 100644 index 00000000..17d419ae --- /dev/null +++ b/SAND_top_opt/source/main.cc @@ -0,0 +1,35 @@ +#include "../include/watchdog.h" +#include "../include/input_information.h" +#include + +///Above are fairly normal files to include. I also use the sparse direct package, which requiresBLAS/LAPACK +/// to perform a direct solve while I work on a fast iterative solver for this problem. + + +int +main(int argc, char *argv[]) { + try + { + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + SAND::NonlinearWatchdog elastic_problem; + elastic_problem.run(); + } + catch (std::exception &exc) { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" << std::endl; + std::cerr << "Exception on processing: " << std::endl << exc.what() + << std::endl << "Aborting!" << std::endl + << "----------------------------------------------------" << std::endl; + + return 1; + } + catch (...) { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" << std::endl; + std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl + << "----------------------------------------------------" << std::endl; + return 1; + } + + return 0; +} diff --git a/SAND_top_opt/source/markov_filter.cc b/SAND_top_opt/source/markov_filter.cc new file mode 100644 index 00000000..b367f200 --- /dev/null +++ b/SAND_top_opt/source/markov_filter.cc @@ -0,0 +1,55 @@ +// +// Created by justin on 2/17/21. +// +#include "../include/markov_filter.h" + +using namespace dealii; + +///Initialized the markov filter with the initial values +void +MarkovFilter::setup(const double objective_value_input, const double barrier_distance_input, + const double feasibility_input, const double barrier_value_input) { + objective_value = objective_value_input; + barrier_distance = barrier_distance_input; + feasibility = feasibility_input; + barrier_value = barrier_value_input; + + filter_barrier_function_value = feasibility + barrier_value * barrier_distance; +} + +///Adds new information to the markov filter +void +MarkovFilter::add_point(const double objective_value_input, const double barrier_distance_input, + const double feasibility_input) +{ + objective_value = objective_value_input; + barrier_distance = barrier_distance_input; + feasibility = feasibility_input; + + filter_barrier_function_value = objective_value + barrier_value * barrier_distance; +} + +///As the barrier always changes, this needs to be taken into account when accepting/rejecting a step. +/// This allows each point to be viewed in comparison to the current barrier value. +void +MarkovFilter::update_barrier_value(const double barrier_value_input) +{ + barrier_value = barrier_value_input; + filter_barrier_function_value = objective_value + barrier_value * barrier_distance; +} + +///Checks if a new point passes the filter. +bool +MarkovFilter::check_filter(const double objective_value_input, const double barrier_distance_input, + const double feasibility_input) const +{ + if ((objective_value_input + barrier_distance_input * barrier_value <= filter_barrier_function_value) || + (feasibility_input <= feasibility)) + { + return true; + } + else + { + return false; + } +} diff --git a/SAND_top_opt/source/matrix_free_elasticity.cc b/SAND_top_opt/source/matrix_free_elasticity.cc new file mode 100644 index 00000000..6da06dca --- /dev/null +++ b/SAND_top_opt/source/matrix_free_elasticity.cc @@ -0,0 +1,213 @@ +// +// Created by justin on 2/17/21. +// +#include "../include/matrix_free_elasticity.h" +#include "../include/input_information.h" +#include "../include/parameters_and_components.h" +namespace SAND { +using namespace dealii; + +///Constructor +template +MF_Elasticity_Operator::MF_Elasticity_Operator() + : MatrixFreeOperators::Base>(), + mpi_communicator(MPI_COMM_WORLD), + pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) +{ +} + +///Clears the objects, removes density information +template +void MF_Elasticity_Operator::clear() +{ + this->cell_data = nullptr; + MatrixFreeOperators::Base>::clear(); +} + +///Computes the diagonal for a preconditioner on the coarsest level. +template +void +MF_Elasticity_Operator::compute_diagonal () +{ + this->inverse_diagonal_entries. + reset(new DiagonalMatrix>()); + this->diagonal_entries. + reset(new DiagonalMatrix>()); + + dealii::LinearAlgebra::distributed::Vector &inverse_diagonal = this->inverse_diagonal_entries->get_vector(); + + dealii::LinearAlgebra::distributed::Vector &diagonal = this->diagonal_entries->get_vector(); + + this->data->initialize_dof_vector(inverse_diagonal); + this->data->initialize_dof_vector(diagonal); + unsigned int dummy = 0; + + + this->data->cell_loop (&MF_Elasticity_Operator::local_compute_diagonal, this, + diagonal, dummy); + this->data->cell_loop (&MF_Elasticity_Operator::local_compute_diagonal, this, + inverse_diagonal, dummy); + + this->set_constrained_entries_to_one(diagonal); + this->set_constrained_entries_to_one(inverse_diagonal); + + // diagonal.compress(VectorOperation::add); + + for (auto &local_element : inverse_diagonal) + { +// Assert(local_element > 0., +// ExcMessage("No diagonal entry in a positive definite operator " +// "should be zero or negative.")); + local_element = 1./local_element; + } + // inverse_diagonal.compress(VectorOperation::insert); + // diagona.print(lstd::cout); + pcout << "diag size: " << diagonal.size() << " with l2 norm " << diagonal.l2_norm() << std::endl; +} + +///Computes the diagonal value locally for a cell +template +void +MF_Elasticity_Operator +::local_compute_diagonal (const MatrixFree &data, + dealii::LinearAlgebra::distributed::Vector &dst, + const unsigned int &, + const std::pair &cell_range) const +{ + // pcout << "local_compute_diagonal" << std::endl; + FEEvaluation displacement (data, 0); + + AlignedVector> diagonal(displacement.dofs_per_cell); + + for (unsigned int cell=cell_range.first; cell cell_density = cell_data->density(cell, 0); + double penalized_density = std::pow(cell_density[0],Input::density_penalty_exponent); + + displacement.reinit(cell); + + for (unsigned int i=0; i(); + + displacement.begin_dof_values()[i] = make_vectorized_array (1.); + + displacement.evaluate (EvaluationFlags::gradients); + + for (unsigned int q=0; q > symgrad_term = penalized_density* 2.0 * Input::material_mu *displacement.get_symmetric_gradient(q); + VectorizedArray div_term = trace(displacement.get_symmetric_gradient(q)); + + for (unsigned int d = 0; d < dim; ++d) + { + symgrad_term[d][d] += penalized_density * Input::material_lambda * div_term; + } + + + displacement.submit_symmetric_gradient( symgrad_term , q); + } + + displacement.integrate (EvaluationFlags::gradients); + + diagonal[i] = displacement.begin_dof_values()[i]; + } + + for (unsigned int i=0; i +void +MF_Elasticity_Operator::local_apply( + const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair & cell_range) const +{ + // pcout << "local_apply" << std::endl; + FEEvaluation displacement(data,0); + + for (unsigned int cell=cell_range.first; cell cell_density = cell_data->density(cell, 0); + double penalized_density = std::pow(cell_density[0],Input::density_penalty_exponent); + + displacement.reinit(cell); +// displacement.read_dof_values(src); + + displacement.gather_evaluate(src, EvaluationFlags::gradients); + + for (unsigned int q = 0; q < displacement.n_q_points; ++q) + { + SymmetricTensor< 2, dim, VectorizedArray > symgrad_term = penalized_density* 2.0 * Input::material_mu *displacement.get_symmetric_gradient(q); + VectorizedArray div_term = penalized_density * Input::material_lambda * trace(displacement.get_symmetric_gradient(q)); + + for (unsigned int d = 0; d < dim; ++d) + { + symgrad_term[d][d] += div_term; + } + + + + displacement.submit_symmetric_gradient(symgrad_term, q); + } + displacement.integrate_scatter(EvaluationFlags::gradients, dst); + } + +} + +///Nothing is applied on a face on the LHS, so left blank. +template +void +MF_Elasticity_Operator +::local_apply_face(const dealii::MatrixFree &, + dealii::LinearAlgebra::distributed::Vector &, + const dealii::LinearAlgebra::distributed::Vector &, + const std::pair &) const +{ +} + +///Nothing is applied on a face on the LHS, so left blank. +template +void +MF_Elasticity_Operator +::local_apply_boundary_face(const dealii::MatrixFree &, + dealii::LinearAlgebra::distributed::Vector &, + const dealii::LinearAlgebra::distributed::Vector &, + const std::pair &) const +{ +} + +///Loops over all cells to apply the elasticity operatorto the entire LHS vector +template +void +MF_Elasticity_Operator +::apply_add (dealii::LinearAlgebra::distributed::Vector &dst, + const dealii::LinearAlgebra::distributed::Vector &src) const +{ + MatrixFreeOperators::Base>:: + data->cell_loop(&MF_Elasticity_Operator::local_apply, this, dst, src); +} + +///Sets cell data (density) to be input given. +template +void +MF_Elasticity_Operator::set_cell_data (const OperatorCellData &data) +{ + this->cell_data = &data; +} + + + +} + +template class SAND::MF_Elasticity_Operator<2,1,double>; +template class SAND::MF_Elasticity_Operator<3,1,double>; + diff --git a/SAND_top_opt/source/poly_pre.cc b/SAND_top_opt/source/poly_pre.cc new file mode 100644 index 00000000..c5dad827 --- /dev/null +++ b/SAND_top_opt/source/poly_pre.cc @@ -0,0 +1,7 @@ +#include "../include/poly_pre.h" +#include "../include/schur_preconditioner.h" + +namespace SAND +{ + +} diff --git a/SAND_top_opt/source/schur_preconditioner.cc b/SAND_top_opt/source/schur_preconditioner.cc new file mode 100644 index 00000000..75828980 --- /dev/null +++ b/SAND_top_opt/source/schur_preconditioner.cc @@ -0,0 +1,1434 @@ +// +// Created by justin on 2/17/21. +// +#include +#include +#include +#include +#include +#include +#include +#include +#include "../include/schur_preconditioner.h" +#include "../include/input_information.h" +#include "../include/sand_tools.h" +#include +#include +#include +#include +#include +#include +#include + +namespace SAND { +using MatrixType = dealii::TrilinosWrappers::SparseMatrix; +using VectorType = dealii::TrilinosWrappers::MPI::Vector; +using PayloadType = dealii::TrilinosWrappers::internal::LinearOperatorImplementation::TrilinosPayload; +using PayloadVectorType = typename PayloadType::VectorType; +using size_type = dealii::types::global_dof_index; + +namespace ChangeVectorTypes +{ +template +void copy_from_displacement_to_system_vector(LA::MPI::Vector &out, + const dealii::LinearAlgebra::distributed::Vector &in, + const std::map & displacement_to_system_dof_index_map) +{ + // dealii::LinearAlgebra::ReadWriteVector rwv( + // out.locally_owned_elements()); + // rwv.import(in, VectorOperation::insert); + for (const auto &index_pair : displacement_to_system_dof_index_map) + { + out[index_pair.second] = in[index_pair.first]; + } + // out.import(rwv, VectorOperation::insert); +} + +template +void copy_from_system_to_displacement_vector(dealii::LinearAlgebra::distributed::Vector &out, + const LA::MPI::Vector &in, + const std::map & displacement_to_system_dof_index_map) +{ + // dealii::LinearAlgebra::ReadWriteVector rwv( + // out.locally_owned_elements()); + // rwv.import(in, VectorOperation::insert); + for (const auto &index_pair : displacement_to_system_dof_index_map) + { + out[index_pair.first] = in[index_pair.second]; + } + // out.import(rwv, VectorOperation::insert); +} + +template +void copy(LA::MPI::Vector & out, + const dealii::LinearAlgebra::distributed::Vector &in) +{ + dealii::LinearAlgebra::ReadWriteVector rwv( + out.locally_owned_elements()); + rwv.import(in, VectorOperation::insert); + out.import(rwv, VectorOperation::insert); +} +template +void copy(dealii::LinearAlgebra::distributed::Vector &out, + const LA::MPI::Vector & in) +{ + dealii::LinearAlgebra::ReadWriteVector rwv; + rwv.reinit(in); + out.import(rwv, VectorOperation::insert); +} + + +} // namespace ChangeVectorTypes + + + +using namespace dealii; + + + template + PolyPreJ::PolyPreJ(const JinvMatrixPart &inner_matrix_in, const int degree_in): + degree(degree_in), + inner_matrix(inner_matrix_in) + { + } + + template + void + PolyPreJ::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const + { + LA::MPI::Vector temp_vec_1 = src; + LA::MPI::Vector temp_vec_2 = src; + for(int k=0; k + PolyPreK::PolyPreK(const KinvMatrixPart &inner_matrix_in, const int degree_in): + degree(degree_in), + inner_matrix(inner_matrix_in) + { + } + + template + void + PolyPreK::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const + { + LA::MPI::Vector temp_vec_1 = src; + LA::MPI::Vector temp_vec_2 = src; + for(int k=0; k +TopOptSchurPreconditioner::TopOptSchurPreconditioner(LA::MPI::BlockSparseMatrix &matrix_in, DoFHandler &big_dof_handler_in, MF_Elasticity_Operator &mf_elasticity_operator_in , PreconditionMG ,MGTransferMatrixFree> + &mf_gmg_preconditioner_in, std::map &displacement_to_system_dof_index_map) + : + system_matrix(matrix_in), + mpi_communicator(MPI_COMM_WORLD), + n_rows(0), + n_columns(0), + n_block_rows(0), + n_block_columns(0), + other_solver_control(1000000, 1e-6), + other_bicgstab(other_solver_control), + other_gmres(other_solver_control), + other_cg(other_solver_control), + a_mat(matrix_in.block(SolutionBlocks::displacement, SolutionBlocks::displacement_multiplier)), + b_mat(matrix_in.block(SolutionBlocks::density, SolutionBlocks::density)), + c_mat(matrix_in.block(SolutionBlocks::displacement,SolutionBlocks::density)), + e_mat(matrix_in.block(SolutionBlocks::displacement_multiplier,SolutionBlocks::density)), + f_mat(matrix_in.block(SolutionBlocks::unfiltered_density_multiplier,SolutionBlocks::unfiltered_density)), + f_t_mat(matrix_in.block(SolutionBlocks::unfiltered_density,SolutionBlocks::unfiltered_density_multiplier)), + d_m_mat(matrix_in.block(SolutionBlocks::density_upper_slack_multiplier, SolutionBlocks::density_upper_slack)), + d_1_mat(matrix_in.block(SolutionBlocks::density_lower_slack, SolutionBlocks::density_lower_slack)), + d_2_mat(matrix_in.block(SolutionBlocks::density_upper_slack, SolutionBlocks::density_upper_slack)), + m_vect(matrix_in.block(SolutionBlocks::density, SolutionBlocks::total_volume_multiplier)), + solver_type("Amesos_Klu"), + additional_data(false, solver_type), + direct_solver_control(1, 0), + a_inv_direct(direct_solver_control, additional_data), + a_inv_mf_gmg(mf_elasticity_operator_in, mf_gmg_preconditioner_in, a_mat, displacement_to_system_dof_index_map), + pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), + timer(pcout, TimerOutput::summary, TimerOutput::wall_times), + g_mat(f_mat, f_t_mat, d_8_mat), + h_mat(a_mat, b_mat, c_mat, e_mat, pre_amg, a_inv_direct, a_inv_mf_gmg), + j_inv_mat(h_mat, g_mat, d_m_mat), + k_inv_mat(h_mat, g_mat, d_m_mat), + j_inv_part(h_mat, g_mat, d_m_mat), + k_inv_part(h_mat, g_mat, d_m_mat), + big_dof_handler(big_dof_handler_in) +{ + + +} + +///Initializes the preconditioner with information about the boundary values and DoF Handler. +template +void TopOptSchurPreconditioner::initialize(LA::MPI::BlockSparseMatrix &matrix, const std::map &boundary_values,const DoFHandler &dof_handler, const LA::MPI::BlockVector &distributed_state) +{ + + a_inv_mf_gmg.set_exemplar_vector(distributed_state.block(SolutionBlocks::displacement)); + + TimerOutput::Scope t(timer, "initialize"); + { + TimerOutput::Scope t(timer, "diag stuff"); + const types::global_dof_index disp_start_index = system_matrix.get_row_indices().block_start( + SolutionBlocks::displacement); + const types::global_dof_index disp_mult_start_index = system_matrix.get_row_indices().block_start( + SolutionBlocks::displacement_multiplier); + + for (auto &pair: boundary_values) { + const auto dof_index=pair.first; + const types::global_dof_index n_u = system_matrix.block(SolutionBlocks::displacement, + SolutionBlocks::displacement).m(); + if ((dof_index >= disp_start_index) && (dof_index < disp_start_index + n_u)) { + double diag_val = system_matrix.block(SolutionBlocks::displacement, + SolutionBlocks::displacement).el( + dof_index - disp_start_index, dof_index - disp_start_index); + + system_matrix.block(SolutionBlocks::displacement, SolutionBlocks::displacement_multiplier).set( + dof_index - disp_start_index, dof_index - disp_start_index, diag_val); + + } + else if ((dof_index >= disp_mult_start_index) && (dof_index < disp_mult_start_index + n_u)) + { + double diag_val = system_matrix.block(SolutionBlocks::displacement_multiplier, + SolutionBlocks::displacement_multiplier).el( + dof_index - disp_mult_start_index, dof_index - disp_mult_start_index); + system_matrix.block(SolutionBlocks::displacement_multiplier, SolutionBlocks::displacement).set( + dof_index - disp_mult_start_index, dof_index - disp_mult_start_index, diag_val); + } + } + + + //set diagonal to 0? + for (auto &pair: boundary_values) { + const auto dof_index=pair.first; + const types::global_dof_index disp_start_index = system_matrix.get_row_indices().block_start( + SolutionBlocks::displacement); + const types::global_dof_index disp_mult_start_index = system_matrix.get_row_indices().block_start( + SolutionBlocks::displacement_multiplier); + const types::global_dof_index n_u = system_matrix.block(SolutionBlocks::displacement, + SolutionBlocks::displacement).m(); + if ((dof_index >= disp_start_index) && (dof_index < disp_start_index + n_u)) { + system_matrix.block(SolutionBlocks::displacement, SolutionBlocks::displacement).set( + dof_index - disp_start_index, dof_index - disp_start_index, 0); + } else if ((dof_index >= disp_mult_start_index) && (dof_index < disp_mult_start_index + n_u)) { + system_matrix.block(SolutionBlocks::displacement_multiplier, + SolutionBlocks::displacement_multiplier).set( + dof_index - disp_mult_start_index, dof_index - disp_mult_start_index, 0); + } + } + + + system_matrix.compress(VectorOperation::insert); + } + + { + TimerOutput::Scope t(timer, "reinit diag matrices"); + d_3_mat.reinit(matrix.block(SolutionBlocks::density, SolutionBlocks::density)); + d_4_mat.reinit(matrix.block(SolutionBlocks::density, SolutionBlocks::density)); + d_5_mat.reinit(matrix.block(SolutionBlocks::density_lower_slack_multiplier, SolutionBlocks::density_lower_slack_multiplier)); + d_6_mat.reinit(matrix.block(SolutionBlocks::density, SolutionBlocks::density)); + d_7_mat.reinit(matrix.block(SolutionBlocks::density, SolutionBlocks::density)); + d_8_mat.reinit(matrix.block(SolutionBlocks::density, SolutionBlocks::density)); + d_m_inv_mat.reinit(matrix.block(SolutionBlocks::density, SolutionBlocks::density)); + + + d_3_mat=0; + d_4_mat=0; + d_5_mat=0; + d_6_mat=0; + d_7_mat=0; + d_8_mat=0; + d_m_inv_mat=0; + } + { + const types::global_dof_index n_p = system_matrix.block(SolutionBlocks::density, + SolutionBlocks::density).m(); + + std::vector l_global(n_p); + std::vector lm_global(n_p); + std::vector u_global(n_p); + std::vector um_global(n_p); + std::vector m_global(n_p); + + std::vector l(n_p); + std::vector lm(n_p); + std::vector u(n_p); + std::vector um(n_p); + std::vector m(n_p); + + TimerOutput::Scope t(timer, "build diag matrices"); + for (const auto cell: dof_handler.active_cell_iterators()) + { + if(cell->is_locally_owned()) + { + std::vector i(cell->get_fe().n_dofs_per_cell()); + cell->get_dof_indices(i); + + const int i_ind = cell->get_fe().component_to_system_index(0,0); + + if(distributed_state.block(SolutionBlocks::density_lower_slack_multiplier).in_local_range(i[i_ind])) + { + lm[i[i_ind]] = distributed_state.block(SolutionBlocks::density_lower_slack_multiplier)[i[i_ind]]; + m[i[i_ind]] = cell->measure(); + } + + if(distributed_state.block(SolutionBlocks::density_lower_slack).in_local_range(i[i_ind])) + { + l[i[i_ind]] = distributed_state.block(SolutionBlocks::density_lower_slack)[i[i_ind]]; + } + + if(distributed_state.block(SolutionBlocks::density_upper_slack_multiplier).in_local_range(i[i_ind])) + { + um[i[i_ind]]= distributed_state.block(SolutionBlocks::density_upper_slack_multiplier)[i[i_ind]]; + } + + if(distributed_state.block(SolutionBlocks::density_upper_slack).in_local_range(i[i_ind])) + { + u[i[i_ind]] = distributed_state.block(SolutionBlocks::density_upper_slack)[i[i_ind]]; + } + } + } + + MPI_Allreduce(lm.data(), lm_global.data(), n_p, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(l.data(), l_global.data(), n_p, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(um.data(), um_global.data(), n_p, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(u.data(), u_global.data(), n_p, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(m.data(), m_global.data(), n_p, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + + for (unsigned int k=0; k< n_p; k++) + { + if(distributed_state.block(0).in_local_range(k)) + { + d_3_mat.set(k, k, -1 * lm_global[k]/(m_global[k]*l_global[k])); + d_4_mat.set(k, k, -1 * um_global[k]/(m_global[k]*u_global[k])); + d_5_mat.set(k, k, lm_global[k]/l_global[k]); + d_6_mat.set(k, k, um_global[k]/u_global[k]); + d_7_mat.set(k, k, m_global[k]*(lm_global[k]*u_global[k] + um_global[k]*l_global[k])/(l_global[k]*u_global[k])); + d_8_mat.set(k, k, l_global[k]*u_global[k]/(m_global[k]*(lm_global[k]*u_global[k] + um_global[k]*l_global[k]))); + d_m_inv_mat.set(k, k, 1 / m_global[k]); + } + } + } + d_3_mat.compress(VectorOperation::insert); + d_4_mat.compress(VectorOperation::insert); + d_5_mat.compress(VectorOperation::insert); + d_6_mat.compress(VectorOperation::insert); + d_7_mat.compress(VectorOperation::insert); + d_8_mat.compress(VectorOperation::insert); + d_m_inv_mat.compress(VectorOperation::insert); + + pre_j=distributed_state.block(SolutionBlocks::density); + pre_k=distributed_state.block(SolutionBlocks::density); + g_d_m_inv_density=distributed_state.block(SolutionBlocks::density); + k_g_d_m_inv_density=distributed_state.block(SolutionBlocks::density); + + LA::MPI::Vector density_exemplar = distributed_state.block(SolutionBlocks::density); + LA::MPI::Vector displacement_exemplar = distributed_state.block(SolutionBlocks::displacement); + + density_exemplar = 0.; + displacement_exemplar = 0.; + g_mat.initialize(density_exemplar, d_m_inv_mat); + h_mat.initialize(density_exemplar, displacement_exemplar,d_m_inv_mat); + j_inv_mat.initialize(density_exemplar,d_m_inv_mat); + k_inv_mat.initialize(density_exemplar,d_m_inv_mat); + // j_inv_part.initialize(density_exemplar); + // k_inv_part.initialize(density_exemplar); + + num_mults = 0; +} + +///Application of the preconditioner, broken up into 5 parts. +template +void TopOptSchurPreconditioner::vmult(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const { + + LA::MPI::BlockVector dst_tmp = dst; + dst_tmp = 0.; + LA::MPI::BlockVector src_temp = src; + + // pcout << "size of matrix block: " << system_matrix.block(9,8).linfty_norm() << " " << system_matrix.block(5,6).linfty_norm() << std::endl; + LA::MPI::BlockVector temp_src; + { + dst = 0.; + TimerOutput::Scope t(timer, "part 1"); + vmult_step_1(dst, src); + temp_src = dst; + pcout << "step 1 done: " << dst.l2_norm() << std::endl; + } + + { + dst = 0.; + TimerOutput::Scope t(timer, "part 2"); + vmult_step_2(dst, temp_src); + temp_src = dst; + pcout << "step 2 done: " << dst.l2_norm() << std::endl; + } + + { + dst = 0.; + TimerOutput::Scope t(timer, "part 3"); + vmult_step_3(dst, temp_src); + temp_src = dst; + pcout << "step 3 done: " << dst.l2_norm() << std::endl; + } + + { + dst = 0.; + TimerOutput::Scope t(timer, "part 4"); + vmult_step_4(dst, temp_src); + temp_src = dst; + pcout << "step 4 done: " << dst.l2_norm() << std::endl; + } + + dst = 0.; + vmult_step_5(dst, temp_src); + pcout << "step 5 done: " << dst.l2_norm() << std::endl; + num_mults++; + +} + +///Not implemented +template +void TopOptSchurPreconditioner::Tvmult(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const { + dst = src; +} + +template +void TopOptSchurPreconditioner::vmult_add(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const { + LA::MPI::BlockVector dst_temp = dst; + vmult(dst_temp, src); + dst += dst_temp; +} + +///Not implemented +template +void TopOptSchurPreconditioner::Tvmult_add(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const { + dst = src; +} + +template +void TopOptSchurPreconditioner::vmult_step_1(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const { + dst = src; + auto dst_temp = dst; + auto dst_temp2 = dst; + auto dst_temp3 = dst; + + d_5_mat.vmult(dst_temp2.block(SolutionBlocks::density_lower_slack_multiplier),src.block(SolutionBlocks::density_lower_slack_multiplier)); + d_6_mat.vmult(dst_temp3.block(SolutionBlocks::density_upper_slack_multiplier),src.block(SolutionBlocks::density_upper_slack_multiplier)); + + dst.block(SolutionBlocks::unfiltered_density) = dst_temp.block(SolutionBlocks::unfiltered_density) + - dst_temp2.block(SolutionBlocks::density_lower_slack_multiplier) + + dst_temp3.block(SolutionBlocks::density_upper_slack_multiplier) + + src.block(SolutionBlocks::density_lower_slack) + - src.block(SolutionBlocks::density_upper_slack); +} + +template +void TopOptSchurPreconditioner::vmult_step_2(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const { + dst = src; + auto dst_temp = dst; + auto dst_temp2 = dst; + auto dst_temp3 = dst; + d_8_mat.vmult(dst_temp2.block(SolutionBlocks::unfiltered_density), src.block(SolutionBlocks::unfiltered_density)); + f_mat.vmult(dst_temp3.block(SolutionBlocks::unfiltered_density),dst_temp2.block(SolutionBlocks::unfiltered_density)); + dst.block(SolutionBlocks::unfiltered_density_multiplier) = dst_temp.block(SolutionBlocks::unfiltered_density_multiplier) - dst_temp3.block(SolutionBlocks::unfiltered_density); + +} + +template +void TopOptSchurPreconditioner::vmult_step_3(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const { + + dst = src; + auto dst_temp = dst; + + if (Input::solver_choice==SolverOptions::inexact_K_with_inexact_A_gmres) + { + + a_inv_mf_gmg.set_tol(0); + a_inv_mf_gmg.set_iter(100); + + a_inv_mf_gmg.vmult(dst_temp.block(SolutionBlocks::displacement_multiplier),src.block(SolutionBlocks::displacement_multiplier)); + a_inv_mf_gmg.vmult(dst_temp.block(SolutionBlocks::displacement),src.block(SolutionBlocks::displacement)); + + a_inv_mf_gmg.set_tol(Input::a_rel_tol); + a_inv_mf_gmg.set_iter(Input::a_inv_iterations); + + c_mat.Tvmult(dst_temp.block(SolutionBlocks::density_upper_slack),dst_temp.block(SolutionBlocks::displacement_multiplier)); + e_mat.Tvmult(dst_temp.block(SolutionBlocks::density_lower_slack),dst_temp.block(SolutionBlocks::displacement)); + + dst.block(SolutionBlocks::density) = dst_temp.block(SolutionBlocks::density) + dst_temp.block(SolutionBlocks::density_upper_slack) + dst_temp.block(SolutionBlocks::density_lower_slack); + + } + else + { + a_inv_direct.vmult(dst_temp.block(SolutionBlocks::displacement_multiplier),src.block(SolutionBlocks::displacement_multiplier)); + a_inv_direct.vmult(dst_temp.block(SolutionBlocks::displacement),src.block(SolutionBlocks::displacement)); + c_mat.Tvmult(dst_temp.block(SolutionBlocks::density_upper_slack),dst_temp.block(SolutionBlocks::displacement_multiplier)); + e_mat.Tvmult(dst_temp.block(SolutionBlocks::density_lower_slack),dst_temp.block(SolutionBlocks::displacement)); + dst.block(SolutionBlocks::density) = dst_temp.block(SolutionBlocks::density) - dst_temp.block(SolutionBlocks::density_upper_slack) - dst_temp.block(SolutionBlocks::density_lower_slack); + } + +} + +template +void TopOptSchurPreconditioner::vmult_step_4(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const { + dst = src; + auto dst_temp = dst; + auto k_density_mult = src.block(SolutionBlocks::density); + + + TrilinosWrappers::PreconditionIdentity preconditioner; + preconditioner.initialize(b_mat); + + // PolyPreK poly_pre_k(k_inv_part,0); + // PolyPreJ poly_pre_j(j_inv_part,0); + + auto d_m_inv_density = g_d_m_inv_density; + d_m_inv_mat.vmult(d_m_inv_density,src.block(SolutionBlocks::density)); + g_mat.vmult(g_d_m_inv_density,d_m_inv_density); + + SolverControl step_4_gmres_control_1 (Input::k_inv_iterations, g_d_m_inv_density.l2_norm()*Input::k_rel_tol); + SolverFGMRES step_4_gmres_1 (step_4_gmres_control_1); + try { + pcout << "SOLVE 1" << std::endl; + k_density_mult = 0.; + step_4_gmres_1.solve(k_inv_mat,k_g_d_m_inv_density,g_d_m_inv_density, PreconditionIdentity() ); + } + catch (std::exception &exc) + { + pcout << "Failure of linear solver 4-1 with convergence of " << step_4_gmres_control_1.last_value()/step_4_gmres_control_1.initial_value()<< std::endl; + // throw; + } + SolverControl step_4_gmres_control_2 (Input::k_inv_iterations, src.block(SolutionBlocks::unfiltered_density_multiplier).l2_norm()*Input::k_rel_tol); + SolverFGMRES step_4_gmres_2 (step_4_gmres_control_2); + try { + pcout << "SOLVE 2" << std::endl; + k_density_mult = 0.; + step_4_gmres_2.solve(k_inv_mat,k_density_mult,src.block(SolutionBlocks::unfiltered_density_multiplier), PreconditionIdentity()); + } catch (std::exception &exc) + { + + pcout << "Failure of linear solver 4-2 with convergence of " << step_4_gmres_control_2.last_value()/step_4_gmres_control_2.initial_value()<< std::endl; + // throw; + } + + if (num_mults == -10) + { + TimerOutput::Scope t(timer, "k_mult"); + SolverControl step_5_gmres_control_4 (100000, 1e-12*g_d_m_inv_density.l2_norm()); + step_5_gmres_control_4.enable_history_data(); + SolverFGMRES step_5_gmres_4(step_5_gmres_control_4); + try { + TimerOutput::Scope t(timer, "actual inverse 5.3"); + pcout << "SOLVE TEST" << std::endl; + dst_temp.block(SolutionBlocks::density) = 0.; + step_5_gmres_4.solve(k_inv_mat,dst_temp.block(SolutionBlocks::density), g_d_m_inv_density , PreconditionIdentity()); + } catch (std::exception &exc) + { + pcout << "Failure of linear solver step_5_gmres_4" << std::endl; + pcout << "first residual: " << step_5_gmres_control_4.initial_value() << std::endl; + pcout << "last residual: " << step_5_gmres_control_4.last_value() << std::endl; + // throw; + } + auto history = step_5_gmres_control_4.get_history_data(); + for (int i=0; i< step_5_gmres_control_4.last_step()+1; ++i) + { + pcout << i << " " << history[i]/step_5_gmres_control_4.initial_value() << std::endl; + } + + + } + dst.block(SolutionBlocks::total_volume_multiplier) = transpose_operator(m_vect)*k_g_d_m_inv_density + - transpose_operator(m_vect)*k_density_mult + +dst_temp.block(SolutionBlocks::total_volume_multiplier); + +} + +///The only non-triangular vmult step, applies inverses down the block diagonal. +template +void TopOptSchurPreconditioner::vmult_step_5(LA::MPI::BlockVector &dst, const LA::MPI::BlockVector &src) const +{ + { + //First Block Inverse + TimerOutput::Scope t(timer, "inverse 1"); + dst.block(SolutionBlocks::density_lower_slack_multiplier) = linear_operator(d_3_mat) * src.block(SolutionBlocks::density_lower_slack_multiplier) + + linear_operator(d_m_inv_mat) * src.block(SolutionBlocks::density_lower_slack); + dst.block(SolutionBlocks::density_upper_slack_multiplier) = linear_operator(d_4_mat) * src.block(SolutionBlocks::density_upper_slack_multiplier) + + linear_operator(d_m_inv_mat) * src.block(SolutionBlocks::density_upper_slack); + dst.block(SolutionBlocks::density_lower_slack) = linear_operator(d_m_inv_mat) * src.block(SolutionBlocks::density_lower_slack_multiplier); + dst.block(SolutionBlocks::density_upper_slack) = linear_operator(d_m_inv_mat) * src.block(SolutionBlocks::density_upper_slack_multiplier); + } + + { + //Second Block Inverse + TimerOutput::Scope t(timer, "inverse 2"); + dst.block(SolutionBlocks::unfiltered_density) = + linear_operator(d_8_mat) * src.block(SolutionBlocks::unfiltered_density); + } + + + { + //Third Block Inverse + TimerOutput::Scope t(timer, "inverse 3"); + if (Input::solver_choice==SolverOptions::inexact_K_with_inexact_A_gmres) + { + LA::MPI::BlockVector dst_temp = dst; + + // SolverControl solver_control(100000, 1e-6); + // SolverCG a_solver_cg(solver_control); + + // auto a_inv_op = inverse_operator(linear_operator(a_mat),a_solver_cg, pre_amg); + + // dst.block(SolutionBlocks::displacement) = a_inv_op * src.block(SolutionBlocks::displacement_multiplier); + // dst.block(SolutionBlocks::displacement_multiplier) = a_inv_op * src.block(SolutionBlocks::displacement); + + a_inv_mf_gmg.set_tol(0); + a_inv_mf_gmg.set_iter(100); + a_inv_mf_gmg.vmult( dst_temp.block(SolutionBlocks::displacement), src.block(SolutionBlocks::displacement_multiplier)); + a_inv_mf_gmg.vmult( dst_temp.block(SolutionBlocks::displacement_multiplier), src.block(SolutionBlocks::displacement)); + a_inv_mf_gmg.set_tol(Input::a_rel_tol); + a_inv_mf_gmg.set_iter(Input::a_inv_iterations); + + dst.block(SolutionBlocks::displacement) = -1 * dst_temp.block(SolutionBlocks::displacement); + dst.block(SolutionBlocks::displacement_multiplier) = -1 * dst_temp.block(SolutionBlocks::displacement_multiplier); + } + else + { + a_inv_direct.vmult( dst.block(SolutionBlocks::displacement), src.block(SolutionBlocks::displacement_multiplier)); + a_inv_direct.vmult( dst.block(SolutionBlocks::displacement_multiplier), src.block(SolutionBlocks::displacement)); + } + + + } + { + //Fourth (ugly) Block Inverse + TimerOutput::Scope t(timer, "inverse 4"); + + + + if (Input::solver_choice == SolverOptions::exact_preconditioner_with_gmres) + { + // pre_j = src.block(SolutionBlocks::density) + linear_operator(h_mat) * linear_operator(d_m_inv_mat) * src.block(SolutionBlocks::unfiltered_density_multiplier); + // pre_k = -1* linear_operator(g_mat) * linear_operator(d_m_inv_mat) * src.block(SolutionBlocks::density) + src.block(SolutionBlocks::unfiltered_density_multiplier); + // dst.block(SolutionBlocks::unfiltered_density_multiplier) = transpose_operator(k_mat) * pre_j; + // dst.block(SolutionBlocks::density) = linear_operator(k_mat) * pre_k; + } + + else if (Input::solver_choice == SolverOptions::inexact_K_with_inexact_A_gmres) + { + + // TrilinosWrappers::PreconditionIdentity preconditioner; + // preconditioner.initialize(b_mat); + auto pre_pre_k = pre_k; + auto pre_pre_pre_k = pre_k; + auto d_m_inv_unfil_density_mult = pre_k; + auto h_d_m_inv_unfil_density_mult = pre_k; + + { + TimerOutput::Scope t(timer, "not inverse 5.1"); + d_m_inv_mat.vmult(d_m_inv_unfil_density_mult, src.block(SolutionBlocks::unfiltered_density_multiplier)); + h_mat.vmult(h_d_m_inv_unfil_density_mult,d_m_inv_unfil_density_mult); + pre_j = src.block(SolutionBlocks::density) + h_d_m_inv_unfil_density_mult; + d_m_inv_mat.vmult(pre_pre_pre_k,src.block(SolutionBlocks::density)); + g_mat.vmult(pre_pre_k,pre_pre_pre_k); + pre_k = -1 * pre_pre_k + src.block(SolutionBlocks::unfiltered_density_multiplier); + } + SolverControl step_5_gmres_control_1 (Input::k_inv_iterations, Input::k_rel_tol*pre_j.l2_norm()); + SolverFGMRES step_5_gmres_1 (step_5_gmres_control_1); + try { + TimerOutput::Scope t(timer, "actual inverse 5.1"); + dst.block(SolutionBlocks::unfiltered_density_multiplier) = 0.; + step_5_gmres_1.solve(j_inv_mat, dst.block(SolutionBlocks::unfiltered_density_multiplier), pre_j , PreconditionIdentity()); + } catch (std::exception &exc) + { + pcout << "Failure of linear solver 5-1 with convergence of " << step_5_gmres_control_1.last_value()/step_5_gmres_control_1.initial_value()<< std::endl; + // throw; + } + + + SolverControl step_5_gmres_control_2 (Input::k_inv_iterations, Input::k_rel_tol*pre_k.l2_norm()); + SolverFGMRES step_5_gmres_2 (step_5_gmres_control_2); + try { + TimerOutput::Scope t(timer, "actual inverse 5.2"); + dst.block(SolutionBlocks::density) = 0.; + step_5_gmres_2.solve(k_inv_mat,dst.block(SolutionBlocks::density), pre_k , PreconditionIdentity()); + } catch (std::exception &exc) + { + pcout << "Failure of linear solver 5-2 with convergence of " << step_5_gmres_control_2.last_value()/step_5_gmres_control_2.initial_value()<< std::endl; + // throw; + } + } + else if (Input::solver_choice == SolverOptions::inexact_K_with_exact_A_gmres) + { + + // TrilinosWrappers::PreconditionIdentity preconditioner; + // preconditioner.initialize(b_mat); + auto pre_pre_k = pre_k; + auto pre_pre_pre_k = pre_k; + auto d_m_inv_unfil_density_mult = pre_k; + auto h_d_m_inv_unfil_density_mult = pre_k; + + { + TimerOutput::Scope t(timer, "not inverse 5.1"); + d_m_inv_mat.vmult(d_m_inv_unfil_density_mult, src.block(SolutionBlocks::unfiltered_density_multiplier)); + h_mat.vmult(h_d_m_inv_unfil_density_mult,d_m_inv_unfil_density_mult); + pre_j = src.block(SolutionBlocks::density) + h_d_m_inv_unfil_density_mult; + d_m_inv_mat.vmult(pre_pre_pre_k,src.block(SolutionBlocks::density)); + g_mat.vmult(pre_pre_k,pre_pre_pre_k); + pre_k = -1 * pre_pre_k + src.block(SolutionBlocks::unfiltered_density_multiplier); + } + SolverControl step_5_gmres_control_1 (Input::k_inv_iterations, Input::k_rel_tol*pre_j.l2_norm()); + SolverFGMRES step_5_gmres_1 (step_5_gmres_control_1); + try { + TimerOutput::Scope t(timer, "actual inverse 5.1"); + dst.block(SolutionBlocks::unfiltered_density_multiplier) = 0.; + pcout << "SOLVE 3" << std::endl; + step_5_gmres_1.solve(j_inv_mat, dst.block(SolutionBlocks::unfiltered_density_multiplier), pre_j , PreconditionIdentity()); + } catch (std::exception &exc) + { + pcout << "Failure of linear solver step_5_gmres_1" << std::endl; + pcout << "first residual: " << step_5_gmres_control_1.initial_value() << std::endl; + pcout << "last residual: " << step_5_gmres_control_1.last_value() << std::endl; + // throw; + } + + + SolverControl step_5_gmres_control_2 (Input::k_inv_iterations, Input::k_rel_tol*pre_k.l2_norm()); + SolverFGMRES step_5_gmres_2 (step_5_gmres_control_2); + try { + TimerOutput::Scope t(timer, "actual inverse 5.2"); + dst.block(SolutionBlocks::density) = 0.; + pcout << "SOLVE 4" << std::endl; + step_5_gmres_2.solve(k_inv_mat,dst.block(SolutionBlocks::density), pre_k , PreconditionIdentity()); + } catch (std::exception &exc) + { + pcout << "Failure of linear solver step_5_gmres_2" << std::endl; + pcout << "first residual: " << step_5_gmres_control_2.initial_value() << std::endl; + pcout << "last residual: " << step_5_gmres_control_2.last_value() << std::endl; + // throw; + } + } + else + { + pcout << "shouldn't get here"; + throw; + } + + } + { + + dst.block(SolutionBlocks::total_volume_multiplier) = src.block(SolutionBlocks::total_volume_multiplier); + + + } + +} + +///I used to use this function to output parts of the preconditioner for debugging. +template +void TopOptSchurPreconditioner::print_stuff() +{ + + +} + +//******************************* Direct Trilinos Solver *********************************** + +///Constructor +VmultTrilinosSolverDirect::VmultTrilinosSolverDirect(SolverControl &cn, + const TrilinosWrappers::SolverDirect::AdditionalData &data) + : solver_direct(cn, data) +{ + +} + +///Initialize a direct solver - works well up to a point of refinement, and then refuses to solve. +void VmultTrilinosSolverDirect::initialize(LA::MPI::SparseMatrix &a_mat) +{ + reinit(a_mat); + solver_direct.initialize(a_mat); + size = a_mat.n(); +} + +///rephrases the solve as a vmult for easier use. +void +VmultTrilinosSolverDirect::vmult(LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const +{ + solver_direct.solve(dst, src); +} + +///rephrases the solve as a vmult for easier use. +void VmultTrilinosSolverDirect::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + solver_direct.solve(dst, src); +} + +///rephrases the solve as a vmult for easier use - note this is a symmetric matrix +void +VmultTrilinosSolverDirect::Tvmult(LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const +{ + solver_direct.solve(dst, src); +} + +///rephrases the solve as a vmult for easier use - note this is a symmetric matrix. +void VmultTrilinosSolverDirect::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + solver_direct.solve(dst, src); +} + + + +// **************** A inv MF GMG ********************** + +///Initializes an object that performs the A inverse operation using matrix free GMG +template +AInvMatMFGMG::AInvMatMFGMG(MF_Elasticity_Operator &mf_elasticity_operator_in , PreconditionMG, MGTransferMatrixFree > &mf_gmg_preconditioner_in, LA::MPI::SparseMatrix &a_mat, std::map &displacement_to_system_dof_index_map_in) + : mf_elasticity_operator(mf_elasticity_operator_in), + mf_gmg_preconditioner(mf_gmg_preconditioner_in), + a_mat_wrapped(a_mat), + displacement_to_system_dof_index_map(displacement_to_system_dof_index_map_in) + +{ + mf_elasticity_operator.initialize_dof_vector(temp_dst); + mf_elasticity_operator.initialize_dof_vector(temp_src); +} + +///Performs the A inverse operation +template +void AInvMatMFGMG::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + + SolverControl solver_control(iterations, src.l2_norm()*tolerance); + SolverCG> a_solver_cg(solver_control); + + ChangeVectorTypes::copy_from_system_to_displacement_vector(temp_src, src, displacement_to_system_dof_index_map); + try{ + temp_dst = 0.; + a_solver_cg.solve(mf_elasticity_operator,temp_dst,temp_src, mf_gmg_preconditioner); + } + catch (std::exception &exc) + { + // std::cout << "failed with a reduction of: " << solver_control.initial_value()/solver_control.last_value() << std::endl; + } + + ChangeVectorTypes::copy_from_displacement_to_system_vector(dst,temp_dst, displacement_to_system_dof_index_map); + dst.compress(VectorOperation::insert); + +} + +///Performs the A inverse operation - note A is symmetric +template +void AInvMatMFGMG::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + SolverControl solver_control(iterations, src.l2_norm()*tolerance); + SolverCG> a_solver_cg(solver_control); + + ChangeVectorTypes::copy_from_system_to_displacement_vector(temp_src, src, displacement_to_system_dof_index_map); + temp_dst = 0.; + a_solver_cg.solve(mf_elasticity_operator,temp_dst,temp_src, mf_gmg_preconditioner); + ChangeVectorTypes::copy_from_displacement_to_system_vector(dst,temp_dst, displacement_to_system_dof_index_map); +} +///Change the solver tolerance if needed +template +void AInvMatMFGMG::set_tol(double tol_in) +{ + tolerance = tol_in; +} + +///Change the maximum number of iterations if needed. +template +void AInvMatMFGMG::set_iter(unsigned int iterations_in) +{ + iterations = iterations_in; +} + + +// ****************** GMATRIX *********************** + +GMatrix::GMatrix(const LA::MPI::SparseMatrix &f_mat_in, const LA::MPI::SparseMatrix &f_t_mat_in, LA::MPI::SparseMatrix &d_8_mat_in) + : + f_mat(f_mat_in), + f_t_mat(f_t_mat_in), + d_8_mat(d_8_mat_in) +{ + +} + +void +GMatrix::initialize(LA::MPI::Vector &exemplar_density_vector, LA::MPI::SparseMatrix &d_m_inv_mat_in) +{ + temp_vect_1 = exemplar_density_vector; + temp_vect_2 = exemplar_density_vector; + temp_vect_1 = 0.; + temp_vect_2 = 0.; + + d_m_inv_mat.copy_from(d_m_inv_mat_in); +} + + +void GMatrix::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + temp_vect_1 = 0.; + temp_vect_2 = 0.; + dst = 0.; + f_t_mat.vmult(temp_vect_1,src); + d_8_mat.vmult(temp_vect_2, temp_vect_1); + f_mat.vmult(dst,temp_vect_2); +} + + + +void GMatrix::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + temp_vect_1 = 0.; + temp_vect_2 = 0.; + dst = 0.; + f_t_mat.Tvmult(temp_vect_1,src); + d_8_mat.vmult(temp_vect_2, temp_vect_1); + f_mat.vmult(dst,temp_vect_2); +} + + +// ****************** HMatrix *********************** + +template +HMatrix::HMatrix(LA::MPI::SparseMatrix &a_mat_in, const LA::MPI::SparseMatrix &b_mat_in, const LA::MPI::SparseMatrix &c_mat_in, const LA::MPI::SparseMatrix &e_mat_in,TrilinosWrappers::PreconditionAMG &pre_amg_in,VmultTrilinosSolverDirect &a_inv_direct_in, AInvMatMFGMG &a_inv_mf_gmg_in) + : + a_mat(a_mat_in), + b_mat(b_mat_in), + c_mat(c_mat_in), + e_mat(e_mat_in), + pre_amg(pre_amg_in), + a_inv_direct(a_inv_direct_in), + a_inv_mf_gmg(a_inv_mf_gmg_in) +{ + +} + +template +void +HMatrix::initialize(LA::MPI::Vector &exemplar_density_vector, LA::MPI::Vector &exemplar_displacement_vector, LA::MPI::SparseMatrix &d_m_inv_mat_in) +{ + temp_vect_1 = exemplar_displacement_vector; + temp_vect_2 = exemplar_displacement_vector; + temp_vect_3 = exemplar_displacement_vector; + temp_vect_4 = exemplar_displacement_vector; + temp_vect_5 = exemplar_density_vector; + temp_vect_6 = exemplar_density_vector; + temp_vect_7 = exemplar_density_vector; + + temp_vect_1 = 0.; + temp_vect_2 = 0.; + temp_vect_3 = 0.; + temp_vect_4 = 0.; + temp_vect_5 = 0.; + temp_vect_6 = 0.; + temp_vect_7 = 0.; + + d_m_inv_mat.copy_from(d_m_inv_mat_in); + + +} + +template +void HMatrix::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + if (Input::solver_choice==SolverOptions::inexact_K_with_inexact_A_gmres) + { + c_mat.vmult(temp_vect_1,src); + e_mat.vmult(temp_vect_2,src); + + try + { + a_inv_mf_gmg.vmult(temp_vect_3,temp_vect_1); + } catch (std::exception &exc) + { + + } + try + { + a_inv_mf_gmg.vmult(temp_vect_4,temp_vect_2); + } catch (std::exception &exc) + { + + } + + c_mat.Tvmult(temp_vect_6,temp_vect_4); + e_mat.Tvmult(temp_vect_5,temp_vect_3); + + b_mat.vmult(temp_vect_7,src); + dst = temp_vect_7 + temp_vect_6 + temp_vect_5; + } + else + { + c_mat.vmult(temp_vect_1,src); + e_mat.vmult(temp_vect_2,src); + + a_inv_direct.vmult(temp_vect_3,temp_vect_1); + a_inv_direct.vmult(temp_vect_4,temp_vect_2); + + c_mat.Tvmult(temp_vect_6,temp_vect_4); + e_mat.Tvmult(temp_vect_5,temp_vect_3); + + b_mat.vmult(temp_vect_7,src); + dst = temp_vect_7 - temp_vect_6 - temp_vect_5; + } + +} + + +template +void HMatrix::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + if (Input::solver_choice==SolverOptions::inexact_K_with_inexact_A_gmres) + { + c_mat.vmult(temp_vect_1,src); + e_mat.vmult(temp_vect_2,src); + + + try + { + a_inv_mf_gmg.vmult(temp_vect_3,temp_vect_1); + + } catch (std::exception &exc) + { + + } + try + { + a_inv_mf_gmg.vmult(temp_vect_4,temp_vect_2); + + } catch (std::exception &exc) + { + + } + + c_mat.Tvmult(temp_vect_6,temp_vect_4); + e_mat.Tvmult(temp_vect_5,temp_vect_3); + + b_mat.vmult(temp_vect_7,src); + dst = temp_vect_7 + temp_vect_6 + temp_vect_5; + } + else + { + c_mat.vmult(temp_vect_1,src); + e_mat.vmult(temp_vect_2,src); + + a_inv_direct.vmult(temp_vect_3,temp_vect_1); + a_inv_direct.vmult(temp_vect_4,temp_vect_2); + + c_mat.Tvmult(temp_vect_6,temp_vect_4); + e_mat.Tvmult(temp_vect_5,temp_vect_3); + + b_mat.vmult(temp_vect_7,src); + dst = temp_vect_7 - temp_vect_6 - temp_vect_5; + } + +} + + +// ****************** HMatrixDirect *********************** + +template +HMatrixDirect::HMatrixDirect(LA::MPI::SparseMatrix &a_mat_in, const LA::MPI::SparseMatrix &b_mat_in, const LA::MPI::SparseMatrix &c_mat_in, const LA::MPI::SparseMatrix &e_mat_in,TrilinosWrappers::PreconditionAMG &pre_amg_in,VmultTrilinosSolverDirect &a_inv_direct_in, AInvMatMFGMG &a_inv_mf_gmg_in) + : + a_mat(a_mat_in), + b_mat(b_mat_in), + c_mat(c_mat_in), + e_mat(e_mat_in), + pre_amg(pre_amg_in), + a_inv_direct(a_inv_direct_in), + a_inv_mf_gmg(a_inv_mf_gmg_in) +{ + +} + +template +void +HMatrixDirect::initialize(LA::MPI::Vector &exemplar_density_vector, LA::MPI::Vector &exemplar_displacement_vector) +{ + temp_vect_1 = exemplar_displacement_vector; + temp_vect_2 = exemplar_displacement_vector; + temp_vect_3 = exemplar_displacement_vector; + temp_vect_4 = exemplar_displacement_vector; + temp_vect_5 = exemplar_density_vector; + temp_vect_6 = exemplar_density_vector; + temp_vect_7 = exemplar_density_vector; + + +} + +template +void HMatrixDirect::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + + c_mat.vmult(temp_vect_1,src); + e_mat.vmult(temp_vect_2,src); + + a_inv_direct.vmult(temp_vect_3,temp_vect_1); + a_inv_direct.vmult(temp_vect_4,temp_vect_2); + + c_mat.Tvmult(temp_vect_6,temp_vect_4); + e_mat.Tvmult(temp_vect_5,temp_vect_3); + + b_mat.vmult(temp_vect_7,src); + dst = temp_vect_7 - temp_vect_6 - temp_vect_5; + + + +} + + +template +void HMatrixDirect::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + + c_mat.vmult(temp_vect_1,src); + e_mat.vmult(temp_vect_2,src); + + a_inv_direct.vmult(temp_vect_3,temp_vect_1); + a_inv_direct.vmult(temp_vect_4,temp_vect_2); + + c_mat.Tvmult(temp_vect_6,temp_vect_4); + e_mat.Tvmult(temp_vect_5,temp_vect_3); + + b_mat.vmult(temp_vect_7,src); + dst = temp_vect_7 - temp_vect_6 - temp_vect_5; + +} + + +// ****************** JinvMatrix *********************** +template +JinvMatrix::JinvMatrix(HMatrix &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in) + : + h_mat(h_mat_in), + g_mat(g_mat_in), + d_m_mat(d_m_mat_in) +{ + +} + +template +void +JinvMatrix::initialize(LA::MPI::Vector &exemplar_density_vector, LA::MPI::SparseMatrix &d_m_inv_mat_in) +{ + temp_vect_1 = exemplar_density_vector; + temp_vect_2 = exemplar_density_vector; + temp_vect_3 = exemplar_density_vector; + temp_vect_4 = exemplar_density_vector; + d_m_inv_mat.copy_from(d_m_inv_mat_in); + +} + +template +void JinvMatrix::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + g_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + h_mat.vmult(temp_vect_3,temp_vect_2); + d_m_mat.vmult(temp_vect_4,src); + + dst = -1*temp_vect_4 - temp_vect_3; +} + +template +void JinvMatrix::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + h_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + g_mat.vmult(temp_vect_3,temp_vect_2); + d_m_mat.vmult(temp_vect_4,src); + + dst = -1*temp_vect_4 - temp_vect_3; +} + +// ****************** JinvMatrixPart *********************** +template +JinvMatrixPart::JinvMatrixPart(HMatrix &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in) + : + h_mat(h_mat_in), + g_mat(g_mat_in), + d_m_mat(d_m_mat_in) +{ +} + +template +void +JinvMatrixPart::initialize(LA::MPI::Vector &exemplar_density_vector) +{ + temp_vect_1 = exemplar_density_vector; + temp_vect_2 = exemplar_density_vector; + temp_vect_3 = exemplar_density_vector; + temp_vect_4 = exemplar_density_vector; + +} + +template +void JinvMatrixPart::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + g_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + h_mat.vmult(temp_vect_3,temp_vect_2); + d_m_inv_mat.vmult(temp_vect_4,temp_vect_3); + dst = -1 * temp_vect_4; +} + +template +void JinvMatrixPart::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + h_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + g_mat.vmult(temp_vect_3,temp_vect_2); + d_m_inv_mat.vmult(temp_vect_4,temp_vect_3); + dst = -1 * temp_vect_4; +} + +// ****************** JinvMatrixDirect *********************** +template +JinvMatrixDirect::JinvMatrixDirect(HMatrixDirect &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in) + : + h_mat(h_mat_in), + g_mat(g_mat_in), + d_m_mat(d_m_mat_in) +{ + +} + +template +void +JinvMatrixDirect::initialize(LA::MPI::Vector &exemplar_density_vector) +{ + temp_vect_1 = exemplar_density_vector; + temp_vect_2 = exemplar_density_vector; + temp_vect_3 = exemplar_density_vector; + temp_vect_4 = exemplar_density_vector; + +} + +template +void JinvMatrixDirect::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + g_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + h_mat.vmult(temp_vect_3,temp_vect_2); + d_m_mat.vmult(temp_vect_4,src); + + dst = -1*temp_vect_4 - temp_vect_3; +} + +template +void JinvMatrixDirect::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + h_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + g_mat.vmult(temp_vect_3,temp_vect_2); + d_m_mat.vmult(temp_vect_4,src); + + dst = -1*temp_vect_4 - temp_vect_3; +} + + +// ****************** KinvMatrix *********************** +template +KinvMatrix::KinvMatrix(HMatrix &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in) + : + h_mat(h_mat_in), + g_mat(g_mat_in), + d_m_mat(d_m_mat_in) +{ + +} + +template +void +KinvMatrix::initialize(LA::MPI::Vector &exemplar_density_vector, LA::MPI::SparseMatrix &d_m_inv_mat_in) +{ + + temp_vect_1 = exemplar_density_vector; + temp_vect_2 = exemplar_density_vector; + temp_vect_3 = exemplar_density_vector; + temp_vect_4 = exemplar_density_vector; + + temp_vect_1 = 0.; + temp_vect_2 = 0.; + temp_vect_3 = 0.; + temp_vect_4 = 0.; + + d_m_inv_mat.copy_from(d_m_inv_mat_in); + +} + +template +void KinvMatrix::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + h_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + g_mat.vmult(temp_vect_3,temp_vect_2); + d_m_mat.vmult(temp_vect_4,src); + + dst = -1*temp_vect_4 - temp_vect_3; +} + +template +void KinvMatrix::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + g_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + h_mat.vmult(temp_vect_3,temp_vect_2); + d_m_mat.vmult(temp_vect_4,src); + + dst = -1*temp_vect_4 - temp_vect_3; +} + +// ****************** KinvMatrixPart *********************** +template +KinvMatrixPart::KinvMatrixPart(HMatrix &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in) + : + h_mat(h_mat_in), + g_mat(g_mat_in), + d_m_mat(d_m_mat_in) +{ + +} + +template +void +KinvMatrixPart::initialize(LA::MPI::Vector &exemplar_density_vector) +{ + temp_vect_1 = exemplar_density_vector; + temp_vect_2 = exemplar_density_vector; + temp_vect_3 = exemplar_density_vector; + temp_vect_4 = exemplar_density_vector; + +} + +template +void KinvMatrixPart::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + + h_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + g_mat.vmult(temp_vect_3,temp_vect_2); + d_m_inv_mat.vmult(temp_vect_4,temp_vect_3); + + dst = temp_vect_4; + //kinv is -Dm-GDmH + //kinv scaled would be I+DminvGDminvH + //kpart is -DminvGDmH + +} + +template +void KinvMatrixPart::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + g_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + h_mat.vmult(temp_vect_3,temp_vect_2); + d_m_inv_mat.vmult(temp_vect_4,temp_vect_3); + + dst = temp_vect_4; + //kinv is -Dm-GDmH + //kinv scaled would be I+DminvGDminvH + //kpart is -DminvGDmH +} + +// ****************** KinvMatrixDirect *********************** +template +KinvMatrixDirect::KinvMatrixDirect(HMatrixDirect &h_mat_in, GMatrix &g_mat_in, const LA::MPI::SparseMatrix &d_m_mat_in) + : + h_mat(h_mat_in), + g_mat(g_mat_in), + d_m_mat(d_m_mat_in) +{ + +} + +template +void +KinvMatrixDirect::initialize(LA::MPI::Vector &exemplar_density_vector) +{ + temp_vect_1 = exemplar_density_vector; + temp_vect_2 = exemplar_density_vector; + temp_vect_3 = exemplar_density_vector; + temp_vect_4 = exemplar_density_vector; + +} + +template +void KinvMatrixDirect::vmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + h_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + g_mat.vmult(temp_vect_3,temp_vect_2); + d_m_mat.vmult(temp_vect_4,src); + + dst = -1*temp_vect_4 - temp_vect_3; +} + +template +void KinvMatrixDirect::Tvmult(LA::MPI::Vector &dst, const LA::MPI::Vector &src) const +{ + g_mat.vmult(temp_vect_1,src); + d_m_inv_mat.vmult(temp_vect_2,temp_vect_1); + h_mat.vmult(temp_vect_3,temp_vect_2); + d_m_mat.vmult(temp_vect_4,src); + + dst = -1*temp_vect_4 - temp_vect_3; +} + + + + +//************************************************** + + +AMatWrapped::AMatWrapped(LA::MPI::SparseMatrix &a_mat_in) + : + a_mat(a_mat_in) +{ +} +void AMatWrapped::vmult(LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const +{ + ChangeVectorTypes::copy(temp_src, src); + a_mat.vmult(temp_dst,temp_src); + ChangeVectorTypes::copy(dst, temp_dst); +} +void AMatWrapped::Tvmult(LinearAlgebra::distributed::Vector &dst, const LinearAlgebra::distributed::Vector &src) const +{ + ChangeVectorTypes::copy(temp_src, src); + a_mat.Tvmult(temp_dst,temp_src); + ChangeVectorTypes::copy(dst, temp_dst); +} +} + + + +template class SAND::TopOptSchurPreconditioner<2>; +template class SAND::TopOptSchurPreconditioner<3>; + +template class SAND::AInvMatMFGMG<2>; +template class SAND::AInvMatMFGMG<3>; + +template class SAND::JinvMatrix<2>; +template class SAND::JinvMatrix<3>; + +template class SAND::KinvMatrix<2>; +template class SAND::KinvMatrix<3>; + +template class SAND::HMatrix<2>; +template class SAND::HMatrix<3>; + diff --git a/SAND_top_opt/source/watchdog.cc b/SAND_top_opt/source/watchdog.cc new file mode 100644 index 00000000..8ba8657f --- /dev/null +++ b/SAND_top_opt/source/watchdog.cc @@ -0,0 +1,478 @@ +#include +#include +#include +#include +#include +#include "../include/markov_filter.h" +#include "../include/kkt_system.h" +#include "../include/input_information.h" +#include "../include/watchdog.h" +#include +#include +#include +#include +#include + +///Above are fairly normal files to include. I also use the sparse direct package, which requiresBLAS/LAPACK +/// to perform a direct solve while I work on a fast iterative solver for this problem. + +namespace SAND { + namespace LA + { + using namespace dealii::LinearAlgebraTrilinos; + } + + + using namespace dealii; + + ///Constructor + template + NonlinearWatchdog::NonlinearWatchdog() + : + mpi_communicator(MPI_COMM_WORLD), + pcout(std::cout,(Utilities::MPI::this_mpi_process(mpi_communicator) == 0)), + overall_timer(pcout, TimerOutput::never, TimerOutput::wall_times) + { + } + + ///A binary search figures out the maximum step that meets the dual feasibility - that s>0 and z>0. The fraction to boundary increases as the barrier size decreases. + + template + std::pair + NonlinearWatchdog::calculate_max_step_size(const LA::MPI::BlockVector &state, const LA::MPI::BlockVector &step) const { + + double step_size_s_low = 0; + double step_size_z_low = 0; + double step_size_s_high = 1; + double step_size_z_high = 1; + double step_size_s, step_size_z; + LA::MPI::BlockVector state_test_s = state; + state_test_s = 0; + LA::MPI::BlockVector state_test_z = state; + state_test_z = 0; + for (unsigned int k = 0; k < 50; k++) + { + + step_size_s = (step_size_s_low + step_size_s_high) / 2; + step_size_z = (step_size_z_low + step_size_z_high) / 2; + const LA::MPI::BlockVector state_test_s = + (Input::fraction_to_boundary * state) + (step_size_s * step); + + const LA::MPI::BlockVector state_test_z = + (Input::fraction_to_boundary * state) + (step_size_z * step); + + const bool accept_s = (state_test_s.block(SolutionBlocks::density_lower_slack).is_non_negative()) + && (state_test_s.block(SolutionBlocks::density_upper_slack).is_non_negative()); + const bool accept_z = (state_test_z.block(SolutionBlocks::density_lower_slack_multiplier).is_non_negative()) + && (state_test_z.block(SolutionBlocks::density_upper_slack_multiplier).is_non_negative()); + + if (accept_s) { + step_size_s_low = step_size_s; + } else { + step_size_s_high = step_size_s; + } + if (accept_z) { + step_size_z_low = step_size_z; + } else { + step_size_z_high = step_size_z; + } + } + pcout << "s step : " << step_size_s_low << " z size : " << step_size_z_low << std::endl; + return {step_size_s_low, step_size_z_low}; + } + + ///Creates a rhs vector that we can use to look at the magnitude of the KKT conditions. This is then used for testing the convergence before shrinking barrier size, as well as in the calculation of the l1 merit. + + template + const LA::MPI::BlockVector + NonlinearWatchdog::find_max_step(const LA::MPI::BlockVector &state) + { + TimerOutput::Scope t(overall_timer, "find step"); + { + TimerOutput::Scope t(overall_timer, "assemble"); + kkt_system.assemble_block_system(state, barrier_size); + } + + pcout << "pre" << std::endl; + LA::MPI::BlockVector step; + { + TimerOutput::Scope t(overall_timer, "solve"); + step = kkt_system.solve(state); + } + + pcout << "post" << std::endl; + const auto max_step_sizes= calculate_max_step_size(state,step); + const double step_size_s = max_step_sizes.first; + const double step_size_z = max_step_sizes.second; + LA::MPI::BlockVector max_step(10); + + max_step.block(SolutionBlocks::density) = step_size_s * step.block(SolutionBlocks::density); + max_step.block(SolutionBlocks::displacement) = step_size_s * step.block(SolutionBlocks::displacement); + max_step.block(SolutionBlocks::unfiltered_density) = step_size_s * step.block(SolutionBlocks::unfiltered_density); + max_step.block(SolutionBlocks::density_lower_slack) = step_size_s * step.block(SolutionBlocks::density_lower_slack); + max_step.block(SolutionBlocks::density_upper_slack) = step_size_s * step.block(SolutionBlocks::density_upper_slack); + max_step.block(SolutionBlocks::unfiltered_density_multiplier) = step_size_z * step.block(SolutionBlocks::unfiltered_density_multiplier); + max_step.block(SolutionBlocks::density_lower_slack_multiplier) = step_size_z * step.block(SolutionBlocks::density_lower_slack_multiplier); + max_step.block(SolutionBlocks::density_upper_slack_multiplier) = step_size_z * step.block(SolutionBlocks::density_upper_slack_multiplier); + max_step.block(SolutionBlocks::displacement_multiplier) = step_size_z * step.block(SolutionBlocks::displacement_multiplier); + max_step.block(SolutionBlocks::total_volume_multiplier) = step_size_z * step.block(SolutionBlocks::total_volume_multiplier); + + pcout << "here" << std::endl; + return max_step; + } + + ///This is a simple back-stepping algorithm for a line search - keeps shrinking step size until it finds a step where the markov filter requirement is met. + + template + LA::MPI::BlockVector + NonlinearWatchdog::take_scaled_step(const LA::MPI::BlockVector &state,const LA::MPI::BlockVector &max_step) const + { + double step_size = 1; + for(unsigned int k = 0; k<10; k++) + { + if(markov_filter.check_filter(kkt_system.calculate_objective_value(state), kkt_system.calculate_barrier_distance(state), kkt_system.calculate_feasibility(state,barrier_size))) + { + break; + } + else + { + step_size = step_size/2; + } + } + return state + (step_size * max_step); + + } + + + + ///Checks to see if the KKT conditions are sufficiently met to lower barrier size. + template + bool + NonlinearWatchdog::check_convergence(const LA::MPI::BlockVector &state) const + { + if (kkt_system.calculate_convergence(state) < Input::required_norm) + { + return true; + } + else + { + return false; + } + } + + ///This updates the barrier value using the selected barrier scheme - more work could be done to optimize + /// the performance of the mixed method + template + void + NonlinearWatchdog::update_barrier(LA::MPI::BlockVector ¤t_state) + { + ///The LOQO scheme uses information about the similarity of the slack/slack multiplier product as a + /// heuristic for decreasing barrier value + if (Input::barrier_reduction == BarrierOptions::loqo) + { + double loqo_min = 1000; + double loqo_average; + double lower_prod; + double full_lower_prod; + double upper_prod; + double full_upper_prod; + unsigned int vect_size = current_state.block(SolutionBlocks::density_lower_slack).size(); + for(unsigned int k = 0; k < vect_size; k++) + { + lower_prod = 1; + if (current_state.block(SolutionBlocks::density_lower_slack).in_local_range(k)) + lower_prod=lower_prod * current_state.block(SolutionBlocks::density_lower_slack)[k]; + if (current_state.block(SolutionBlocks::density_lower_slack_multiplier).in_local_range(k)) + lower_prod=lower_prod * current_state.block(SolutionBlocks::density_lower_slack_multiplier)[k]; + + upper_prod=1; + if (current_state.block(SolutionBlocks::density_upper_slack).in_local_range(k)) + upper_prod=upper_prod * current_state.block(SolutionBlocks::density_upper_slack)[k]; + if (current_state.block(SolutionBlocks::density_upper_slack_multiplier).in_local_range(k)) + upper_prod=upper_prod * current_state.block(SolutionBlocks::density_upper_slack_multiplier)[k]; + + MPI_Allreduce(&lower_prod, &full_lower_prod, 1, MPI_DOUBLE, MPI_PROD, MPI_COMM_WORLD); + MPI_Allreduce(&upper_prod, &full_upper_prod, 1, MPI_DOUBLE, MPI_PROD, MPI_COMM_WORLD); + if (full_lower_prod < loqo_min) + { + loqo_min = full_lower_prod; + } + if (full_upper_prod < loqo_min) + { + loqo_min = full_upper_prod; + } + } + loqo_average = (current_state.block(SolutionBlocks::density_lower_slack)*current_state.block(SolutionBlocks::density_lower_slack_multiplier) + + current_state.block(SolutionBlocks::density_upper_slack)*current_state.block(SolutionBlocks::density_upper_slack_multiplier) + )/(2*vect_size); + double loqo_complimentarity_deviation = loqo_min/loqo_average; + pcout << "loqo cd: " << loqo_complimentarity_deviation << std::endl; + double loqo_multiplier; + if((.05 * (1-loqo_complimentarity_deviation)/loqo_complimentarity_deviation)<2) + { + loqo_multiplier = .1*std::pow((.05 * (1-loqo_complimentarity_deviation)/loqo_complimentarity_deviation),3); + } + else + { + loqo_multiplier = .8; + } + pcout << "loqo mult: " << loqo_multiplier << std::endl; + if (loqo_multiplier< 0) + { + barrier_size = std::abs(loqo_multiplier) * loqo_average; + } + else + { + barrier_size = loqo_multiplier * loqo_average; + } + if (barrier_size < Input::min_barrier_size) + { + barrier_size=Input::min_barrier_size; + } + } + + ///The monotome scheme fully solves the problem with one barrier size before decreasing the + /// barrier and starting again + if (Input::barrier_reduction == BarrierOptions::monotone) + { + if (kkt_system.calculate_rhs_norm(current_state,barrier_size) < barrier_size * 1e-3) + { + barrier_size = barrier_size * .7; + } + if (barrier_size < Input::min_barrier_size) + { + barrier_size=Input::min_barrier_size; + } + } + + ///The mixed method uses LOQO unless it gets stuck, at which point it switches to monotone, allowing for an adaptive method + /// that still globally converges the barrier value to 0. + if (Input::barrier_reduction == BarrierOptions::mixed) + { + if (mixed_barrier_monotone_mode) + { + if (kkt_system.calculate_rhs_norm(current_state,barrier_size) < barrier_size) + { + barrier_size = barrier_size * .8; + mixed_barrier_monotone_mode=false; + pcout << "monotone mode turned off" << std::endl; + } + } + else + { + double loqo_min = 1000; + double loqo_average; + unsigned int vect_size = current_state.block(SolutionBlocks::density_lower_slack).size(); + double lower_prod, full_lower_prod, upper_prod, full_upper_prod; + for(unsigned int k = 0; k < vect_size; k++) + { + lower_prod = 1; + if (current_state.block(SolutionBlocks::density_lower_slack).in_local_range(k)) + lower_prod=lower_prod * current_state.block(SolutionBlocks::density_lower_slack)[k]; + if (current_state.block(SolutionBlocks::density_lower_slack_multiplier).in_local_range(k)) + lower_prod=lower_prod * current_state.block(SolutionBlocks::density_lower_slack_multiplier)[k]; + + upper_prod=1; + if (current_state.block(SolutionBlocks::density_upper_slack).in_local_range(k)) + upper_prod=upper_prod * current_state.block(SolutionBlocks::density_upper_slack)[k]; + if (current_state.block(SolutionBlocks::density_upper_slack_multiplier).in_local_range(k)) + upper_prod=upper_prod * current_state.block(SolutionBlocks::density_upper_slack_multiplier)[k]; + + MPI_Allreduce(&lower_prod, &full_lower_prod, 1, MPI_DOUBLE, MPI_PROD, MPI_COMM_WORLD); + MPI_Allreduce(&upper_prod, &full_upper_prod, 1, MPI_DOUBLE, MPI_PROD, MPI_COMM_WORLD); + + if (full_lower_prod < loqo_min) + { + loqo_min = full_lower_prod; + } + if (full_upper_prod < loqo_min) + { + loqo_min = full_upper_prod; + } + } + loqo_average = (current_state.block(SolutionBlocks::density_lower_slack)*current_state.block(SolutionBlocks::density_lower_slack_multiplier) + + current_state.block(SolutionBlocks::density_upper_slack)*current_state.block(SolutionBlocks::density_upper_slack_multiplier) + )/(2*vect_size); + double loqo_complimentarity_deviation = loqo_min/loqo_average; + double loqo_multiplier; + if((.05 * (1-loqo_complimentarity_deviation)/loqo_complimentarity_deviation)<2) + { + loqo_multiplier = .1*std::pow((.05 * (1-loqo_complimentarity_deviation)/loqo_complimentarity_deviation),3); + } + else + { + loqo_multiplier = 1/.8; + mixed_barrier_monotone_mode = true; + pcout << "monotone mode turned on" << std::endl; + } + if (loqo_multiplier<.01) + { + barrier_size = .01 * loqo_average; + } + else + { + barrier_size = loqo_multiplier * loqo_average; + } + if (barrier_size < Input::min_barrier_size) + { + barrier_size=Input::min_barrier_size; + } + } + } + + } + + template + void + NonlinearWatchdog::perform_initial_setup() + { + barrier_size = Input::initial_barrier_size; + kkt_system.create_triangulation(); + kkt_system.setup_boundary_values(); + pcout << "setup kkt system" << std::endl; + kkt_system.setup_block_system(); + pcout << "setup kkt system" << std::endl; + + if (Input::barrier_reduction==BarrierOptions::mixed) + { + mixed_barrier_monotone_mode = false; + } + } + + + template + void + NonlinearWatchdog::nonlinear_step(LA::MPI::BlockVector ¤t_state, LA::MPI::BlockVector ¤t_step, const unsigned int max_uphill_steps, unsigned int &iteration_number) + { + + bool converged = false; + //while not converged + while(!converged && iteration_number < Input::max_steps) + { + bool found_step = false; + //save current state as watchdog state + + const LA::MPI::BlockVector watchdog_state = current_state; + LA::MPI::BlockVector watchdog_step; + //for 1-8 steps - this is the number of steps away we will let it go uphill before demanding downhill + for(unsigned int k = 0; k + void + NonlinearWatchdog::run() { + overall_timer.enter_subsection("Total Time"); + + perform_initial_setup(); + + const unsigned int max_uphill_steps = 8; + unsigned int iteration_number = 0; + + //while barrier value above minimal value and total iterations under some value + LA::MPI::BlockVector current_state = kkt_system.get_initial_state(); + LA::MPI::BlockVector current_step; + + markov_filter.setup(kkt_system.calculate_objective_value(current_state), kkt_system.calculate_barrier_distance(current_state), kkt_system.calculate_feasibility(current_state,barrier_size), barrier_size); + + // std::cout << "finished setup - beginning watchdog steps" << std::endl; + + while((barrier_size > Input::min_barrier_size || !check_convergence(current_state)) && iteration_number < Input::max_steps) + { + nonlinear_step(current_state, current_step, max_uphill_steps, iteration_number); + } +// kkt_system.output_stl(current_state); + } + +} // namespace SAND + + +template class SAND::NonlinearWatchdog<2>; +template class SAND::NonlinearWatchdog<3>;