Chaste Commit::ca8ccdedf819b6e02855bc0e8e6f50bdecbc5208
AbstractContinuumMechanicsSolver.hpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#ifndef ABSTRACTCONTINUUMMECHANICSSOLVER_HPP_
37#define ABSTRACTCONTINUUMMECHANICSSOLVER_HPP_
38
39#include "ContinuumMechanicsProblemDefinition.hpp"
40#include "CompressibilityType.hpp"
41#include "LinearSystem.hpp"
42#include "OutputFileHandler.hpp"
43#include "ReplicatableVector.hpp"
44#include "AbstractTetrahedralMesh.hpp"
45#include "QuadraticMesh.hpp"
46#include "DistributedQuadraticMesh.hpp"
47#include "Warnings.hpp"
48#include "PetscException.hpp"
49#include "GaussianQuadratureRule.hpp"
50#include "PetscTools.hpp"
51#include "MechanicsEventHandler.hpp"
52#include "CommandLineArguments.hpp"
53#include "VtkMeshWriter.hpp"
54
55
61typedef enum _ApplyDirichletBcsType
62{
63 LINEAR_PROBLEM,
64 NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY,
65 NONLINEAR_PROBLEM_APPLY_TO_EVERYTHING
66} ApplyDirichletBcsType;
67
68//forward declarations
69template<unsigned DIM> class StressRecoveror;
70template<unsigned DIM> class VtkNonlinearElasticitySolutionWriter;
71
86template<unsigned DIM>
88{
89 friend class StressRecoveror<DIM>;
91
92protected:
98
101
104
106 std::string mOutputDirectory;
107
110
114 std::vector<c_vector<double,DIM> > mSpatialSolution;
115
117 std::vector<double> mPressureSolution;
118
119
126 std::vector<double> mCurrentSolution;
127
130
133
138 CompressibilityType mCompressibilityType;
139
149
153 unsigned mNumDofs;
154
161
183
188
193
198
203
208
209
235 void ApplyDirichletBoundaryConditions(ApplyDirichletBcsType type, bool symmetricProblem);
236
269 void AddIdentityBlockForDummyPressureVariables(ApplyDirichletBcsType type);
270
271
289
290public:
300 std::string outputDirectory,
301 CompressibilityType compressibilityType);
302
305
317 void WriteCurrentSpatialSolution(std::string fileName, std::string fileExtension, int counterToAppend=-1);
318
331 void WriteCurrentPressureSolution( int counterToAppend=-1);
332
338 void SetWriteOutput(bool writeOutput=true);
339
347 void CreateVtkOutput(std::string spatialSolutionName="Spatial solution");
348
353 std::vector<double>& rGetCurrentSolution()
354 {
355 return mCurrentSolution;
356 }
357
362 virtual std::vector<c_vector<double,DIM> >& rGetSpatialSolution()=0;
363
372 std::vector<double>& rGetPressures();
373};
374
375template<unsigned DIM>
378 std::string outputDirectory,
379 CompressibilityType compressibilityType)
380 : mrQuadMesh(rQuadMesh),
381 mrProblemDefinition(rProblemDefinition),
382 mOutputDirectory(outputDirectory),
383 mpOutputFileHandler(nullptr),
384 mpQuadratureRule(nullptr),
385 mpBoundaryQuadratureRule(nullptr),
386 mCompressibilityType(compressibilityType),
387 mResidualVector(nullptr),
388 mSystemLhsMatrix(nullptr),
389 mPreconditionMatrix(nullptr)
390{
391 assert(DIM==2 || DIM==3);
392
393 //Check that the mesh is Quadratic
394 QuadraticMesh<DIM>* p_quad_mesh = dynamic_cast<QuadraticMesh<DIM>* >(&rQuadMesh);
395 DistributedQuadraticMesh<DIM>* p_distributed_quad_mesh = dynamic_cast<DistributedQuadraticMesh<DIM>* >(&rQuadMesh);
396
397 if ((p_quad_mesh == nullptr) && (p_distributed_quad_mesh == nullptr))
398 {
399 EXCEPTION("Continuum mechanics solvers require a quadratic mesh");
400 }
401
402
403 mVerbose = (mrProblemDefinition.GetVerboseDuringSolve() ||
404 CommandLineArguments::Instance()->OptionExists("-mech_verbose") ||
405 CommandLineArguments::Instance()->OptionExists("-mech_very_verbose") );
406
408 if (mWriteOutput)
409 {
411 }
412
413 // See dox for mProblemDimension
414 mProblemDimension = mCompressibilityType==COMPRESSIBLE ? DIM : DIM+1;
416
418
419 // In general the Jacobian for a mechanics problem is non-polynomial.
420 // We therefore use the highest order integration rule available.
422 // The boundary forcing terms (or tractions) are also non-polynomial in general.
423 // Again, we use the highest order integration rule available.
425
426 mCurrentSolution.resize(mNumDofs, 0.0);
427}
428
429template<unsigned DIM>
431{
432 if (mpOutputFileHandler)
433 {
434 delete mpOutputFileHandler;
435 }
436
437 if (mpQuadratureRule)
438 {
439 delete mpQuadratureRule;
440 delete mpBoundaryQuadratureRule;
441 }
442
443 if (mResidualVector)
444 {
445 PetscTools::Destroy(mResidualVector);
446 PetscTools::Destroy(mLinearSystemRhsVector);
447 PetscTools::Destroy(mSystemLhsMatrix);
448 PetscTools::Destroy(mPreconditionMatrix);
449 }
450
451 if (mDirichletBoundaryConditionsVector)
452 {
453 PetscTools::Destroy(mDirichletBoundaryConditionsVector);
454 }
455}
456
457template<unsigned DIM>
459 std::string fileExtension,
460 int counterToAppend)
461{
462 // Only write output if the flag mWriteOutput has been set
463 if (!mWriteOutput)
464 {
465 return;
466 }
467
469 {
470 std::stringstream file_name;
471
472 file_name << fileName;
473 if (counterToAppend >= 0)
474 {
475 file_name << "_" << counterToAppend;
476 }
477 file_name << "." << fileExtension;
478
479 out_stream p_file = mpOutputFileHandler->OpenOutputFile(file_name.str());
480
481 std::vector<c_vector<double,DIM> >& r_spatial_solution = rGetSpatialSolution();
482 for (unsigned i=0; i<r_spatial_solution.size(); i++)
483 {
484 // for (unsigned j=0; j<DIM; j++)
485 // {
486 // *p_file << mrQuadMesh.GetNode(i)->rGetLocation()[j] << " ";
487 // }
488 for (unsigned j=0; j<DIM; j++)
489 {
490 *p_file << r_spatial_solution[i](j) << " ";
491 }
492 *p_file << "\n";
493 }
494 p_file->close();
495 }
496 PetscTools::Barrier("WriteSpatial");
497}
498
499template<unsigned DIM>
501{
502 // Only write output if the flag mWriteOutput has been set
503 if (mWriteOutput)
504 {
505 // Only the master writes
506 if (PetscTools::AmMaster() && mWriteOutput)
507 {
508 std::stringstream file_name;
509
510 file_name << "pressure";
511 if (counterToAppend >= 0)
512 {
513 file_name << "_" << counterToAppend;
514 }
515 file_name << ".txt";
516
517 out_stream p_file = mpOutputFileHandler->OpenOutputFile(file_name.str());
518
519 std::vector<double> &r_pressure = rGetPressures();
520 for (unsigned i = 0; i < r_pressure.size(); i++)
521 {
522 for (unsigned j = 0; j < DIM; j++)
523 {
524 *p_file << mrQuadMesh.GetNode(i)->rGetLocation()[j] << " ";
525 }
526
527 *p_file << r_pressure[i] << "\n";
528 }
529 p_file->close();
530 }
531 PetscTools::Barrier("WritePressure");
532 }
533}
534
535template<unsigned DIM>
537{
538 if (writeOutput && (mOutputDirectory==""))
539 {
540 EXCEPTION("Can't write output if no output directory was given in constructor");
541 }
542 mWriteOutput = writeOutput;
543}
544
545// ** TO BE DEPRECATED - see #2321 **
546template<unsigned DIM>
548{
549 if (this->mOutputDirectory=="")
550 {
551 EXCEPTION("No output directory was given so no output was written, cannot convert to VTK format");
552 }
553#ifdef CHASTE_VTK
554 VtkMeshWriter<DIM, DIM> mesh_writer(this->mOutputDirectory + "/vtk", "solution", true);
555
556 mesh_writer.AddPointData(spatialSolutionName, this->rGetSpatialSolution());
557
558 if (mCompressibilityType==INCOMPRESSIBLE)
559 {
560 mesh_writer.AddPointData("Pressure", rGetPressures());
561 }
562
563 //Output the element attribute as cell data.
564 std::vector<double> element_attribute;
565 for (typename QuadraticMesh<DIM>::ElementIterator iter = this->mrQuadMesh.GetElementIteratorBegin();
566 iter != this->mrQuadMesh.GetElementIteratorEnd();
567 ++iter)
568 {
569 element_attribute.push_back(iter->GetAttribute());
570 }
571 mesh_writer.AddCellData("Attribute", element_attribute);
572
573 mesh_writer.WriteFilesUsingMesh(this->mrQuadMesh);
574#endif
575}
576
577template<unsigned DIM>
579{
580 assert(mProblemDimension==DIM+1);
581
582 mPressureSolution.clear();
583 mPressureSolution.resize(mrQuadMesh.GetNumNodes());
584
585 for (unsigned i=0; i<mrQuadMesh.GetNumNodes(); i++)
586 {
587 mPressureSolution[i] = mCurrentSolution[mProblemDimension*i + DIM];
588 }
589 return mPressureSolution;
590}
591
592template<unsigned DIM>
594{
595 assert(mProblemDimension==DIM+1);
596
597 // For quadratic triangles, node 3 is between nodes 1 and 2, node 4 is between 0 and 2, etc
598 unsigned internal_nodes_2d[3] = {3,4,5};
599 unsigned neighbouring_vertices_2d[3][2] = { {1,2}, {2,0}, {0,1} };
600
601 // ordering for quadratic tetrahedra
602 unsigned internal_nodes_3d[6] = {4,5,6,7,8,9};
603 unsigned neighbouring_vertices_3d[6][2] = { {0,1}, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} };
604
605 unsigned num_internal_nodes_per_element = DIM==2 ? 3 : 6;
606
607 // loop over elements, then loop over edges.
608 for (typename AbstractTetrahedralMesh<DIM,DIM>::ElementIterator iter = mrQuadMesh.GetElementIteratorBegin();
609 iter != mrQuadMesh.GetElementIteratorEnd();
610 ++iter)
611 {
612 for (unsigned i=0; i<num_internal_nodes_per_element; i++)
613 {
614 unsigned global_index;
615 double left_val;
616 double right_val;
617
618 if (DIM == 2)
619 {
620 global_index = iter->GetNodeGlobalIndex( internal_nodes_2d[i] );
621 unsigned vertex_0_global_index =iter->GetNodeGlobalIndex( neighbouring_vertices_2d[i][0] );
622 unsigned vertex_1_global_index =iter->GetNodeGlobalIndex( neighbouring_vertices_2d[i][1] );
623 left_val = mCurrentSolution[mProblemDimension*vertex_0_global_index + DIM];
624 right_val = mCurrentSolution[mProblemDimension*vertex_1_global_index + DIM];
625 }
626 else
627 {
628 global_index = iter->GetNodeGlobalIndex( internal_nodes_3d[i] );
629 unsigned vertex_0_global_index =iter->GetNodeGlobalIndex( neighbouring_vertices_3d[i][0] );
630 unsigned vertex_1_global_index =iter->GetNodeGlobalIndex( neighbouring_vertices_3d[i][1] );
631 left_val = mCurrentSolution[mProblemDimension*vertex_0_global_index + DIM];
632 right_val = mCurrentSolution[mProblemDimension*vertex_1_global_index + DIM];
633 }
634
635 // this line assumes the internal node is midway between the two vertices
636 mCurrentSolution[mProblemDimension*global_index + DIM] = 0.5 * (left_val + right_val);
637 }
638 }
639}
640
641/*
642 * This method applies the appropriate BCs to the residual and/or linear system
643 *
644 * For the latter, and second input parameter==true, the BCs are imposed in such a way as
645 * to ensure that a symmetric linear system remains symmetric. For each row with boundary condition
646 * applied, both the row and column are zero'd and the RHS vector modified to take into account the
647 * zero'd column.
648 *
649 * Suppose we have a matrix
650 * [a b c] [x] = [ b1 ]
651 * [d e f] [y] [ b2 ]
652 * [g h i] [z] [ b3 ]
653 * and we want to apply the boundary condition x=v without losing symmetry if the matrix is
654 * symmetric. We apply the boundary condition
655 * [1 0 0] [x] = [ v ]
656 * [d e f] [y] [ b2 ]
657 * [g h i] [z] [ b3 ]
658 * and then zero the column as well, adding a term to the RHS to take account for the
659 * zero-matrix components
660 * [1 0 0] [x] = [ v ] - v[ 0 ]
661 * [0 e f] [y] [ b2 ] [ d ]
662 * [0 h i] [z] [ b3 ] [ g ]
663 * Note the last term is the first column of the matrix, with one component zeroed, and
664 * multiplied by the boundary condition value. This last term is then stored in
665 * mDirichletBoundaryConditionsVector, and in general is equal to:
666 * SUM_{d=1..D} v_d a'_d
667 * where v_d is the boundary value of boundary condition d (d an index into the matrix),
668 * and a'_d is the dth-column of the matrix but with the d-th component zeroed, and where
669 * there are D boundary conditions
670 */
671template<unsigned DIM>
672void AbstractContinuumMechanicsSolver<DIM>::ApplyDirichletBoundaryConditions(ApplyDirichletBcsType type, bool symmetricProblem)
673{
674 std::vector<unsigned> rows;
675 std::vector<double> values;
676
677 // Whether to apply symmetrically, ie alter columns as well as rows (see comment above)
678 bool applySymmetrically = (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY) && symmetricProblem;
679
680 if (applySymmetrically)
681 {
682 if (mDirichletBoundaryConditionsVector == nullptr)
683 {
684 VecDuplicate(mResidualVector, &mDirichletBoundaryConditionsVector);
685 }
686
687 PetscVecTools::Zero(mDirichletBoundaryConditionsVector);
688 PetscMatTools::Finalise(mSystemLhsMatrix);
689 }
690
692 // collect the entries to be altered
694
695 for (unsigned i=0; i<mrProblemDefinition.rGetDirichletNodes().size(); i++)
696 {
697 unsigned node_index = mrProblemDefinition.rGetDirichletNodes()[i];
698
699 for (unsigned j=0; j<DIM; j++)
700 {
701 double dirichlet_val = mrProblemDefinition.rGetDirichletNodeValues()[i](j);
702
704 {
705 double val;
706 unsigned dof_index = mProblemDimension*node_index+j;
707
708 if (type == LINEAR_PROBLEM)
709 {
710 val = dirichlet_val;
711 }
712 else
713 {
714 // The boundary conditions on the NONLINEAR SYSTEM are x=boundary_values
715 // on the boundary nodes. However:
716 // The boundary conditions on the LINEAR SYSTEM at each Newton step (Ju=f,
717 // where J is the Jacobian, u the negative update vector and f is the residual) is:
718 // u=current_soln-boundary_values on the boundary nodes
719 val = mCurrentSolution[dof_index] - dirichlet_val;
720 }
721 rows.push_back(dof_index);
722 values.push_back(val);
723 }
724 }
725 }
726
728 // do the alterations
730
731 if (applySymmetrically)
732 {
733 // Modify the matrix columns
734 for (unsigned i=0; i<rows.size(); i++)
735 {
736 unsigned col = rows[i];
737 double minus_value = -values[i];
738
739 // Get a vector which will store the column of the matrix (column d, where d is
740 // the index of the row (and column) to be altered for the boundary condition.
741 // Since the matrix is symmetric when get row number "col" and treat it as a column.
742 // PETSc uses compressed row format and therefore getting rows is far more efficient
743 // than getting columns.
744 Vec matrix_col = PetscMatTools::GetMatrixRowDistributed(mSystemLhsMatrix,col);
745
746 // Zero the correct entry of the column
747 PetscVecTools::SetElement(matrix_col, col, 0.0);
748
749 // Set up the RHS Dirichlet boundary conditions vector
750 // E.g. for a boundary node at the zeroth node (x_0 = value), this is equal to
751 // -value*[0 a_21 a_31 .. a_N1]
752 // and will be added to the RHS.
753 PetscVecTools::AddScaledVector(mDirichletBoundaryConditionsVector, matrix_col, minus_value);
754 PetscTools::Destroy(matrix_col);
755 }
756 }
757
758 if (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY) // i.e. doing a whole linear system
759 {
760 // Now zero the appropriate rows and columns of the matrix. If the matrix is symmetric we apply the
761 // boundary conditions in a way the symmetry isn't lost (rows and columns). If not only the row is
762 // zeroed.
763 if (applySymmetrically)
764 {
765 PetscMatTools::ZeroRowsAndColumnsWithValueOnDiagonal(mSystemLhsMatrix, rows, 1.0);
766 PetscMatTools::ZeroRowsAndColumnsWithValueOnDiagonal(mPreconditionMatrix, rows, 1.0);
767
768 // Apply the RHS boundary conditions modification if required.
769 // (AddScaledVector now needs assembly to be finalised.)
770 PetscVecTools::Finalise(mLinearSystemRhsVector);
771 PetscVecTools::AddScaledVector(mLinearSystemRhsVector, mDirichletBoundaryConditionsVector, 1.0);
772 }
773 else
774 {
775 PetscMatTools::ZeroRowsWithValueOnDiagonal(mSystemLhsMatrix, rows, 1.0);
776 PetscMatTools::ZeroRowsWithValueOnDiagonal(mPreconditionMatrix, rows, 1.0);
777 }
778 }
779
780 if (type!=LINEAR_PROBLEM)
781 {
782 for (unsigned i=0; i<rows.size(); i++)
783 {
784 PetscVecTools::SetElement(mResidualVector, rows[i], values[i]);
785 }
786 }
787
788 if (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY)
789 {
790 for (unsigned i=0; i<rows.size(); i++)
791 {
792 PetscVecTools::SetElement(mLinearSystemRhsVector, rows[i], values[i]);
793 }
794 }
795}
796
797template<unsigned DIM>
799{
800 assert(mCompressibilityType==INCOMPRESSIBLE);
801
802 int lo, hi;
803 VecGetOwnershipRange(mResidualVector, &lo, &hi);
804
805 for (unsigned i=0; i<mrQuadMesh.GetNumNodes(); i++)
806 {
807 if (mrQuadMesh.GetNode(i)->IsInternal())
808 {
809 unsigned row = (DIM+1)*i + DIM; // DIM+1 is the problem dimension
810 if (lo <= (int)row && (int)row < hi)
811 {
812 if (type!=LINEAR_PROBLEM)
813 {
814 PetscVecTools::SetElement(mResidualVector, row, mCurrentSolution[row]-0.0);
815 }
816 if (type!=NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY) // ie doing a whole linear system
817 {
818 double rhs_vector_val = type==LINEAR_PROBLEM ? 0.0 : mCurrentSolution[row]-0.0;
819 PetscVecTools::SetElement(mLinearSystemRhsVector, row, rhs_vector_val);
820 // This assumes the row is already zero, which is should be..
821 PetscMatTools::SetElement(mSystemLhsMatrix, row, row, 1.0);
822 PetscMatTools::SetElement(mPreconditionMatrix, row, row, 1.0);
823 }
824 }
825 }
826 }
827}
828
829template<unsigned DIM>
831{
832 Vec template_vec = mrQuadMesh.GetDistributedVectorFactory()->CreateVec(mProblemDimension);
833
835 // three vectors
837 VecDuplicate(template_vec, &mResidualVector);
838 VecDuplicate(mResidualVector, &mLinearSystemRhsVector);
839 // the one is only allocated if it will be needed (in ApplyDirichletBoundaryConditions),
840 // depending on whether the matrix is kept symmetric.
841 mDirichletBoundaryConditionsVector = nullptr;
842 PetscTools::Destroy(template_vec);
843
845 // two matrices
847
848 int lo, hi;
849 VecGetOwnershipRange(mResidualVector, &lo, &hi);
850 PetscInt local_size = hi - lo;
851
852
853 if (DIM==2)
854 {
855 // 2D: N elements around a point => 7N+3 non-zeros in that row? Assume N<=10 (structured mesh would have N_max=6) => 73.
856 unsigned num_non_zeros = std::min(75u, mNumDofs);
857
858 PetscTools::SetupMat(mSystemLhsMatrix, mNumDofs, mNumDofs, num_non_zeros, local_size, local_size);
859 PetscTools::SetupMat(mPreconditionMatrix, mNumDofs, mNumDofs, num_non_zeros, local_size, local_size);
860 }
861 else
862 {
863 assert(DIM==3);
864
865 // in 3d we get the number of containing elements for each node and use that to obtain an upper bound
866 // for the number of non-zeros for each DOF associated with that node.
867
868 int* num_non_zeros_each_row = new int[mNumDofs];
869 for (unsigned i=0; i<mNumDofs; i++)
870 {
871 num_non_zeros_each_row[i] = 0;
872 }
873
874 for (typename AbstractMesh<DIM,DIM>::NodeIterator iter = mrQuadMesh.GetNodeIteratorBegin();
875 iter != mrQuadMesh.GetNodeIteratorEnd();
876 ++iter)
877 {
878 // this upper bound neglects the fact that two containing elements will share the same nodes..
879 // 4 = max num dofs associated with this node
880 // 30 = 3*9+3 = 3 dimensions x 9 other nodes on this element + 3 vertices with a pressure unknown
881 unsigned num_non_zeros_upper_bound = 4 + 30*iter->GetNumContainingElements();
882
883 num_non_zeros_upper_bound = std::min(num_non_zeros_upper_bound, mNumDofs);
884
885 unsigned i = iter->GetIndex();
886
887 num_non_zeros_each_row[mProblemDimension*i + 0] = num_non_zeros_upper_bound;
888 num_non_zeros_each_row[mProblemDimension*i + 1] = num_non_zeros_upper_bound;
889 num_non_zeros_each_row[mProblemDimension*i + 2] = num_non_zeros_upper_bound;
890
891 if (mCompressibilityType==INCOMPRESSIBLE)
892 {
893 if (!iter->IsInternal())
894 {
895 num_non_zeros_each_row[mProblemDimension*i + 3] = num_non_zeros_upper_bound;
896 }
897 else
898 {
899 num_non_zeros_each_row[mProblemDimension*i + 3] = 1;
900 }
901 }
902 }
903
904 // NOTE: PetscTools::SetupMat() or the below creates a MATAIJ matrix, which means the matrix will
905 // be of type MATSEQAIJ if num_procs=1 and MATMPIAIJ otherwise. In the former case
906 // MatSeqAIJSetPreallocation MUST be called [MatMPIAIJSetPreallocation will have
907 // no effect (silently)], and vice versa in the latter case
908
911 //PetscTools::SetupMat(mSystemLhsMatrix, mNumDofs, mNumDofs, 0, PETSC_DECIDE, PETSC_DECIDE);
912 //PetscTools::SetupMat(mPreconditionMatrix, mNumDofs, mNumDofs, 0, PETSC_DECIDE, PETSC_DECIDE);
914
915 // possible todo: create a PetscTools::SetupMatNoAllocation()
916
917
918#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
919 MatCreate(PETSC_COMM_WORLD,local_size,local_size,mNumDofs,mNumDofs,&mSystemLhsMatrix);
920 MatCreate(PETSC_COMM_WORLD,local_size,local_size,mNumDofs,mNumDofs,&mPreconditionMatrix);
921#else //New API
922 MatCreate(PETSC_COMM_WORLD,&mSystemLhsMatrix);
923 MatCreate(PETSC_COMM_WORLD,&mPreconditionMatrix);
924 MatSetSizes(mSystemLhsMatrix,local_size,local_size,mNumDofs,mNumDofs);
925 MatSetSizes(mPreconditionMatrix,local_size,local_size,mNumDofs,mNumDofs);
926#endif
927
929 {
930 MatSetType(mSystemLhsMatrix, MATSEQAIJ);
931 MatSetType(mPreconditionMatrix, MATSEQAIJ);
932 MatSeqAIJSetPreallocation(mSystemLhsMatrix, PETSC_DEFAULT, num_non_zeros_each_row);
933 MatSeqAIJSetPreallocation(mPreconditionMatrix, PETSC_DEFAULT, num_non_zeros_each_row);
934 }
935 else
936 {
937 int* num_non_zeros_each_row_in_diag = new int[local_size];
938 int* num_non_zeros_each_row_off_diag = new int[local_size];
939 for (unsigned i=0; i<unsigned(local_size); i++)
940 {
941 num_non_zeros_each_row_in_diag[i] = num_non_zeros_each_row[lo+i];
942 num_non_zeros_each_row_off_diag[i] = num_non_zeros_each_row[lo+i];
943 // In the on process ("diagonal block") there cannot be more non-zero columns specified than there are rows
944 if (num_non_zeros_each_row_in_diag[i] > local_size)
945 {
946 num_non_zeros_each_row_in_diag[i] = local_size;
947 }
948 }
949
950 MatSetType(mSystemLhsMatrix, MATMPIAIJ);
951 MatSetType(mPreconditionMatrix, MATMPIAIJ);
952 MatMPIAIJSetPreallocation(mSystemLhsMatrix, PETSC_DEFAULT, num_non_zeros_each_row_in_diag, PETSC_DEFAULT, num_non_zeros_each_row_off_diag);
953 MatMPIAIJSetPreallocation(mPreconditionMatrix, PETSC_DEFAULT, num_non_zeros_each_row_in_diag, PETSC_DEFAULT, num_non_zeros_each_row_off_diag);
954 }
955
956 MatSetFromOptions(mSystemLhsMatrix);
957 MatSetFromOptions(mPreconditionMatrix);
958#if (PETSC_VERSION_MAJOR == 3) //PETSc 3.x.x
959 MatSetOption(mSystemLhsMatrix, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
960 MatSetOption(mPreconditionMatrix, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
961#else
962 MatSetOption(mSystemLhsMatrix, MAT_IGNORE_OFF_PROC_ENTRIES);
963 MatSetOption(mPreconditionMatrix, MAT_IGNORE_OFF_PROC_ENTRIES);
964#endif
965
966 //unsigned total_non_zeros = 0;
967 //for (unsigned i=0; i<mNumDofs; i++)
968 //{
969 // total_non_zeros += num_non_zeros_each_row[i];
970 //}
971 //std::cout << total_non_zeros << " versus " << 500*mNumDofs << "\n" << std::flush;
972
973 delete [] num_non_zeros_each_row;
974 }
975}
976#endif // ABSTRACTCONTINUUMMECHANICSSOLVER_HPP_
#define EXCEPTION(message)
std::vector< c_vector< double, DIM > > mSpatialSolution
AbstractTetrahedralMesh< DIM, DIM > & mrQuadMesh
GaussianQuadratureRule< DIM-1 > * mpBoundaryQuadratureRule
AbstractContinuumMechanicsSolver(AbstractTetrahedralMesh< DIM, DIM > &rQuadMesh, ContinuumMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory, CompressibilityType compressibilityType)
void WriteCurrentPressureSolution(int counterToAppend=-1)
void WriteCurrentSpatialSolution(std::string fileName, std::string fileExtension, int counterToAppend=-1)
virtual std::vector< c_vector< double, DIM > > & rGetSpatialSolution()=0
void AddIdentityBlockForDummyPressureVariables(ApplyDirichletBcsType type)
ContinuumMechanicsProblemDefinition< DIM > & mrProblemDefinition
void CreateVtkOutput(std::string spatialSolutionName="Spatial solution")
void ApplyDirichletBoundaryConditions(ApplyDirichletBcsType type, bool symmetricProblem)
virtual unsigned GetNumNodes() const
bool OptionExists(const std::string &rOption)
static CommandLineArguments * Instance()
static void ZeroRowsWithValueOnDiagonal(Mat matrix, std::vector< unsigned > &rRows, double diagonalValue)
static void ZeroRowsAndColumnsWithValueOnDiagonal(Mat matrix, std::vector< unsigned > rowColIndices, double diagonalValue)
static Vec GetMatrixRowDistributed(Mat matrix, unsigned rowIndex)
static void Finalise(Mat matrix)
static void SetElement(Mat matrix, PetscInt row, PetscInt col, double value)
static void Destroy(Vec &rVec)
static bool AmMaster()
static void Barrier(const std::string callerId="")
static bool IsSequential()
static void SetupMat(Mat &rMat, int numRows, int numColumns, unsigned rowPreallocation, int numLocalRows=PETSC_DECIDE, int numLocalColumns=PETSC_DECIDE, bool ignoreOffProcEntries=true, bool newAllocationError=true)
static void AddScaledVector(Vec y, Vec x, double scaleFactor)
static void Finalise(Vec vector)
static void SetElement(Vec vector, PetscInt row, double value)
static void Zero(Vec vector)
void AddPointData(std::string name, std::vector< double > data)
void AddCellData(std::string name, std::vector< double > data)
void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)