Commit d30bf58a authored by jonathan.froehlich's avatar jonathan.froehlich
Browse files

[local-move-cardiac-transfers] Moved special transfer types to project

parent c57a5a8d
......@@ -39,7 +39,6 @@ set(fem_src
transfer/LinearTransfer.cpp
transfer/LagrangeTransfer.cpp
transfer/MixedLagrangeTransfer.cpp
transfer/MultiPartTransfer.cpp
transfer/Transfers.cpp
transfer/TransferWrapper.cpp
plot/Plotting.cpp
......
#ifndef ITRANSFER_HPP
#define ITRANSFER_HPP
#include "Config.hpp"
#include "Rowlist.hpp"
#include "RMatrix.hpp"
#include "Operator.hpp"
......
#include <exception>
#include <LagrangeDoF.hpp>
#include "MultiPartTransfer.hpp"
#include "LinearTransfer.hpp"
MultiPartTransfer::MultiPartTransfer(const Vector &coarse, const Vector &fine) : rowList(coarse,
fine),
discs(
(const LagrangeDiscretization &) coarse.GetDisc(),
(const MultiPartDiscretization &) fine.GetDisc()) {
constructRowList(coarse, fine);
}
void MultiPartTransfer::Prolongate(const Vector &coarse, Vector &fine) const {
if (&coarse.GetDisc() != &discs.first || &fine.GetDisc() != &discs.second) {
THROW("Discretizations of transfer and vectors do not match")
}
fine = 0.0;
for (int i = 0; i < rowList.IdSize(); ++i) {
for (int j = 0; j < fine.Dof(i); ++j) {
fine(i, j) = rowList.CalculateFineEntry(i, coarse, 0);
}
}
fine.ClearDirichletValues();
}
void MultiPartTransfer::ProlongateTransposed(Vector &coarse, const Vector &fine) const {
if (&coarse.GetDisc() != &discs.first || &fine.GetDisc() != &discs.second) {
THROW("Discretizations of transfer and vectors do not match")
}
coarse = 0.0;
for (int i = 0; i < rowList.IdSize(); ++i) {
for (const auto &cRow: rowList.CoarseRows(i)) {
for (int j = 0; j < fine.Dof(i); ++j) {
coarse(cRow.Id, 0) += cRow.Weight * fine(i, j) / fine.Dof(i);
}
}
}
coarse.ClearDirichletValues();
coarse.Collect();
}
void MultiPartTransfer::Restrict(Vector &coarse, const Vector &fine) const {
if (&coarse.GetDisc() != &discs.first || &fine.GetDisc() != &discs.second) {
THROW("Discretizations of transfer and vectors do not match")
}
coarse = 0.0;
for (int i = 0; i < rowList.IdSize(); ++i) {
for (const auto &cRow: rowList.CoarseRows(i)) {
for (int j = 0; j < fine.Dof(i); ++j) {
coarse(cRow.Id, 0) +=
cRow.Weight / rowList.CoarseWeight(cRow.Id) * fine(i, j) / fine.Dof(i);
}
}
}
coarse.ClearDirichletValues();
}
void MultiPartTransfer::Project(Vector &coarse, const Vector &fine) const {
if (&coarse.GetDisc() != &discs.first || &fine.GetDisc() != &discs.second) {
THROW("Discretizations of transfer and vectors do not match")
}
for (row r = coarse.rows(); r != coarse.rows_end(); ++r) {
coarse(r, 0) = 0.0;
row fineR = fine.find_row(r());
if (fineR == fine.rows_end()) continue;
for (int j = 0; j < fineR.n(); ++j) {
coarse(r, 0) += fine(fineR, j) / fineR.n();
}
}
coarse.ClearDirichletValues();
}
double checkRefinement(std::pair<int, int> degrees, std::pair<int, int> levels) {
int dlevel = (int) log2(degrees.first);
if (dlevel < log2(degrees.first)) {
THROW("degree " + std::to_string(degrees.first) + " transfer not implemented")
}
if (((int) log2(degrees.first)) + levels.second < dlevel + levels.first) {
THROW("Can't interpolate: Coarser mesh is finer than fine mesh.")
}
return pow(2, levels.second - (levels.first + dlevel)) * degrees.second;
}
void MultiPartTransfer::constructRowList(
const Vector &coarse, const Vector &fine) {
int dlevel = (int) log2(discs.first.Degree());
double lvl = checkRefinement({discs.first.Degree(), discs.second.Degree()},
{coarse.SpaceLevel(), fine.SpaceLevel()});
std::vector<int> coarseSD = coarse.GetMesh().IncludedSubdomains();
std::vector<int> fineSD = fine.GetMesh().IncludedSubdomains();
std::vector<int> intersectSD;
std::set_intersection(coarseSD.begin(), coarseSD.end(), fineSD.begin(), fineSD.end(),
std::back_inserter(intersectSD));
std::vector<Point> childrenVerteces{};
std::vector<Point> childVerteces{};
std::vector<Point> coarseNodalPoints{};
for (cell c = coarse.cells(); c != coarse.cells_end(); ++c) {
if (std::find(intersectSD.begin(), intersectSD.end(), c.Subdomain()) ==
intersectSD.end()) {
continue;
}
std::vector<std::vector<std::unique_ptr<Cell>>> recursiveCells(dlevel + 1);
recursiveCells[0].emplace_back(CreateCell(c.Type(), c.Subdomain(), c.AsVector(), false));
for (int l = 1; l <= dlevel; ++l) {
for (const auto &parent: recursiveCells[l - 1]) {
const vector<Rule> &R = parent->Refine(childrenVerteces);
for (int i = 0; i < R.size(); ++i) {
R[i](childrenVerteces, childVerteces);
recursiveCells[l].emplace_back(
CreateCell(R[i].type(), c.Subdomain(), childVerteces, false));
}
}
}
for (const auto &child: recursiveCells[dlevel]) {
LagrangeNodalPoints(*child, coarseNodalPoints, 1);
vector<int> coarseIds = mpp::nodalIds(coarse, coarseNodalPoints);
mpp::linearTransfer(child->ReferenceType(), coarseNodalPoints, coarseIds, lvl, fine, rowList);
}
}
}
RMatrix MultiPartTransfer::AsMatrix() const {
RMatrix transferMatrix((int) rowList.IdSize(), (int) rowList.WeightSize());
for (int i = 0; i < rowList.IdSize(); ++i) {
for (auto r: rowList.CoarseRows(i)) {
transferMatrix(i, r.Id) = r.Weight;
}
}
return transferMatrix;
}
#ifndef MULTIPARTTRANSFER_HPP
#define MULTIPARTTRANSFER_HPP
#include "MultiPartDiscretization.hpp"
#include "ITransfer.hpp"
class MultiPartTransfer {
std::pair<const LagrangeDiscretization &, const MultiPartDiscretization &> discs;
RowList rowList;
void constructRowList(const Vector &coarse, const Vector &fine);
public:
MultiPartTransfer(const Vector &coarse, const Vector &fine);
void Prolongate(const Vector &coarse, Vector &fine) const;
void ProlongateTransposed(Vector &coarse, const Vector &fine) const;
void Restrict(Vector &coarse, const Vector &fine) const;
void Project(Vector &coarse, const Vector &fine) const;
RMatrix AsMatrix() const;
};
#endif //MULTIPARTTRANSFER_HPP
#include "Transfers.hpp"
#include "LagrangeTransfer.hpp"
#include "MultiPartTransfer.hpp"
#include "TransferWrapper.hpp"
TRANSFERTYPE getTransferType(const IDiscretization &disc) {
if (typeid(disc) == typeid(LagrangeDiscretization))
return LAGRANGE;
if (typeid(disc) == typeid(MultiPartDiscretization))
return MULTIPART;
THROW("Transfertype not implemented for " + disc.DiscName())
}
std::unique_ptr<ITransfer>
getTransfer(TRANSFERTYPE type, const Vector &coarse, const Vector &fine) {
switch (type) {
case LAGRANGE:
return std::make_unique<ITransfer>(new LagrangeTransfer(coarse, fine));
case MULTIPART:
return std::make_unique<ITransfer>(new MultiPartTransfer(coarse, fine));
default:
THROW("No Transfer defined.")
}
}
std::unique_ptr<ITransfer> getTransfer(const Vector &coarse, const Vector &fine) {
TRANSFERTYPE coarseType = getTransferType(coarse.GetDisc());
TRANSFERTYPE fineType = getTransferType(fine.GetDisc());
const auto& coarseType = typeid(coarse.GetDisc());
const auto& fineType = typeid(fine.GetDisc());
if (coarseType == LAGRANGE) { //from LagrangeDiscretization to
if (fineType == LAGRANGE)
return std::make_unique<ITransfer>(new LagrangeTransfer(coarse, fine));
if (fineType == MULTIPART)
return std::make_unique<ITransfer>(new MultiPartTransfer(coarse, fine));
if (coarseType != fineType){
THROW("Discretizations of Vectors don't match")
}
if (fineType == typeid(LagrangeDiscretization)) {
return std::make_unique<ITransfer>(new LagrangeTransfer(coarse, fine));
}
return std::make_unique<ITransfer>(new TransferWrapper(coarse, fine));
}
......
......@@ -4,15 +4,6 @@
#include <Config.hpp>
#include "ITransfer.hpp"
//TODO: Move to cpp after deleting getTransfer(type, coarse, fine)
enum TRANSFERTYPE {
LAGRANGE,
MULTIPART
};
[[deprecated("Use getTransfer(coarse, fine) instead")]]
std::unique_ptr<ITransfer> getTransfer(TRANSFERTYPE type, const Vector &coarse, const Vector &fine);
std::unique_ptr<ITransfer> getTransfer(const Vector &coarse, const Vector &fine);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment