#include <sstream>
#include <boost/lexical_cast.hpp>
#include "ProblemVec.h"
#include "RecoveryEstimator.h"
#include "Serializer.h"
#include "AbstractFunction.h"
#include "Operator.h"
#include "SystemVector.h"
#include "DOFMatrix.h"
#include "FiniteElemSpace.h"
#include "Estimator.h"
#include "Marker.h"
#include "AdaptInfo.h"
#include "FileWriter.h"
#include "CoarseningManager.h"
#include "RefinementManager.h"
#include "DualTraverse.h"
#include "Mesh.h"
#include "OEMSolver.h"
#include "DirichletBC.h"
#include "RobinBC.h"
#include "PeriodicBC.h"
#include "Lagrange.h"
#include "Flag.h"
#include "TraverseParallel.h"
#include "VtkWriter.h"
#include "ValueReader.h"
#include "ProblemVecDbg.h"
#include "Debug.h"
namespace AMDiS {
using boost::lexical_cast;
void ProblemVec::initialize(Flag initFlag,
ProblemVec *adoptProblem,
Flag adoptFlag)
{
FUNCNAME("ProblemVec::initialize()");
// === create meshes ===
if (meshes.size() != 0) {
WARNING("meshes already created\n");
} else {
if (initFlag.isSet(CREATE_MESH) ||
(!adoptFlag.isSet(INIT_MESH) &&
(initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE))))
createMesh();
if (adoptProblem &&
(adoptFlag.isSet(INIT_MESH) ||
adoptFlag.isSet(INIT_SYSTEM) ||
adoptFlag.isSet(INIT_FE_SPACE))) {
meshes = adoptProblem->getMeshes();
componentMeshes = adoptProblem->componentMeshes;
refinementManager = adoptProblem->refinementManager;
coarseningManager = adoptProblem->coarseningManager;
// If the adopt problem has fewer components than this problem, but only one
// mesh for all component, than scal up the componentMeshes array.
if (adoptProblem->getNumComponents() < nComponents) {
TEST_EXIT(meshes.size() == 1)("Daran muss ich noch arbeiten!\n");
componentMeshes.resize(nComponents);
for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
componentMeshes[i] = componentMeshes[0];
}
}
}
if (meshes.size() == 0)
WARNING("no mesh created\n");
// === create fespace ===
if (feSpaces.size() != 0) {
WARNING("feSpaces already created\n");
} else {
if (initFlag.isSet(INIT_FE_SPACE) ||
(initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE)))
createFeSpace(NULL);
if (adoptProblem &&
(adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
feSpaces = adoptProblem->getFeSpaces();
componentSpaces = adoptProblem->componentSpaces;
traverseInfo = adoptProblem->traverseInfo;
// If the adopt problem has fewer components than this problem, but only one
// fe space for all component, than scal up the componentSpaces array.
if (adoptProblem->getNumComponents() < nComponents) {
TEST_EXIT(feSpaces.size() == 1)("Daran muss ich noch arbeiten!\n");
componentSpaces.resize(nComponents);
for (int i = adoptProblem->getNumComponents(); i < nComponents; i++)
componentSpaces[i] = componentSpaces[0];
}
}
}
if (feSpaces.size() == 0)
WARNING("no feSpace created\n");
// === create system ===
if (initFlag.isSet(INIT_SYSTEM))
createMatricesAndVectors();
if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
solution = adoptProblem->getSolution();
rhs = adoptProblem->getRhs();
systemMatrix = adoptProblem->getSystemMatrix();
}
// === create solver ===
if (solver) {
WARNING("solver already created\n");
} else {
if (initFlag.isSet(INIT_SOLVER)) {
createSolver();
}
if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
TEST_EXIT(!solver)("solver already created\n");
solver = adoptProblem->getSolver();
}
}
if (!solver)
WARNING("no solver created\n");
// === create estimator ===
if (initFlag.isSet(INIT_ESTIMATOR))
createEstimator();
if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR))
estimator = adoptProblem->getEstimators();
// === create marker ===
if (initFlag.isSet(INIT_MARKER))
createMarker();
if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
marker = adoptProblem->getMarkers();
// === create file writer ===
if (initFlag.isSet(INIT_FILEWRITER))
createFileWriter();
// === read serialization and init mesh ===
// There are two possiblities where the user can define a serialization
// to be read from disk. Either by providing the parameter -rs when executing
// the program or in the init file. The -rs parameter is always checked first,
// because it can be added automatically when rescheduling the program
// before timeout of the runqueue.
int readSerialization = 0;
std::string serializationFilename = "";
GET_PARAMETER(0, "argv->rs", &serializationFilename);
// If the parameter -rs is set, we do nothing here, because the problem will be
// deserialized in the constructor of a following AdaptInstationary initialization.
if (!serializationFilename.compare("")) {
int readSerializationWithAdaptInfo = 0;
GET_PARAMETER(0, name + "->input->read serialization", "%d", &readSerialization);
GET_PARAMETER(0, name + "->input->serialization with adaptinfo", "%d",
&readSerializationWithAdaptInfo);
// The serialization file is only read, if the adaptInfo part should not be used.
// If the adaptInfo part should be also read, the serialization file will be read
// in the constructor of the AdaptInstationary problem, because we do not have here
// the adaptInfo object.
if (readSerialization && !readSerializationWithAdaptInfo) {
GET_PARAMETER(0, name + "->input->serialization filename",
&serializationFilename);
TEST_EXIT(serializationFilename != "")("no serialization file\n");
// If AMDiS is compiled for parallel computations, the deserialization is
// controlled by the parallel problem object.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
std::ifstream in(serializationFilename.c_str());
deserialize(in);
in.close();
MSG("Deserialization from file: %s\n", serializationFilename.c_str());
#endif
} else {
int globalRefinements = 0;
// If AMDiS is compied for parallel computations, the global refinements are
// ignored here. Later, each rank will add the global refinements to its
// private mesh.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
GET_PARAMETER(0, meshes[0]->getName() + "->global refinements", "%d",
&globalRefinements);
#endif
// Initialize the meshes if there is no serialization file.
for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
bool initMesh = initFlag.isSet(INIT_MESH) ||
(adoptProblem && adoptFlag.isSet(INIT_MESH));
if (initMesh && meshes[i] && !(meshes[i]->isInitialized()))
meshes[i]->initialize();
}
// === read value file and use it for the mesh values ===
std::string valueFilename("");
GET_PARAMETER(0, meshes[0]->getName() + "->value file name", &valueFilename);
if (valueFilename.length()) {
ValueReader::readValue(valueFilename,
meshes[0],
solution->getDOFVector(0),
meshes[0]->getMacroFileInfo());
meshes[0]->clearMacroFileInfo();
}
// === do global refinements ===
for (unsigned int i = 0; i < meshes.size(); i++)
if (initFlag.isSet(INIT_MESH) && meshes[i])
refinementManager->globalRefine(meshes[i], globalRefinements);
}
}
doOtherStuff();
}
void ProblemVec::createMesh()
{
FUNCNAME("ProblemVec::createMesh()");
componentMeshes.resize(nComponents);
std::map<int, Mesh*> meshForRefinementSet;
std::string meshName("");
GET_PARAMETER(0, name + "->mesh", &meshName);
TEST_EXIT(meshName != "")("no mesh name specified\n");
int dim = 0;
GET_PARAMETER(0, name + "->dim", "%d", &dim);
TEST_EXIT(dim)("no problem dimension specified!\n");
for (int i = 0; i < nComponents; i++) {
int refSet = -1;
GET_PARAMETER(0, name + "->refinement set[" +
lexical_cast<std::string>(i) + "]", "%d", &refSet);
if (refSet < 0)
refSet = 0;
if (meshForRefinementSet[refSet] == NULL) {
Mesh *newMesh = new Mesh(meshName, dim);
meshForRefinementSet[refSet] = newMesh;
meshes.push_back(newMesh);
nMeshes++;
}
componentMeshes[i] = meshForRefinementSet[refSet];
}
switch(dim) {
case 1:
coarseningManager = new CoarseningManager1d();
refinementManager = new RefinementManager1d();
break;
case 2:
coarseningManager = new CoarseningManager2d();
refinementManager = new RefinementManager2d();
break;
case 3:
coarseningManager = new CoarseningManager3d();
refinementManager = new RefinementManager3d();
break;
default:
ERROR_EXIT("invalid dim!\n");
}
}
void ProblemVec::createFeSpace(DOFAdmin *admin)
{
FUNCNAME("ProblemVec::createFeSpace()");
std::map<std::pair<Mesh*, int>, FiniteElemSpace*> feSpaceMap;
int dim = -1;
GET_PARAMETER(0, name + "->dim", "%d", &dim);
TEST_EXIT(dim != -1)("no problem dimension specified!\n");
componentSpaces.resize(nComponents, NULL);
traverseInfo.resize(nComponents);
for (int i = 0; i < nComponents; i++) {
int degree = 1;
GET_PARAMETER(0, name + "->polynomial degree[" +
boost::lexical_cast<std::string>(i) + "]","%d", °ree);
TEST_EXIT(componentSpaces[i] == NULL)("feSpace already created\n");
if (feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] == NULL) {
stringstream s;
s << name << "->feSpace[" << i << "]";
FiniteElemSpace *newFeSpace =
FiniteElemSpace::provideFeSpace(admin, Lagrange::getLagrange(dim, degree),
componentMeshes[i], s.str());
feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)] = newFeSpace;
feSpaces.push_back(newFeSpace);
}
componentSpaces[i] = feSpaceMap[std::pair<Mesh*, int>(componentMeshes[i], degree)];
}
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++)
traverseInfo.getMatrix(i, j).setFeSpace(componentSpaces[i], componentSpaces[j]);
traverseInfo.getVector(i).setFeSpace(componentSpaces[i]);
}
// create dof admin for vertex dofs if neccessary
for (int i = 0; i < static_cast<int>(meshes.size()); i++) {
if (meshes[i]->getNumberOfDOFs(VERTEX) == 0) {
DimVec<int> ln_dof(meshes[i]->getDim(), DEFAULT_VALUE, 0);
ln_dof[VERTEX] = 1;
meshes[i]->createDOFAdmin("vertex dofs", ln_dof);
}
}
}
void ProblemVec::createMatricesAndVectors()
{
FUNCNAME("ProblemVec::createMatricesAndVectors()");
// === create vectors and system matrix ===
systemMatrix = new Matrix<DOFMatrix*>(nComponents, nComponents);
systemMatrix->set(NULL);
rhs = new SystemVector("rhs", componentSpaces, nComponents);
solution = new SystemVector("solution", componentSpaces, nComponents);
for (int i = 0; i < nComponents; i++) {
(*systemMatrix)[i][i] = new DOFMatrix(componentSpaces[i],
componentSpaces[i], "A_ii");
(*systemMatrix)[i][i]->setCoupleMatrix(false);
std::string numberedName = "rhs[" + boost::lexical_cast<std::string>(i) + "]";
rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], numberedName));
numberedName = "solution[" + boost::lexical_cast<std::string>(i) + "]";
solution->setDOFVector(i, new DOFVector<double>(componentSpaces[i],
numberedName));
solution->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL);
solution->getDOFVector(i)->set(0.0);
}
}
void ProblemVec::createSolver()
{
FUNCNAME("ProblemVec::createSolver()");
// === create solver ===
std::string solverType("0");
GET_PARAMETER(0, name + "->solver", &solverType);
OEMSolverCreator *solverCreator =
dynamic_cast<OEMSolverCreator*>(CreatorMap<OEMSolver>::getCreator(solverType));
TEST_EXIT(solverCreator)("no solver type\n");
solverCreator->setName(name + "->solver");
solver = solverCreator->create();
solver->initParameters();
}
void ProblemVec::createEstimator()
{
FUNCNAME("ProblemVec::createEstimator()");
// create and set leaf data prototype
for (unsigned int i = 0; i < meshes.size(); i++)
meshes[i]->setElementDataPrototype
(new LeafDataEstimatableVec(new LeafDataCoarsenableVec));
for (int i = 0; i < nComponents; i++) {
TEST_EXIT(estimator[i] == NULL)("estimator already created\n");
std::string estName =
name + "->estimator[" + boost::lexical_cast<std::string>(i) + "]";
// === create estimator ===
std::string estimatorType("0");
GET_PARAMETER(0, estName, &estimatorType);
EstimatorCreator *estimatorCreator =
dynamic_cast<EstimatorCreator*>(CreatorMap<Estimator>::getCreator(estimatorType));
if (estimatorCreator) {
estimatorCreator->setName(estName);
estimatorCreator->setRow(i);
if (estimatorType == "recovery")
dynamic_cast<RecoveryEstimator::Creator*>(estimatorCreator)->
setSolution(solution->getDOFVector(i));
estimator[i] = estimatorCreator->create();
}
if (estimator[i])
for (int j = 0; j < nComponents; j++)
estimator[i]->addSystem((*systemMatrix)[i][j],
solution->getDOFVector(j),
rhs->getDOFVector(j));
}
}
void ProblemVec::createMarker()
{
FUNCNAME("ProblemVec::createMarker()");
int nMarkersCreated = 0;
for (int i = 0; i < nComponents; i++) {
marker[i] = Marker::createMarker
(name + "->marker[" + boost::lexical_cast<std::string>(i) + "]", i);
if (marker[i]) {
nMarkersCreated++;
// If there is more than one marker, and all components are defined
// on the same mesh, the maximum marking has to be enabled.
if (nMarkersCreated > 1 && nMeshes == 1)
marker[i]->setMaximumMarking(true);
}
}
}
void ProblemVec::createFileWriter()
{
FUNCNAME("ProblemVec::createFileWriter()");
// Create one filewriter for all components of the problem
std::string numberedName = name + "->output";
std::string filename = "";
GET_PARAMETER(0, numberedName + "->filename", &filename);
if (filename != "") {
std::vector< DOFVector<double>* > solutionList(nComponents);
for (int i = 0; i < nComponents; i++) {
TEST_EXIT(componentMeshes[0] == componentMeshes[i])
("All Meshes have to be equal to write a vector file.\n");
solutionList[i] = solution->getDOFVector(i);
}
fileWriters.push_back(new FileWriter(numberedName,
componentMeshes[0],
solutionList));
}
// Create own filewriters for each components of the problem
for (int i = 0; i < nComponents; i++) {
numberedName = name + "->output[" + boost::lexical_cast<std::string>(i) + "]";
filename = "";
GET_PARAMETER(0, numberedName + "->filename", &filename);
if (filename != "")
fileWriters.push_back(new FileWriter(numberedName,
componentMeshes[i],
solution->getDOFVector(i)));
}
// Check for serializer
int writeSerialization = 0;
GET_PARAMETER(0, name + "->write serialization", "%d", &writeSerialization);
if (writeSerialization) {
MSG("Use are using the obsolete parameter: %s->write serialization\n", name.c_str());
MSG("Please use instead the following parameter: %s->output->write serialization\n", name.c_str());
ERROR_EXIT("Usage of an obsolete parameter (see message above)!\n");
}
GET_PARAMETER(0, name + "->output->write serialization", "%d", &writeSerialization);
// Serialization is not allowed to be done by the problem, if its part of a parallel
// problem definition. Than, the parallel problem object must be serialized.
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
if (writeSerialization)
fileWriters.push_back(new Serializer<ProblemVec>(this));
#endif
}
void ProblemVec::doOtherStuff()
{
}
void ProblemVec::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
{
FUNCNAME("Problem::solve()");
if (!solver) {
WARNING("no solver\n");
return;
}
#ifdef _OPENMP
double wtime = omp_get_wtime();
#endif
clock_t first = clock();
solver->solveSystem(solverMatrix, *solution, *rhs);
#ifdef _OPENMP
INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
TIME_USED(first, clock()), omp_get_wtime() - wtime);
#else
INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
TIME_USED(first, clock()));
#endif
adaptInfo->setSolverIterations(solver->getIterations());
adaptInfo->setMaxSolverIterations(solver->getMaxIterations());
adaptInfo->setSolverTolerance(solver->getTolerance());
adaptInfo->setSolverResidual(solver->getResidual());
}
void ProblemVec::estimate(AdaptInfo *adaptInfo)
{
FUNCNAME("ProblemVec::estimate()");
clock_t first = clock();
#ifdef _OPENMP
double wtime = omp_get_wtime();
#endif
if (computeExactError) {
computeError(adaptInfo);
} else {
for (int i = 0; i < nComponents; i++) {
Estimator *scalEstimator = estimator[i];
if (scalEstimator) {
traverseInfo.updateStatus();
scalEstimator->setTraverseInfo(traverseInfo);
scalEstimator->estimate(adaptInfo->getTimestep());
adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);
adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
} else {
WARNING("No estimator for component %d\n" , i);
}
}
}
#ifdef _OPENMP
INFO(info, 8)("estimation of the error needed %.5f seconds system time / %.5f seconds wallclock time\n",
TIME_USED(first, clock()), omp_get_wtime() - wtime);
#else
INFO(info, 8)("estimation of the error needed %.5f seconds\n",
TIME_USED(first, clock()));
#endif
}
Flag ProblemVec::markElements(AdaptInfo *adaptInfo)
{
FUNCNAME("ProblemVec::markElements()");
TEST_EXIT_DBG(static_cast<unsigned int>(nComponents) == marker.size())
("Wrong number of markers!\n");
Flag markFlag = 0;
for (int i = 0; i < nComponents; i++) {
if (marker[i])
markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]);
else
WARNING("No marker for component %d\n", i);
}
return markFlag;
}
Flag ProblemVec::refineMesh(AdaptInfo *adaptInfo)
{
FUNCNAME("ProblemVec::refineMesh()");
int nMeshes = static_cast<int>(meshes.size());
Flag refineFlag = 0;
for (int i = 0; i < nMeshes; i++)
if (adaptInfo->isRefinementAllowed(i))
refineFlag |= refinementManager->refineMesh(meshes[i]);
return refineFlag;
}
Flag ProblemVec::coarsenMesh(AdaptInfo *adaptInfo)
{
FUNCNAME("ProblemVec::coarsenMesh()");
int nMeshes = static_cast<int>(meshes.size());
Flag coarsenFlag = 0;
for (int i = 0; i < nMeshes; i++)
if (adaptInfo->isCoarseningAllowed(i))
coarsenFlag |= coarseningManager->coarsenMesh(meshes[i]);
return coarsenFlag;
}
Flag ProblemVec::oneIteration(AdaptInfo *adaptInfo, Flag toDo)
{
FUNCNAME("ProblemVec::oneIteration()");
for (int i = 0; i < nComponents; i++)
if (adaptInfo->spaceToleranceReached(i))
adaptInfo->allowRefinement(false, i);
else
adaptInfo->allowRefinement(true, i);
return StandardProblemIteration::oneIteration(adaptInfo, toDo);
}
void ProblemVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
bool asmMatrix, bool asmVector)
{
FUNCNAME("ProblemVec::buildAfterCoarsen()");
if (dualMeshTraverseRequired()) {
dualAssemble(adaptInfo, flag);
return;
}
// printOpenmpTraverseInfo(this, true);
// std::cout << "ElInfo = " << ElInfo::subElemMatrices.size() << std::endl;
for (unsigned int i = 0; i < meshes.size(); i++)
meshes[i]->dofCompress();
clock_t first = clock();
#ifdef _OPENMP
double wtime = omp_get_wtime();
#endif
Flag assembleFlag =
flag |
(*systemMatrix)[0][0]->getAssembleFlag() |
rhs->getDOFVector(0)->getAssembleFlag() |
Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS |
Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_NEIGH;
if (useGetBound)
assembleFlag |= Mesh::FILL_BOUND;
traverseInfo.updateStatus();
// Used to calculate the overall number of non zero entries.
int nnz = 0;
for (int i = 0; i < nComponents; i++) {
MSG("%d DOFs for %s\n",
componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[i]->getName().c_str());
rhs->getDOFVector(i)->set(0.0);
for (int j = 0; j < nComponents; j++) {
if (writeAsmInfo)
std::cout << "-------" << i << " " << j << "----------------" << std::endl;
// Only if this variable is true, the current matrix will be assembled.
bool assembleMatrix = true;
// The DOFMatrix which should be assembled (or not, if assembleMatrix
// will be set to false).
DOFMatrix *matrix = (asmMatrix ? (*systemMatrix)[i][j] : NULL);
if (writeAsmInfo && matrix) {
for (std::vector<Operator*>::iterator it = matrix->getOperatorsBegin();
it != matrix->getOperatorsEnd(); ++it) {
Assembler *assembler = (*it)->getAssembler();
if (assembler->getZeroOrderAssembler())
std::cout << "ZOA: " << assembler->getZeroOrderAssembler()->getName() << std::endl;
if (assembler->getFirstOrderAssembler(GRD_PSI))
std::cout << "FOA GRD_PSI: " << assembler->getFirstOrderAssembler(GRD_PSI)->getName() << std::endl;
if (assembler->getFirstOrderAssembler(GRD_PHI))
std::cout << "FOA GRD_PHI: " << assembler->getFirstOrderAssembler(GRD_PHI)->getName() << std::endl;
if (assembler->getSecondOrderAssembler())
std::cout << "SOA: " << assembler->getSecondOrderAssembler()->getName() << std::endl;
}
}
if (matrix)
matrix->calculateNnz();
// If the matrix was assembled before and it is marked to be assembled
// only once, it will not be assembled.
if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
assembleMatrix = false;
} else if (matrix) {
matrix->getBaseMatrix().
change_dim(componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[j]->getAdmin()->getUsedSize());
set_to_zero(matrix->getBaseMatrix());
}
// If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
if (!matrix || !assembleMatrix)
assembleMatrix = false;
// If the matrix should not be assembled, the rhs vector has to be considered.
// This will be only done, if i == j. So, if both is not true, we can jump
// to the next matrix.
if (!assembleMatrix && i != j) {
if (matrix)
nnz += matrix->getBaseMatrix().nnz();
continue;
}
if (assembleMatrix && matrix->getBoundaryManager())
matrix->getBoundaryManager()->initMatrix(matrix);
// The simplest case: either the right hand side has no operaters, no aux
// fe spaces, or all aux fe spaces are equal to the row and col fe space.
assembleOnOneMesh(componentSpaces[i],
assembleFlag,
assembleMatrix ? matrix : NULL,
((i == j) && asmVector) ? rhs->getDOFVector(i) : NULL);
assembledMatrix[i][j] = true;
if (assembleMatrix)
matrix->finishInsertion();
if (assembleMatrix && matrix->getBoundaryManager())
matrix->getBoundaryManager()->exitMatrix(matrix);
if (matrix)
nnz += matrix->getBaseMatrix().nnz();
}
// And now assemble boundary conditions on the vectors
assembleBoundaryConditions(rhs->getDOFVector(i),
solution->getDOFVector(i),
componentMeshes[i],
assembleFlag);
}
if (asmMatrix) {
solverMatrix.setMatrix(*systemMatrix);
createPrecon();
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
}
#ifdef _OPENMP
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
TIME_USED(first, clock()), omp_get_wtime() - wtime);
#else
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
TIME_USED(first, clock()));
#endif
}
bool ProblemVec::dualMeshTraverseRequired()
{
FUNCNAME("ProblemVec::dualMeshTraverseRequired()");
TEST_EXIT(meshes.size() <= 2)("More than two mneshes are not supported yet!\n");
return (meshes.size() == 2);
}
void ProblemVec::dualAssemble(AdaptInfo *adaptInfo, Flag flag)
{
FUNCNAME("ProblemVec::dualAssemble()");
for (unsigned int i = 0; i < meshes.size(); i++)
meshes[i]->dofCompress();
clock_t first = clock();
#ifdef _OPENMP
double wtime = omp_get_wtime();
#endif
Flag assembleFlag =
flag |
(*systemMatrix)[0][0]->getAssembleFlag() |
rhs->getDOFVector(0)->getAssembleFlag() |
Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS |
Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_NEIGH;
if (useGetBound)
assembleFlag |= Mesh::FILL_BOUND;
traverseInfo.updateStatus();
// Used to calculate the overall number of non zero entries.
int nnz = 0;
for (int i = 0; i < nComponents; i++) {
MSG("%d DOFs for %s\n",
componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[i]->getName().c_str());
rhs->getDOFVector(i)->set(0.0);
for (int j = 0; j < nComponents; j++) {
if (writeAsmInfo)
std::cout << "-------" << i << " " << j << "----------------" << std::endl;
// Only if this variable is true, the current matrix will be assembled.
bool assembleMatrix = true;
// The DOFMatrix which should be assembled (or not, if assembleMatrix
// will be set to false).
DOFMatrix *matrix = (*systemMatrix)[i][j];
if (writeAsmInfo && matrix) {
for (std::vector<Operator*>::iterator it = matrix->getOperatorsBegin();
it != matrix->getOperatorsEnd(); ++it) {
Assembler *assembler = (*it)->getAssembler();
if (assembler->getZeroOrderAssembler())
std::cout << "ZOA: " << assembler->getZeroOrderAssembler()->getName() << std::endl;
if (assembler->getFirstOrderAssembler(GRD_PSI))
std::cout << "FOA GRD_PSI: " << assembler->getFirstOrderAssembler(GRD_PSI)->getName() << std::endl;
if (assembler->getFirstOrderAssembler(GRD_PHI))
std::cout << "FOA GRD_PHI: " << assembler->getFirstOrderAssembler(GRD_PHI)->getName() << std::endl;
if (assembler->getSecondOrderAssembler())
std::cout << "SOA: " << assembler->getSecondOrderAssembler()->getName() << std::endl;
}
}
if (matrix)
matrix->calculateNnz();
// If the matrix was assembled before and it is marked to be assembled
// only once, it will not be assembled.
if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
assembleMatrix = false;
} else if (matrix) {
matrix->getBaseMatrix().
change_dim(componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[j]->getAdmin()->getUsedSize());
set_to_zero(matrix->getBaseMatrix());
matrix->startInsertion(matrix->getNnz());
if (matrix->getBoundaryManager())
matrix->getBoundaryManager()->initMatrix(matrix);
}
// If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
if (!matrix || !assembleMatrix)
assembleMatrix = false;
// If the matrix should not be assembled, the rhs vector has to be considered.
// This will be only done, if i == j. So, if both is not true, we can jump
// to the next matrix.
if (!assembleMatrix && i != j)
if (matrix)
nnz += matrix->getBaseMatrix().nnz();
}
}
TEST_EXIT(meshes.size() == 2)("There is something wrong!\n");
const BasisFunction *basisFcts = componentSpaces[0]->getBasisFcts();
BoundaryType *bound =
useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
DualTraverse dualTraverse;
DualElInfo dualElInfo;
int oldElIndex0 = -1;
int oldElIndex1 = -1;
dualTraverse.setFillSubElemMat(true, basisFcts);
bool cont = dualTraverse.traverseFirst(meshes[0], meshes[1], -1, -1,
assembleFlag, assembleFlag, dualElInfo);
while (cont) {
bool newEl0 = (dualElInfo.rowElInfo->getElement()->getIndex() != oldElIndex0);
bool newEl1 = (dualElInfo.colElInfo->getElement()->getIndex() != oldElIndex1);
oldElIndex0 = dualElInfo.rowElInfo->getElement()->getIndex();
oldElIndex1 = dualElInfo.colElInfo->getElement()->getIndex();
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
DOFMatrix *matrix = (*systemMatrix)[i][j];
if (!matrix)
continue;
if (traverseInfo.eqSpaces(i, j)) {
ElInfo *elInfo = NULL;
if (componentMeshes[i] == meshes[0] && newEl0)
elInfo = dualElInfo.rowElInfo;
if (componentMeshes[i] == meshes[1] && newEl1)
elInfo = dualElInfo.colElInfo;
if (elInfo != NULL) {
if (useGetBound)
basisFcts->getBound(elInfo, bound);
if (matrix)
matrix->assemble(1.0, elInfo, bound);
if (i == j)
rhs->getDOFVector(i)->assemble(1.0, elInfo, bound);
}
ElInfo *mainElInfo, *auxElInfo;
if (traverseInfo.getRowFeSpace(i)->getMesh() == meshes[0]) {
mainElInfo = dualElInfo.rowElInfo;
auxElInfo = dualElInfo.colElInfo;
} else {
mainElInfo = dualElInfo.colElInfo;
auxElInfo = dualElInfo.rowElInfo;
}
if (useGetBound && mainElInfo != elInfo)
basisFcts->getBound(mainElInfo, bound);
if (traverseInfo.difAuxSpace(i) && i == j)
rhs->getDOFVector(i)->assemble2(1.0, mainElInfo, auxElInfo,
dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);
if (traverseInfo.difAuxSpace(i, j) && matrix)
matrix->assemble2(1.0, mainElInfo, auxElInfo,
dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);
if (matrix && matrix->getBoundaryManager())
matrix->getBoundaryManager()->fillBoundaryConditions(mainElInfo, matrix);
} else {
TEST_EXIT_DBG(traverseInfo.getStatus(i, j) != SingleComponentInfo::DIF_SPACES_WITH_DIF_AUX)
("Not yet supported!\n");
ElInfo *rowElInfo, *colElInfo;
if (componentMeshes[i] == meshes[0]) {
rowElInfo = dualElInfo.rowElInfo;
colElInfo = dualElInfo.colElInfo;
} else {
rowElInfo = dualElInfo.colElInfo;
colElInfo = dualElInfo.rowElInfo;
}
if (useGetBound)
basisFcts->getBound(rowElInfo, bound);
if (matrix) {
matrix->assemble(1.0, rowElInfo, colElInfo,
dualElInfo.smallElInfo, dualElInfo.largeElInfo, bound);
if (matrix->getBoundaryManager())
matrix->getBoundaryManager()->fillBoundaryConditions(rowElInfo, matrix);
}
if (i == j)
rhs->getDOFVector(i)->assemble(1.0, rowElInfo, bound);
}
}
}
cont = dualTraverse.traverseNext(dualElInfo);
}
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
DOFMatrix *matrix = (*systemMatrix)[i][j];
if (matrix)
matrix->removeRowsWithDBC(matrix->getApplyDBCs());
if (matrix)
matrix->finishInsertion();
if (matrix && matrix->getBoundaryManager())
matrix->getBoundaryManager()->exitMatrix(matrix);
if (matrix)
nnz += matrix->getBaseMatrix().nnz();
}
// And now assemble boundary conditions on the vectors
assembleBoundaryConditions(rhs->getDOFVector(i),
solution->getDOFVector(i),
componentMeshes[i],
assembleFlag);
}
solverMatrix.setMatrix(*systemMatrix);
createPrecon();
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
#ifdef _OPENMP
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds system time / %.5f seconds wallclock time\n",
TIME_USED(first, clock()), omp_get_wtime() - wtime);
#else
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
TIME_USED(first, clock()));
#endif
}
void ProblemVec::createPrecon()
{
std::string preconType("no");
GET_PARAMETER(0, name + "->solver->left precon", &preconType);
CreatorInterface<ITL_BasePreconditioner> *preconCreator =
CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
solver->setLeftPrecon(preconCreator->create(solverMatrix.getMatrix()));
preconType= "no";
GET_PARAMETER(0, name + "->solver->right precon", &preconType);
preconCreator = CreatorMap<ITL_BasePreconditioner>::getCreator(preconType);
solver->setRightPrecon(preconCreator->create(solverMatrix.getMatrix()));
}
void ProblemVec::writeFiles(AdaptInfo *adaptInfo, bool force)
{
FUNCNAME("ProblemVec::writeFiles()");
clock_t first = clock();
#ifdef _OPENMP
double wtime = omp_get_wtime();
#endif
int i;
int size = static_cast<int>(fileWriters.size());
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 1)
#endif
for (i = 0; i < size; i++) {
fileWriters[i]->writeFiles(adaptInfo, force);
}
#ifdef _OPENMP
INFO(info, 8)("writeFiles needed %.5f seconds system time / %.5f seconds wallclock time\n",
TIME_USED(first, clock()),
omp_get_wtime() - wtime);
#else
INFO(info, 8)("writeFiles needed %.5f seconds\n",
TIME_USED(first, clock()));
#endif
}
void ProblemVec::writeFiles(AdaptInfo &adaptInfo, bool force)
{
writeFiles(&adaptInfo, force);
}
void ProblemVec::interpolInitialSolution(std::vector<AbstractFunction<double, WorldVector<double> >*> *fct)
{
FUNCNAME("ProblemVec::interpolInitialSolution()");
solution->interpol(fct);
}
void ProblemVec::addMatrixOperator(Operator *op, int i, int j,
double *factor, double *estFactor)
{
FUNCNAME("ProblemVec::addMatrixOperator()");
TEST_EXIT(!boundaryConditionSet)
("Do not add operators after boundary conditions were set!\n");
if (!(*systemMatrix)[i][j]) {
TEST_EXIT(i != j)("should have been created already\n");
(*systemMatrix)[i][j] = new DOFMatrix(componentSpaces[i], componentSpaces[j], "");
(*systemMatrix)[i][j]->setCoupleMatrix(true);
(*systemMatrix)[i][j]->getBoundaryManager()->
setBoundaryConditionMap((*systemMatrix)[i][i]->getBoundaryManager()->
getBoundaryConditionMap());
if (estimator[i])
estimator[i]->setNewMatrix(j, (*systemMatrix)[i][j]);
}
(*systemMatrix)[i][j]->addOperator(op, factor, estFactor);
traverseInfo.getMatrix(i, j).setAuxFeSpaces(op->getAuxFeSpaces());
for (std::set<const FiniteElemSpace*>::iterator it = op->getAuxFeSpaces().begin();
it != op->getAuxFeSpaces().end(); ++it) {
if ((*it)->getMesh() != componentSpaces[i]->getMesh() ||
(*it)->getMesh() != componentSpaces[j]->getMesh()) {
op->setNeedDualTraverse(true);
break;
}
}
OperatorPos opPos = {i, j, factor, estFactor, Operator::MATRIX_OPERATOR};
operators[op].push_back(opPos);
opFlags[op].setFlag(Operator::MATRIX_OPERATOR);
}
void ProblemVec::addMatrixOperator(Operator &op, int i, int j,
double *factor, double *estFactor)
{
addMatrixOperator(&op, i, j, factor, estFactor);
}
void ProblemVec::addVectorOperator(Operator *op, int i,
double *factor, double *estFactor)
{
FUNCNAME("ProblemVec::addVectorOperator()");
TEST_EXIT(!boundaryConditionSet)
("Do not add operators after boundary conditions were set!\n");
rhs->getDOFVector(i)->addOperator(op, factor, estFactor);
traverseInfo.getVector(i).setAuxFeSpaces(op->getAuxFeSpaces());
for (std::set<const FiniteElemSpace*>::iterator it = op->getAuxFeSpaces().begin();
it != op->getAuxFeSpaces().end(); ++it) {
if ((*it)->getMesh() != componentSpaces[i]->getMesh()) {
op->setNeedDualTraverse(true);
break;
}
}
OperatorPos opPos = {i, -1, factor, estFactor, Operator::VECTOR_OPERATOR};
operators[op].push_back(opPos);
opFlags[op].setFlag(Operator::VECTOR_OPERATOR);
}
void ProblemVec::addVectorOperator(Operator &op, int i,
double *factor, double *estFactor)
{
addVectorOperator(&op, i, factor, estFactor);
}
void ProblemVec::addDirichletBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> >* b)
{
FUNCNAME("ProblemVec::addDirichletBC()");
TEST_EXIT(row >= 0 && row < nComponents)("Wrong row number: %d\n", row);
TEST_EXIT(col >= 0 && col < nComponents)("Wrong col number: %d\n", col);
boundaryConditionSet = true;
DirichletBC *dirichletApply =
new DirichletBC(type, b, componentSpaces[row], componentSpaces[col], true);
DirichletBC *dirichletNotApply =
new DirichletBC(type, b, componentSpaces[row], componentSpaces[col], false);
for (int i = 0; i < nComponents; i++)
if (systemMatrix && (*systemMatrix)[row][i])
if (i == col)
(*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletApply);
else
(*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletNotApply);
if (rhs)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
if (solution)
solution->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
}
void ProblemVec::addDirichletBC(BoundaryType type, int row, int col,
DOFVector<double> *vec)
{
FUNCNAME("ProblemVec::addDirichletBC()");
TEST_EXIT(row >= 0 && row < nComponents)("Wrong row number: %d\n", row);
TEST_EXIT(col >= 0 && col < nComponents)("Wrong col number: %d\n", col);
boundaryConditionSet = true;
DirichletBC *dirichletApply = new DirichletBC(type, vec, true);
DirichletBC *dirichletNotApply = new DirichletBC(type, vec, false);
for (int i = 0; i < nComponents; i++)
if (systemMatrix && (*systemMatrix)[row][i])
if (i == col)
(*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletApply);
else
(*systemMatrix)[row][i]->getBoundaryManager()->addBoundaryCondition(dirichletNotApply);
if (rhs)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
if (solution)
solution->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(dirichletApply);
}
void ProblemVec::addNeumannBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *n)
{
FUNCNAME("ProblemVec::addNeumannBC()");
boundaryConditionSet = true;
NeumannBC *neumann =
new NeumannBC(type, n, componentSpaces[row], componentSpaces[col]);
if (rhs)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(neumann);
}
void ProblemVec::addRobinBC(BoundaryType type, int row, int col,
AbstractFunction<double, WorldVector<double> > *n,
AbstractFunction<double, WorldVector<double> > *r)
{
FUNCNAME("ProblemVec::addRobinBC()");
boundaryConditionSet = true;
RobinBC *robin =
new RobinBC(type, n, r, componentSpaces[row], componentSpaces[col]);
if (systemMatrix && (*systemMatrix)[row][col])
(*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(robin);
if (rhs)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(robin);
}
void ProblemVec::addPeriodicBC(BoundaryType type, int row, int col)
{
FUNCNAME("ProblemVec::addPeriodicBC()");
boundaryConditionSet = true;
FiniteElemSpace *feSpace = componentSpaces[row];
PeriodicBC *periodic = new PeriodicBC(type, feSpace);
if (systemMatrix && (*systemMatrix)[row][col])
(*systemMatrix)[row][col]->getBoundaryManager()->addBoundaryCondition(periodic);
if (rhs && row == col)
rhs->getDOFVector(row)->getBoundaryManager()->addBoundaryCondition(periodic);
}
void ProblemVec::assembleOnOneMesh(FiniteElemSpace *feSpace,
Flag assembleFlag,
DOFMatrix *matrix, DOFVector<double> *vector)
{
Mesh *mesh = feSpace->getMesh();
const BasisFunction *basisFcts = feSpace->getBasisFcts();
#ifdef _OPENMP
TraverseParallelStack stack(0, 1);
#else
TraverseStack stack;
#endif
// == Initialize matrix and vector. If we have to assemble in parallel, ===
// == temporary thread owned matrix and vector must be created. ===
#ifdef _OPENMP
#pragma omp parallel
#endif
{
BoundaryType *bound =
useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
// Create for every thread its private matrix and vector, on that
// the thread will assemble its part of the mesh.
DOFMatrix *tmpMatrix = NULL;
DOFVector<double> *tmpVector = NULL;
#ifdef _OPENMP
if (matrix) {
tmpMatrix = new DOFMatrix(matrix->getRowFeSpace(), matrix->getColFeSpace(), "tmp");
// Copy the global matrix to the private matrix, because we need the
// operators defined on the global matrix in the private one. Only the
// values have to be set to zero.
*tmpMatrix = *matrix;
tmpMatrix->clear();
tmpMatrix->getBaseMatrix().change_dim(matrix->getRowFeSpace()->getAdmin()->getUsedSize(),
matrix->getColFeSpace()->getAdmin()->getUsedSize());
tmpMatrix->startInsertion(10);
}
if (vector) {
tmpVector = new DOFVector<double>(vector->getFeSpace(), "tmp");
// Copy the global vector to the private vector, because we need the
// operatirs defined on the global vector in the private one. But set
// the values to zero of the private vector after copying.
*tmpVector = *vector;
tmpVector->set(0.0);
}
#else
if (matrix) {
tmpMatrix = matrix;
tmpMatrix->startInsertion(matrix->getNnz());
}
if (vector) {
tmpVector = vector;
tmpVector->set(0.0);
}
#endif
// == Traverse mesh (either sequentially or in parallel) and assemble. ==
// Because we are using the parallel traverse stack, each thread will
// traverse only a part of the mesh.
ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
// After creating privat copies of the DOFMatrix and the DOFVector, all threads
// have to wait at this barrier. Especially for small problems this is required,
// because otherwise one thread may be finished with assembling, before another
// has made his private copy.
#ifdef _OPENMP
#pragma omp barrier
#endif
while (elInfo) {
if (useGetBound)
basisFcts->getBound(elInfo, bound);
if (matrix) {
tmpMatrix->assemble(1.0, elInfo, bound);
// Take the matrix boundary manager from the public matrix,
// but assemble the boundary conditions on the thread private matrix.
if (matrix->getBoundaryManager())
matrix->getBoundaryManager()->fillBoundaryConditions(elInfo, tmpMatrix);
}
if (vector)
tmpVector->assemble(1.0, elInfo, bound, NULL);
elInfo = stack.traverseNext(elInfo);
}
// == Finally, if we have assembled in parallel, we have to add the thread ==
// == private matrix and vector to the global one. ==
#ifdef _OPENMP
tmpMatrix->finishInsertion();
// After mesh traverse, all thread have to added their private matrices and
// vectors to the global public matrix and public vector. Therefore, this is
// a critical section, which is allowed to be executed by on thread only at
// the same time.
if (matrix) {
#pragma omp critical
{
matrix->getBaseMatrix() += tmpMatrix->getBaseMatrix();
}
}
#pragma omp barrier
#pragma omp master
{
if (matrix)
matrix->startInsertion();
}
#pragma omp barrier
if (matrix) {
// Remove rows corresponding to DOFs on a Dirichlet boundary.
#pragma omp critical
matrix->removeRowsWithDBC(tmpMatrix->getApplyDBCs());
delete tmpMatrix;
}
if (vector) {
#pragma omp critical
*vector += *tmpVector;
delete tmpVector;
}
#else
if (matrix)
matrix->removeRowsWithDBC(matrix->getApplyDBCs());
#endif
if (useGetBound)
delete [] bound;
} // pragma omp parallel
}
void ProblemVec::assembleBoundaryConditions(DOFVector<double> *rhs,
DOFVector<double> *solution,
Mesh *mesh,
Flag assembleFlag)
{
FUNCNAME("ProblemVec::assembleBoundaryConditions()");
// === Initialization of vectors ===
if (rhs->getBoundaryManager())
rhs->getBoundaryManager()->initVector(rhs);
if (solution->getBoundaryManager())
solution->getBoundaryManager()->initVector(solution);
#ifdef _OPENMP
// === Parallel boundary assemblage ===
TraverseParallelStack stack;
#pragma omp parallel
{
// Each thread assembles on its own dof-vectors.
DOFVector<double> *tmpRhsVec = new DOFVector<double>(rhs->getFeSpace(), "tmpRhs");
DOFVector<double> *tmpSolVec = new DOFVector<double>(solution->getFeSpace(), "tmpSol");
tmpRhsVec->set(0.0);
tmpSolVec->set(0.0);
// (Parallel) traverse of mesh.
ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
while (elInfo) {
if (rhs->getBoundaryManager())
rhs->getBoundaryManager()-> fillBoundaryConditions(elInfo, tmpRhsVec);
if (solution->getBoundaryManager())
solution->getBoundaryManager()->fillBoundaryConditions(elInfo, tmpSolVec);
elInfo = stack.traverseNext(elInfo);
}
// After (parallel) mesh traverse, the result is applied to the final
// vectors. This section is not allowed to be executed by more than one
// thread at the same time.
#pragma omp critical
{
DOFVector<double>::Iterator rhsIt(rhs, USED_DOFS);
DOFVector<double>::Iterator solIt(solution, USED_DOFS);
DOFVector<double>::Iterator tmpRhsIt(tmpRhsVec, USED_DOFS);
DOFVector<double>::Iterator tmpSolIt(tmpSolVec, USED_DOFS);
for (rhsIt.reset(), solIt.reset(), tmpRhsIt.reset(), tmpSolIt.reset();
!rhsIt.end();
++rhsIt, ++solIt, ++tmpRhsIt, ++tmpSolIt) {
*rhsIt += *tmpRhsIt;
*solIt += *tmpSolIt;
}
} // pragma omp critical
delete tmpRhsVec;
delete tmpSolVec;
} // pragma omp parallel
#else
// === Sequentiell boundary assemblage ===
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
while (elInfo) {
if (rhs->getBoundaryManager())
rhs->getBoundaryManager()->fillBoundaryConditions(elInfo, rhs);
if (solution->getBoundaryManager())
solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution);
elInfo = stack.traverseNext(elInfo);
}
#endif
// === Finalize vectors ===
if (rhs->getBoundaryManager())
rhs->getBoundaryManager()->exitVector(rhs);
if (solution->getBoundaryManager())
solution->getBoundaryManager()->exitVector(solution);
}
void ProblemVec::writeResidualMesh(int comp, AdaptInfo *adaptInfo, std::string name)
{
FUNCNAME("ProblemVec::writeResidualMesh()");
std::map<int, double> vec;
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(this->getMesh(comp), -1,
Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
while (elInfo) {
vec[elInfo->getElement()->getIndex()] = elInfo->getElement()->getEstimation(comp);
elInfo = stack.traverseNext(elInfo);
}
ElementFileWriter::writeFile(vec, this->getFeSpace(comp), name);
}
void ProblemVec::serialize(std::ostream &out)
{
FUNCNAME("ProblemVec::serialize()");
for (int i = 0; i < static_cast<int>(meshes.size()); i++)
meshes[i]->serialize(out);
solution->serialize(out);
}
void ProblemVec::deserialize(std::istream &in)
{
FUNCNAME("ProblemVec::deserialize()");
for (int i = 0; i < static_cast<int>(meshes.size()); i++)
meshes[i]->deserialize(in);
solution->deserialize(in);
}
void ProblemVec::computeError(AdaptInfo *adaptInfo)
{
FUNCNAME("ProblemVec::computeError()");
for (int i = 0; i < nComponents; i++) {
TEST_EXIT(exactSolutionFcts[i])("No solution function given!\n");
// Compute the difference between exact and computed solution
DOFVector<double> *tmp = new DOFVector<double>(componentSpaces[i], "tmp");
tmp->interpol(exactSolutionFcts[i]);
double solMax = tmp->absMax();
*tmp -= *(solution->getDOFVector(i));
MSG("L2 error = %.8e\n", tmp->L2Norm());
MSG("L-inf error = %.8e\n", tmp->absMax() / solMax);
adaptInfo->setEstSum(tmp->absMax() / solMax, i);
adaptInfo->setEstMax(tmp->absMax() / solMax, i);
// To set element estimates, compute a vector with the difference
// between exact and computed solution for each DOF.
DOFVector<double> *sol = new DOFVector<double>(componentSpaces[i], "tmp");
sol->interpol(exactSolutionFcts[i]);
DOFVector<double>::Iterator it1(sol, USED_DOFS);
DOFVector<double>::Iterator it2(tmp, USED_DOFS);
for (it1.reset(), it2.reset(); !it1.end(); ++it1, ++it2) {
if (abs(*it1) <= DBL_TOL || abs(*it2) <= DBL_TOL)
*it2 = 0.0;
else
*it2 = abs(*it2 / *it1);
}
// Compute estimate for every mesh element
std::vector<DegreeOfFreedom> locInd(componentSpaces[i]->getBasisFcts()->getNumber());
TraverseStack stack;
ElInfo *elInfo = stack.traverseFirst(componentMeshes[i], -1, Mesh::CALL_LEAF_EL);
while (elInfo) {
componentSpaces[i]->getBasisFcts()->getLocalIndices(elInfo->getElement(),
componentSpaces[i]->getAdmin(),
locInd);
double estimate = 0.0;
for (int j = 0; j < componentSpaces[i]->getBasisFcts()->getNumber(); j++)
estimate += (*tmp)[locInd[j]];
elInfo->getElement()->setEstimation(estimate, i);
elInfo->getElement()->setMark(0);
elInfo = stack.traverseNext(elInfo);
}
delete tmp;
delete sol;
}
}
}