Veselin Dobrev
27ffe35369
general/globals.[ch]. Renamed class WrappedOStream to OutStream and use a different implementation in order to support the case when the OutStream is passed as a parameter of type std::ostream. In class MPI_Session, added an MPI_Comm parameter to both constructors, with a default value of MPI_COMM_WORLD. Hide the global variable MFEM_COMM_WORLD and allow access (read and write) with global functions: GetGlobalMPI_Comm() and SetGlobalMPI_Comm() - in the future, we may want to perform some tasks when changing the "global" MPI communicator. In mfem_error(), use GetGlobalMPI_Comm() instead of MPI_COMM_WORLD. Minor formatting changes and #include reordering.
150 lines
3.5 KiB
C++
150 lines
3.5 KiB
C++
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
|
|
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
|
|
// reserved. See file COPYRIGHT for details.
|
|
//
|
|
// This file is part of the MFEM library. For more information and source code
|
|
// availability see http://mfem.org.
|
|
//
|
|
// MFEM is free software; you can redistribute it and/or modify it under the
|
|
// terms of the GNU Lesser General Public License (as published by the Free
|
|
// Software Foundation) version 2.1 dated February 1999.
|
|
|
|
#include "vector.hpp"
|
|
#include "operator.hpp"
|
|
|
|
#include <iostream>
|
|
#include <iomanip>
|
|
|
|
namespace mfem
|
|
{
|
|
|
|
void Operator::FormLinearSystem(const Array<int> &ess_tdof_list,
|
|
Vector &x, Vector &b,
|
|
Operator* &Aout, Vector &X, Vector &B,
|
|
int copy_interior)
|
|
{
|
|
const Operator *P = this->GetProlongation();
|
|
const Operator *R = this->GetRestriction();
|
|
Operator *rap;
|
|
|
|
if (P)
|
|
{
|
|
// Variational restriction with P
|
|
B.SetSize(P->Width());
|
|
P->MultTranspose(b, B);
|
|
X.SetSize(R->Height());
|
|
R->Mult(x, X);
|
|
rap = new RAPOperator(*P, *this, *P);
|
|
}
|
|
else
|
|
{
|
|
// rap, X and B point to the same data as this, x and b
|
|
X.NewDataAndSize(x.GetData(), x.Size());
|
|
B.NewDataAndSize(b.GetData(), b.Size());
|
|
rap = this;
|
|
}
|
|
|
|
if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
|
|
|
|
// Impose the boundary conditions through a ConstrainedOperator, which owns
|
|
// the rap operator when P and R are non-trivial
|
|
ConstrainedOperator *A = new ConstrainedOperator(rap, ess_tdof_list,
|
|
rap != this);
|
|
A->EliminateRHS(X, B);
|
|
Aout = A;
|
|
}
|
|
|
|
void Operator::RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
|
|
{
|
|
const Operator *P = this->GetProlongation();
|
|
if (P)
|
|
{
|
|
// Apply conforming prolongation
|
|
x.SetSize(P->Height());
|
|
P->Mult(X, x);
|
|
}
|
|
else
|
|
{
|
|
// X and x point to the same data
|
|
}
|
|
}
|
|
|
|
void Operator::PrintMatlab(std::ostream & out, int n, int m) const
|
|
{
|
|
using namespace std;
|
|
if (n == 0) { n = width; }
|
|
if (m == 0) { m = height; }
|
|
|
|
Vector x(n), y(m);
|
|
x = 0.0;
|
|
|
|
out << setiosflags(ios::scientific | ios::showpos);
|
|
for (int i = 0; i < n; i++)
|
|
{
|
|
x(i) = 1.0;
|
|
Mult(x, y);
|
|
for (int j = 0; j < m; j++)
|
|
{
|
|
if (y(j))
|
|
{
|
|
out << j+1 << " " << i+1 << " " << y(j) << '\n';
|
|
}
|
|
}
|
|
x(i) = 0.0;
|
|
}
|
|
}
|
|
|
|
|
|
ConstrainedOperator::ConstrainedOperator(Operator *A, const Array<int> &list,
|
|
bool _own_A)
|
|
: Operator(A->Height(), A->Width()), A(A), own_A(_own_A)
|
|
{
|
|
constraint_list.MakeRef(list);
|
|
z.SetSize(height);
|
|
w.SetSize(height);
|
|
}
|
|
|
|
void ConstrainedOperator::EliminateRHS(const Vector &x, Vector &b) const
|
|
{
|
|
w = 0.0;
|
|
|
|
for (int i = 0; i < constraint_list.Size(); i++)
|
|
{
|
|
w(constraint_list[i]) = x(constraint_list[i]);
|
|
}
|
|
|
|
A->Mult(w, z);
|
|
|
|
b -= z;
|
|
|
|
for (int i = 0; i < constraint_list.Size(); i++)
|
|
{
|
|
b(constraint_list[i]) = x(constraint_list[i]);
|
|
}
|
|
}
|
|
|
|
void ConstrainedOperator::Mult(const Vector &x, Vector &y) const
|
|
{
|
|
if (constraint_list.Size() == 0)
|
|
{
|
|
A->Mult(x, y);
|
|
return;
|
|
}
|
|
|
|
z = x;
|
|
|
|
for (int i = 0; i < constraint_list.Size(); i++)
|
|
{
|
|
z(constraint_list[i]) = 0.0;
|
|
}
|
|
|
|
A->Mult(z, y);
|
|
|
|
for (int i = 0; i < constraint_list.Size(); i++)
|
|
{
|
|
y(constraint_list[i]) = x(constraint_list[i]);
|
|
}
|
|
}
|
|
|
|
}
|