From: janssen Date: Sun, 28 Mar 2010 17:10:52 +0000 (+0000) Subject: the Simple classes now use constraints as well; some small modifications X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=427e127a4942eb7e2ec865b39188a8fc31a64c07;p=dealii-svn.git the Simple classes now use constraints as well; some small modifications git-svn-id: https://svn.dealii.org/trunk@20898 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h index abd9c729e1..4ec6764d87 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h @@ -265,6 +265,10 @@ namespace MeshWorker { public: void initialize(NamedData& results); + /** + * Initialize the constraints. + */ + void initialize(const ConstraintMatrix& constraints); /** * Initialize the local data * in the @@ -306,6 +310,10 @@ namespace MeshWorker * filled by assemble(). */ NamedData > > residuals; + /** + * A pointer to the object containing constraints. + */ + SmartPointer > constraints; }; /** @@ -394,16 +402,19 @@ namespace MeshWorker * stored as a vector of * pointers. */ - NamedData > > residuals; + NamedData > > residuals; /** * A pointer to the object containing the block structure. */ - SmartPointer block_info; + SmartPointer > block_info; /** * A pointer to the object containing constraints. */ - SmartPointer constraints; + SmartPointer > constraints; }; @@ -440,6 +451,10 @@ namespace MeshWorker * for later assembling. */ void initialize(MATRIX& m); + /** + * Initialize the constraints. + */ + void initialize(const ConstraintMatrix& constraints); /** * Initialize the local data * in the @@ -492,6 +507,10 @@ namespace MeshWorker * assembled. */ SmartPointer > matrix; + /** + * A pointer to the object containing constraints. + */ + SmartPointer > constraints; /** * The smallest positive @@ -753,12 +772,13 @@ namespace MeshWorker /** * A pointer to the object containing the block structure. */ - SmartPointer block_info; + SmartPointer > block_info; /** * A pointer to the object containing constraints. */ SmartPointer > constraints; + MatrixLocalBlocksToGlobalBlocks > constraints; /** * The smallest positive @@ -911,7 +931,7 @@ namespace MeshWorker /** * A pointer to the object containing the block structure. */ - SmartPointer block_info; + SmartPointer > block_info; /** * The smallest positive @@ -1114,6 +1134,13 @@ namespace MeshWorker residuals = results; } + template + inline void + ResidualSimple::initialize(const ConstraintMatrix& c) + { + constraints = &c; + } + template template @@ -1130,8 +1157,16 @@ namespace MeshWorker ResidualSimple::assemble(const DoFInfo& info) { for (unsigned int k=0;kdistribute_local_to_global( + info.vector(k).block(0), info.indices, (*residuals(k))); + } } @@ -1143,10 +1178,20 @@ namespace MeshWorker { for (unsigned int k=0;kdistribute_local_to_global( + info1.vector(k).block(0), info1.indices, (*residuals(k))); + constraints->distribute_local_to_global( + info2.vector(k).block(0), info2.indices, (*residuals(k))); + } } } @@ -1189,20 +1234,20 @@ namespace MeshWorker { if(constraints == 0) { - for (unsigned int b=0;bblock_info->local().local_to_global(b, j); - global(dof[jcell]) += local.block(b)(j); - } + for (unsigned int b=0;bblock_info->local().local_to_global(b, j); + global(dof[jcell]) += local.block(b)(j); + } } else constraints->distribute_local_to_global(local, dof, global); @@ -1251,6 +1296,14 @@ namespace MeshWorker { matrix = &m; } + + + template + inline void + MatrixSimple::initialize(const ConstraintMatrix& c) + { + constraints = &c; + } template @@ -1272,11 +1325,17 @@ namespace MeshWorker { AssertDimension(M.m(), i1.size()); AssertDimension(M.n(), i2.size()); - + + if(constraints == 0) + { for (unsigned int j=0; j= threshold) matrix->add(i1[j], i2[k], M(j,k)); + } + else + constraints->distribute_local_to_global( + M, i1, i2, *matrix); } diff --git a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h index 41f709fd65..2a01c5839e 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h @@ -68,6 +68,18 @@ namespace MeshWorker bool gradients = false, bool hessians = false); + /** + * Does the same as the function + * above but it is possible to + * select a block of the global + * vector. + */ +// void add(const std::string& name, +// const unsigned int selected_block, +// bool values = true, +// bool gradients = false, +// bool hessians = false); + /** * Initialize the selection * field with a data @@ -449,6 +461,17 @@ namespace MeshWorker } + //inline void + //VectorSelector::add(const std::string& name, + // const unsigned int block, + // bool values, bool gradients, bool hessians) + //{ + // if (values) value_selection.add(name, block); + // if (gradients) gradient_selection.add(name, block); + // if (hessians) hessian_selection.add(name, block); + //} + + template inline void VectorSelector::initialize(const NamedData& src) diff --git a/deal.II/deal.II/include/numerics/newton.h b/deal.II/deal.II/include/numerics/newton.h index 5192221d94..5d03cea5a4 100644 --- a/deal.II/deal.II/include/numerics/newton.h +++ b/deal.II/deal.II/include/numerics/newton.h @@ -82,6 +82,12 @@ namespace Algorithms * Read the parameters. */ void initialize (ParameterHandler& param); + + /** + * Initialize the pointer + * data_out for debugging. + */ + void initialize (OutputOperator& output); /** * The actual Newton @@ -116,13 +122,20 @@ namespace Algorithms /** * The operator computing the residual. */ - SmartPointer > residual; + SmartPointer, Newton > residual; /** * The operator applying the * inverse derivative to the residual. */ - SmartPointer > inverse_derivative; + SmartPointer, Newton > inverse_derivative; + + /** + * The operator handling the output + * in case the debug_vectors is true. + * Call the initialize function first. + */ + SmartPointer, Newton > data_out; /** * This flag is set by the @@ -166,6 +179,8 @@ namespace Algorithms * every Newton step. */ double assemble_threshold; + + public: /** * Print residual, update and * updated solution after each diff --git a/deal.II/deal.II/include/numerics/newton.templates.h b/deal.II/deal.II/include/numerics/newton.templates.h index 6f1e0173f9..3728c1fe8c 100644 --- a/deal.II/deal.II/include/numerics/newton.templates.h +++ b/deal.II/deal.II/include/numerics/newton.templates.h @@ -60,6 +60,12 @@ namespace Algorithms param.leave_subsection (); } + template + void + Newton::initialize (OutputOperator& output) + { + data_out = &output; + } template void @@ -125,6 +131,16 @@ namespace Algorithms << e.last_residual << std::endl; } + if (debug_vectors) + { + NamedData out; + VECTOR* p = &u; out.add(p, "solution"); + p = Du; out.add(p, "update"); + p = res; out.add(p, "residual"); + *data_out << step; + *data_out << out; + } + u.add(-1., *Du); old_residual = resnorm; (*residual)(out1, src1); @@ -147,19 +163,6 @@ namespace Algorithms (*residual)(out1, src1); resnorm = res->l2_norm(); } - - if (debug_vectors) - { - std::ostringstream streamOut; - streamOut << "Newton_" << std::setw(3) << std::setfill('0') << step; - std::string name = streamOut.str(); - NamedData out; - out.add(&u, "solution"); - out.add(Du, "update"); - out.add(res, "residual"); -//TODO: Implement! -// app->data_out(name, out); - } } deallog.pop(); diff --git a/deal.II/deal.II/include/numerics/operator.h b/deal.II/deal.II/include/numerics/operator.h index fd6d6ae980..cad77faf3b 100644 --- a/deal.II/deal.II/include/numerics/operator.h +++ b/deal.II/deal.II/include/numerics/operator.h @@ -16,6 +16,9 @@ #include #include #include +#include + +#include DEAL_II_NAMESPACE_OPEN @@ -94,11 +97,21 @@ namespace Algorithms template class OutputOperator : public Subscriptor { + OutputOperator(const OutputOperator&); public: + OutputOperator (); /** * Empty virtual destructor. */ virtual ~OutputOperator(); + /** + * Set the stream #os to + * which data is written. If + * no stream is selected with + * this function, data goes + * to #deallog. + */ + void initialize_stream(std::ostream& stream); /** * Set the current step. */ @@ -107,9 +120,11 @@ namespace Algorithms /** * Output all the vectors in NamedData. */ - virtual OutputOperator& operator<< (const NamedData vectors); + virtual OutputOperator& operator<< (const NamedData& vectors); protected: unsigned int step; + private: + std::ostream* os; }; template @@ -119,6 +134,27 @@ namespace Algorithms step = s; return *this; } + + template + class DoFOutputOperator : public OutputOperator + { + public: + void initialize (DoFHandler& dof_handler); + + virtual OutputOperator& + operator<<(const NamedData& vectors); + + private: + SmartPointer, + DoFOutputOperator > dof; + }; + + template + void DoFOutputOperator::initialize( + DoFHandler& dof_handler) + { + dof = &dof_handler; + } } diff --git a/deal.II/deal.II/include/numerics/operator.templates.h b/deal.II/deal.II/include/numerics/operator.templates.h index 20b7bc96da..a2ac9e57f6 100644 --- a/deal.II/deal.II/include/numerics/operator.templates.h +++ b/deal.II/deal.II/include/numerics/operator.templates.h @@ -11,6 +11,7 @@ //--------------------------------------------------------------------------- #include +#include DEAL_II_NAMESPACE_OPEN @@ -19,6 +20,7 @@ namespace Algorithms template Operator::~Operator() {} + template @@ -40,20 +42,58 @@ namespace Algorithms OutputOperator::~OutputOperator() {} + template + OutputOperator::OutputOperator() + : + os(0) + {} + + template + void OutputOperator::initialize_stream(std::ostream& stream) + { + os =&stream; + } template OutputOperator& - OutputOperator::operator<< (const NamedData vectors) + OutputOperator::operator<< (const NamedData& vectors) { - std::cout << ' ' << step; - for (unsigned int i=0;iprint(std::cout); - std::cout << std::endl; + if(os == 0) + { + //TODO: make this possible + //deallog << ' ' << step; + //for (unsigned int i=0;iprint(deallog); + //deallog << std::endl; + } + else + { + (*os) << ' ' << step; + for (unsigned int i=0;iprint(*os); + (*os) << std::endl; + } + return *this; + } + + template + OutputOperator& + DoFOutputOperator::operator<<( + const NamedData& vectors) + { + Assert ((dof!=0), ExcNotInitialized()); + DataOut out; + out.attach_dof_handler (*dof); + out.add_data_vector (*vectors(vectors.find("solution")), "solution"); + out.add_data_vector (*vectors(vectors.find("update")), "update"); + out.add_data_vector (*vectors(vectors.find("residual")), "residual"); + std::ostringstream streamOut; + streamOut << "Newton_" << std::setw(3) << std::setfill('0') << this->step; + std::ofstream out_filename (streamOut.str().c_str()); + out.build_patches (2); + out.write_gnuplot (out_filename); return *this; } - - - } DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/numerics/data_out.cc b/deal.II/deal.II/source/numerics/data_out.cc index f23595084e..668eeb2a6d 100644 --- a/deal.II/deal.II/source/numerics/data_out.cc +++ b/deal.II/deal.II/source/numerics/data_out.cc @@ -1130,7 +1130,9 @@ add_data_vector (const V &, \ INSTANTIATE(DH,D0,D1,D2,D3,Vector); \ INSTANTIATE(DH,D0,D1,D2,D3,Vector); \ INSTANTIATE(DH,D0,D1,D2,D3,BlockVector) ; \ - INSTANTIATE(DH,D0,D1,D2,D3,BlockVector) + INSTANTIATE(DH,D0,D1,D2,D3,BlockVector); \ + INSTANTIATE(DH,D0,D1,D2,D3,Vector); \ + INSTANTIATE(DH,D0,D1,D2,D3,BlockVector) #else # define INSTANTIATE_VECTORS(DH,D0,D1,D2,D3) \ INSTANTIATE(DH,D0,D1,D2,D3,Vector); \ diff --git a/deal.II/deal.II/source/numerics/operator.all_dimensions.cc b/deal.II/deal.II/source/numerics/operator.all_dimensions.cc index fe946a6256..90a137b58f 100644 --- a/deal.II/deal.II/source/numerics/operator.all_dimensions.cc +++ b/deal.II/deal.II/source/numerics/operator.all_dimensions.cc @@ -31,6 +31,6 @@ DEAL_II_NAMESPACE_OPEN using namespace Algorithms; -#include "operator.inst" +#include "operator.all_dimensions.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/numerics/operator.all_dimensions.inst.in b/deal.II/deal.II/source/numerics/operator.all_dimensions.inst.in new file mode 100644 index 0000000000..d4cf0f7536 --- /dev/null +++ b/deal.II/deal.II/source/numerics/operator.all_dimensions.inst.in @@ -0,0 +1,20 @@ +//--------------------------------------------------------------------------- +// $Id: vectors.inst.in 19046 2009-07-08 19:30:23Z bangerth $ +// Version: $Name$ +// +// Copyright (C) 2010 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +for (VEC : SERIAL_VECTORS) +{ + template class Operator; + template class OutputOperator; + template class Newton; + template class ThetaTimestepping; +} diff --git a/deal.II/deal.II/source/numerics/operator.cc b/deal.II/deal.II/source/numerics/operator.cc new file mode 100644 index 0000000000..fe946a6256 --- /dev/null +++ b/deal.II/deal.II/source/numerics/operator.cc @@ -0,0 +1,36 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2010 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +using namespace Algorithms; + +#include "operator.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/deal.II/source/numerics/operator.inst.in b/deal.II/deal.II/source/numerics/operator.inst.in index 7b5086980f..631ec8e983 100644 --- a/deal.II/deal.II/source/numerics/operator.inst.in +++ b/deal.II/deal.II/source/numerics/operator.inst.in @@ -13,7 +13,5 @@ for (VEC : SERIAL_VECTORS) { - template class Operator; - template class Newton; - template class ThetaTimestepping; + template class DoFOutputOperator; }