diff --git a/dune/gdt/local/assembler/bilinear-form-assemblers.hh b/dune/gdt/local/assembler/bilinear-form-assemblers.hh index 5a0db390f9e3789cda72070f98a45840d90a8f62..816cf0fc6165cd0556c59df484b8c36e58766692 100644 --- a/dune/gdt/local/assembler/bilinear-form-assemblers.hh +++ b/dune/gdt/local/assembler/bilinear-form-assemblers.hh @@ -57,6 +57,7 @@ public: , local_bilinear_form_(local_two_form.copy()) , global_matrix_(global_matrix) , param_(param) + , scaling_(param_.has_key("matrixoperator.scaling") ? param_.get("matrixoperator.scaling").at(0) : 1.) , local_matrix_(test_space_.mapper().max_local_size(), ansatz_space_.mapper().max_local_size()) , global_test_indices_(test_space_.mapper().max_local_size()) , global_ansatz_indices_(ansatz_space_.mapper().max_local_size()) @@ -82,6 +83,7 @@ public: , local_bilinear_form_(other.local_bilinear_form_->copy()) , global_matrix_(other.global_matrix_) , param_(other.param_) + , scaling_(other.scaling_) , local_matrix_(test_space_.mapper().max_local_size(), ansatz_space_.mapper().max_local_size()) , global_test_indices_(test_space_.mapper().max_local_size()) , global_ansatz_indices_(ansatz_space_.mapper().max_local_size()) @@ -108,7 +110,8 @@ public: ansatz_space_.mapper().global_indices(element, global_ansatz_indices_); for (size_t ii = 0; ii < test_basis_->size(param_); ++ii) for (size_t jj = 0; jj < ansatz_basis_->size(param_); ++jj) - global_matrix_.add_to_entry(global_test_indices_[ii], global_ansatz_indices_[jj], local_matrix_[ii][jj]); + global_matrix_.add_to_entry( + global_test_indices_[ii], global_ansatz_indices_[jj], scaling_ * local_matrix_[ii][jj]); } // ... apply_local(...) private: @@ -117,6 +120,7 @@ private: const std::unique_ptr<LocalBilinearFormType> local_bilinear_form_; MatrixType& global_matrix_; XT::Common::Parameter param_; + const double scaling_; DynamicMatrix<FieldType> local_matrix_; DynamicVector<size_t> global_test_indices_; DynamicVector<size_t> global_ansatz_indices_; @@ -164,6 +168,7 @@ public: , local_bilinear_form_(local_two_form.copy()) , global_matrix_(global_matrix) , param_(param) + , scaling_(param_.has_key("matrixoperator.scaling") ? param_.get("matrixoperator.scaling").at(0) : 1.) , local_matrix_in_in_(test_space_.mapper().max_local_size(), ansatz_space_.mapper().max_local_size()) , local_matrix_in_out_(test_space_.mapper().max_local_size(), ansatz_space_.mapper().max_local_size()) , local_matrix_out_in_(test_space_.mapper().max_local_size(), ansatz_space_.mapper().max_local_size()) @@ -196,6 +201,7 @@ public: , local_bilinear_form_(other.local_bilinear_form_->copy()) , global_matrix_(other.global_matrix_) , param_(other.param_) + , scaling_(other.scaling_) , local_matrix_in_in_(test_space_.mapper().max_local_size(), ansatz_space_.mapper().max_local_size()) , local_matrix_in_out_(test_space_.mapper().max_local_size(), ansatz_space_.mapper().max_local_size()) , local_matrix_out_in_(test_space_.mapper().max_local_size(), ansatz_space_.mapper().max_local_size()) @@ -245,18 +251,18 @@ public: for (size_t ii = 0; ii < test_basis_inside_->size(param_); ++ii) { for (size_t jj = 0; jj < ansatz_basis_inside_->size(param_); ++jj) global_matrix_.add_to_entry( - global_test_indices_in_[ii], global_ansatz_indices_in_[jj], local_matrix_in_in_[ii][jj]); + global_test_indices_in_[ii], global_ansatz_indices_in_[jj], scaling_ * local_matrix_in_in_[ii][jj]); for (size_t jj = 0; jj < ansatz_basis_outside_->size(param_); ++jj) global_matrix_.add_to_entry( - global_test_indices_in_[ii], global_ansatz_indices_out_[jj], local_matrix_in_out_[ii][jj]); + global_test_indices_in_[ii], global_ansatz_indices_out_[jj], scaling_ * local_matrix_in_out_[ii][jj]); } for (size_t ii = 0; ii < test_basis_outside_->size(param_); ++ii) { for (size_t jj = 0; jj < ansatz_basis_inside_->size(param_); ++jj) global_matrix_.add_to_entry( - global_test_indices_out_[ii], global_ansatz_indices_in_[jj], local_matrix_out_in_[ii][jj]); + global_test_indices_out_[ii], global_ansatz_indices_in_[jj], scaling_ * local_matrix_out_in_[ii][jj]); for (size_t jj = 0; jj < ansatz_basis_outside_->size(param_); ++jj) global_matrix_.add_to_entry( - global_test_indices_out_[ii], global_ansatz_indices_out_[jj], local_matrix_out_out_[ii][jj]); + global_test_indices_out_[ii], global_ansatz_indices_out_[jj], scaling_ * local_matrix_out_out_[ii][jj]); } } // ... apply_local(...) @@ -266,6 +272,7 @@ private: const std::unique_ptr<LocalBilinearFormType> local_bilinear_form_; MatrixType& global_matrix_; XT::Common::Parameter param_; + const double scaling_; DynamicMatrix<FieldType> local_matrix_in_in_; DynamicMatrix<FieldType> local_matrix_in_out_; DynamicMatrix<FieldType> local_matrix_out_in_; diff --git a/dune/gdt/local/assembler/operator-fd-jacobian-assemblers.hh b/dune/gdt/local/assembler/operator-fd-jacobian-assemblers.hh index 5920a575c558d4f307ff61b72a966cafcd9567d4..509adb00eb959a1117d77a11f431f41d8fe1d544 100644 --- a/dune/gdt/local/assembler/operator-fd-jacobian-assemblers.hh +++ b/dune/gdt/local/assembler/operator-fd-jacobian-assemblers.hh @@ -82,6 +82,7 @@ public: , source_vector_(source_vector) , local_op_(local_operator.copy()) , param_(param) + , scaling_(param_.has_key("matrixoperator.scaling") ? param_.get("matrixoperator.scaling").at(0) : 1.) , eps_(param_.has_key("finite-difference-jacobians.eps") ? param_.get("finite-difference-jacobians.eps").at(0) : 1e-7) , source_(source_space_) // This is a full vector, and intended! @@ -99,6 +100,7 @@ public: , source_vector_(other.source_vector_) , local_op_(other.local_op_->copy()) , param_(other.param_) + , scaling_(other.scaling_) , eps_(other.eps_) , source_(source_space_) // This is a full vector, and intended! , range_(range_space_) // This is a full vector, and intended! @@ -143,7 +145,7 @@ public: auto derivative = (local_range_->dofs()[ii] - range_DoFs_[ii]) / eps; if (XT::Common::FloatCmp::eq(derivative, eps)) derivative = 0; - matrix_.add_to_entry(global_range_indices_[ii], global_source_indices_[jj], derivative); + matrix_.add_to_entry(global_range_indices_[ii], global_source_indices_[jj], scaling_ * derivative); } // restore source local_source_->dofs()[jj] = jjth_source_DoF; @@ -157,6 +159,7 @@ private: const VectorType& source_vector_; const std::unique_ptr<LocalElementOperatorType> local_op_; const XT::Common::Parameter param_; + const double scaling_; const real_t<F> eps_; DiscreteFunction<V, SGV, s_r, s_rC, F> source_; DiscreteFunction<V, RGV, r_r, r_rC, F> range_; @@ -226,6 +229,7 @@ public: , source_vector_(source_vector) , local_op_(local_operator.copy()) , param_(param) + , scaling_(param_.has_key("matrixoperator.scaling") ? param_.get("matrixoperator.scaling").at(0) : 1.) , eps_(eps) , source_(source_space_) // This is a full vector, and intended! , range_(range_space_) // This is a full vector, and intended! @@ -244,6 +248,7 @@ public: , source_vector_(other.source_vector_) , local_op_(other.local_op_->copy()) , param_(other.param_) + , scaling_(other.scaling_) , eps_(other.eps_) , source_(source_space_) // This is a full vector, and intended! , range_(range_space_) // This is a full vector, and intended! @@ -309,7 +314,8 @@ public: auto derivative = (local_range_inside_->dofs()[ii] - range_DoFs_inside_[ii]) / eps; if (XT::Common::FloatCmp::eq(derivative, eps)) derivative = 0; - matrix_.add_to_entry(global_range_indices_inside_[ii], global_source_indices_inside_[jj], derivative); + matrix_.add_to_entry( + global_range_indices_inside_[ii], global_source_indices_inside_[jj], scaling_ * derivative); } // observe perturbation in outside range DoFs if (treat_outside) { @@ -317,7 +323,8 @@ public: auto derivative = (local_range_outside_->dofs()[ii] - range_DoFs_outside_[ii]) / eps; if (XT::Common::FloatCmp::eq(derivative, eps)) derivative = 0; - matrix_.add_to_entry(global_range_indices_outside_[ii], global_source_indices_inside_[jj], derivative); + matrix_.add_to_entry( + global_range_indices_outside_[ii], global_source_indices_inside_[jj], scaling_ * derivative); } } // restore source @@ -340,14 +347,16 @@ public: auto derivative = (local_range_inside_->dofs()[ii] - range_DoFs_inside_[ii]) / eps; if (XT::Common::FloatCmp::eq(derivative, eps)) derivative = 0; - matrix_.add_to_entry(global_range_indices_inside_[ii], global_source_indices_outside_[jj], derivative); + matrix_.add_to_entry( + global_range_indices_inside_[ii], global_source_indices_outside_[jj], scaling_ * derivative); } // observe perturbation in outside range DoFs for (size_t ii = 0; ii < local_range_outside_size; ++ii) { auto derivative = (local_range_outside_->dofs()[ii] - range_DoFs_outside_[ii]) / eps; if (XT::Common::FloatCmp::eq(derivative, eps)) derivative = 0; - matrix_.add_to_entry(global_range_indices_outside_[ii], global_source_indices_outside_[jj], derivative); + matrix_.add_to_entry( + global_range_indices_outside_[ii], global_source_indices_outside_[jj], scaling_ * derivative); } // restore source local_source_outside_->dofs()[jj] = jjth_source_DoF; @@ -362,6 +371,7 @@ private: const VectorType& source_vector_; const std::unique_ptr<LocalIntersectionOperatorType> local_op_; const XT::Common::Parameter param_; + const double scaling_; const real_t<F> eps_; DiscreteFunction<V, SGV, s_r, s_rC, F> source_; DiscreteFunction<V, RGV, r_r, r_rC, F> range_; diff --git a/dune/gdt/operators/matrix-based.hh b/dune/gdt/operators/matrix-based.hh index f66257a3faeaaa37a838ffcaa93f86baa3194713..1b8885fcd7ed60e3fdd3ed94782adac52f68ba46 100644 --- a/dune/gdt/operators/matrix-based.hh +++ b/dune/gdt/operators/matrix-based.hh @@ -242,6 +242,7 @@ public: using typename OperatorBaseType::SourceSpaceType; using typename OperatorBaseType::RangeSpaceType; using typename OperatorBaseType::MatrixOperatorType; + using typename OperatorBaseType::FieldType; using typename OperatorBaseType::V; using typename OperatorBaseType::F; @@ -265,6 +266,7 @@ public: , OperatorBaseType(source_spc, range_spc, MatrixStorage::access()) , WalkerBaseType(assembly_grid_view) , assembled_(false) + , scaling(1.) { // to detect assembly this->append( @@ -283,6 +285,7 @@ public: , OperatorBaseType(source_spc, range_spc, MatrixStorage::access()) , WalkerBaseType(assembly_grid_view) , assembled_(false) + , scaling(1.) { // to detect assembly this->append( @@ -296,6 +299,8 @@ public: return MatrixStorage::access(); } + FieldType scaling; + using WalkerBaseType::append; ThisType& append(const LocalElementBilinearFormInterface<E, r_r, r_rC, F, F, s_r, s_rC, F>& local_bilinear_form, @@ -303,8 +308,11 @@ public: const ElementFilterType& filter = ApplyOnAllElements()) { using LocalAssemblerType = LocalElementBilinearFormAssembler<M, SGV, r_r, r_rC, F, RGV, SGV, s_r, s_rC, F>; - this->append(new LocalAssemblerType( - this->range_space(), this->source_space(), local_bilinear_form, MatrixStorage::access(), param), + this->append(new LocalAssemblerType(this->range_space(), + this->source_space(), + local_bilinear_form, + MatrixStorage::access(), + param + XT::Common::Parameter("matrixoperator.scaling", scaling)), filter); return *this; } @@ -315,8 +323,11 @@ public: { using LocalAssemblerType = LocalIntersectionBilinearFormAssembler<MatrixType, AssemblyGridViewType, r_r, r_rC, F, RGV, SGV, s_r, s_rC, F>; - this->append(new LocalAssemblerType( - this->range_space(), this->source_space(), local_bilinear_form, MatrixStorage::access(), param), + this->append(new LocalAssemblerType(this->range_space(), + this->source_space(), + local_bilinear_form, + MatrixStorage::access(), + param + XT::Common::Parameter("matrixoperator.scaling", scaling)), filter); return *this; } // ... append(...) @@ -330,7 +341,12 @@ public: const ElementFilterType& filter = ApplyOnAllElements()) { this->append(new LocalElementOperatorFiniteDifferenceJacobianAssembler<M, SGV, s_r, s_rC, F, r_r, r_rC>( - this->source_space(), this->range_space(), MatrixStorage::access(), source, local_operator, param), + this->source_space(), + this->range_space(), + MatrixStorage::access(), + source, + local_operator, + param + XT::Common::Parameter("matrixoperator.scaling", scaling)), filter); return *this; } @@ -341,7 +357,12 @@ public: const IntersectionFilterType& filter = ApplyOnAllIntersections()) { this->append(new LocalIntersectionOperatorFiniteDifferenceJacobianAssembler<M, SGV, s_r, s_rC, F, r_r, r_rC>( - this->source_space(), this->range_space(), MatrixStorage::access(), source, local_operator, param), + this->source_space(), + this->range_space(), + MatrixStorage::access(), + source, + local_operator, + param + XT::Common::Parameter("matrixoperator.scaling", scaling)), filter); return *this; }