Skip to content

GitLab

Commit 3236301f authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* AMDiS changed to compile Andreas' vesicle code

parent dec7dac6
Hide whitespace changes
Inline Side-by-side
Showing with 590 additions and 45 deletions
AMDiS/src/Assembler.cc
View file @ 3236301f
......@@ -275,7 +275,8 @@ namespace AMDiS {
}
ElementMatrix *Assembler::initElementMatrix(ElementMatrix *elMat,
const ElInfo *elInfo)
const ElInfo *rowElInfo,
const ElInfo *colElInfo)
{
if (!elMat) {
elMat = NEW ElementMatrix(nRow, nCol);
......@@ -283,17 +284,23 @@ namespace AMDiS {
elMat->set(0.0);
rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
rowFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
rowFESpace->getAdmin(),
&(elMat->rowIndices));
if (rowFESpace == colFESpace) {
elMat->colIndices = elMat->rowIndices;
} else {
colFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
colFESpace->getAdmin(),
&(elMat->colIndices));
if (colElInfo) {
colFESpace->getBasisFcts()->getLocalIndicesVec(colElInfo->getElement(),
colFESpace->getAdmin(),
&(elMat->colIndices));
} else {
// If there is no colElInfo pointer, use rowElInfo the get the indices.
colFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
colFESpace->getAdmin(),
&(elMat->colIndices));
}
}
return elMat;
......@@ -309,7 +316,6 @@ namespace AMDiS {
elVec->set(0.0);
DOFAdmin *rowAdmin = rowFESpace->getAdmin();
Element *element = elInfo->getElement();
rowFESpace->getBasisFcts()->getLocalIndicesVec(element, rowAdmin, &(elVec->dofIndices));
......
AMDiS/src/Assembler.h
View file @ 3236301f
......@@ -70,7 +70,8 @@ namespace AMDiS {
virtual ~Assembler() {};
ElementMatrix *initElementMatrix(ElementMatrix *elMat,
const ElInfo *elInfo);
const ElInfo *rowElInfo,
const ElInfo *colElInfo = NULL);
ElementVector *initElementVector(ElementVector *elVec,
......
AMDiS/src/DOFMatrix.cc
View file @ 3236301f
......@@ -245,7 +245,6 @@ namespace AMDiS {
for (int i = 0; i < nRow; i++) { // for all rows of element matrix
row = elMat.rowIndices[i];
BoundaryCondition *condition =
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
......@@ -402,30 +401,28 @@ namespace AMDiS {
void DOFMatrix::freeDOFContent(int index)
{
int i, j, col=0, col2;
int col = 0;
if (0 < matrix[index].size()) {
// for all columns in this row
int size = static_cast<int>(matrix[index].size());
for (i=0; i<size; i++) {
for (int i = 0; i < size; i++) {
// if entry is used
if (entryUsed(index,i)) {
if (entryUsed(index, i)) {
// get column of this entry
col = matrix[index][i].col;
if (col != index) { // remove symmetric entry if exists
int colsize = static_cast<int>(matrix[col].size());
for (j=0; j< colsize; j++) {
col2 = matrix[col][j].col;
for (int j = 0; j < colsize; j++) {
int col2 = matrix[col][j].col;
if (col2 == index) {
matrix[col][j].col = DOFMatrix::UNUSED_ENTRY;
}
else if (col2 == DOFMatrix::NO_MORE_ENTRIES) {
} else if (col2 == DOFMatrix::NO_MORE_ENTRIES) {
break;
}
}
}
}
else if (col == DOFMatrix::NO_MORE_ENTRIES) {
} else if (col == DOFMatrix::NO_MORE_ENTRIES) {
break;
}
}
......@@ -473,7 +470,8 @@ namespace AMDiS {
}
Operator *operat = op ? op : operators[0];
operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, rowElInfo);
operat->getAssembler(omp_get_thread_num())->
initElementMatrix(elementMatrix, rowElInfo, colElInfo);
if (op) {
op->getElementMatrix(rowElInfo, colElInfo, elementMatrix);
......@@ -488,7 +486,7 @@ namespace AMDiS {
*factorIt ? **factorIt : 1.0);
}
}
addElementMatrix(factor, *elementMatrix, bound);
}
......
AMDiS/src/DOFVector.cc
View file @ 3236301f
......@@ -655,7 +655,9 @@ namespace AMDiS {
}
DOFVectorDOF::DOFVectorDOF() : DOFVector<DegreeOfFreedom>() {};
DOFVectorDOF::DOFVectorDOF() :
DOFVector<DegreeOfFreedom>()
{};
void DOFVectorDOF::freeDOFContent(DegreeOfFreedom dof) {
......
AMDiS/src/DOFVector.h
View file @ 3236301f
......@@ -316,8 +316,8 @@ namespace AMDiS {
*/
DOFVector()
: DOFVectorBase<T>(),
refineInter(false),
feSpace(NULL),
refineInter(false),
coarsenOperation(NO_OPERATION)
{};
......@@ -531,6 +531,11 @@ namespace AMDiS {
inline double l1norm() const {
return asum();
};
/** \brief
* Calculates doublewell of this DOFVector
*/
double DoubleWell(Quadrature* q = NULL) const;
/** \brief
* Calculates the sum of this DOFVector
......@@ -637,6 +642,11 @@ namespace AMDiS {
*/
static int H1Norm_fct(ElInfo* elinfo);
/** \brief
* Used by DoubleWell while mesh traversal
*/
static int DoubleWell_fct(ElInfo* elinfo);
protected:
/** \brief
* Name of this DOFVector
......@@ -859,6 +869,15 @@ namespace AMDiS {
const DOFVector<T>& v2,
DOFVector<T>& result);
template<typename T>
inline const DOFVector<T>& mod(const DOFVector<T>& v,
DOFVector<T>& result);
template<typename T>
inline const DOFVector<T>& Tanh(const DOFVector<T>& v,
DOFVector<T>& result);
template<typename T>
inline void set(DOFVector<T>& vec, T d)
{
......
AMDiS/src/DOFVector.hh
View file @ 3236301f
#include <list>
#include <algorithm>
#include <math.h>
#include "FixVec.h"
#include "Boundary.h"
......@@ -16,6 +17,7 @@
#include "Assembler.h"
#include "OpenMP.h"
#include "Operator.h"
#include "Parameters.h"
namespace AMDiS {
......@@ -1161,9 +1163,43 @@ namespace AMDiS {
*rIterator = (*v1Iterator) + (*v2Iterator);
};
return result;
}
template<typename T>
inline const DOFVector<T>& mod(const DOFVector<T>& v,
DOFVector<T>& result)
{
typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
for (vIterator.reset(), rIterator.reset();
!vIterator.end();
++vIterator, ++rIterator) {
*rIterator = fmod((*vIterator), 1.0);
}
return result;
}
template<typename T>
inline const DOFVector<T>& Tanh(const DOFVector<T>& v,
DOFVector<T>& result)
{
typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
double eps;
GET_PARAMETER(1, "vecphase->epsilon", "%f", &eps);
for (vIterator.reset(), rIterator.reset();
!vIterator.end();
++vIterator, ++rIterator) {
*rIterator = 0.5 * (1.0 - tanh(3.0 * (*vIterator) / eps));
}
return result;
}
template<typename T>
const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
{
......@@ -1248,4 +1284,47 @@ namespace AMDiS {
return const_cast<const T*>(result);
}
// integral u^2(1-u)^2
template<typename T>
double DOFVector<T>::DoubleWell(Quadrature* q) const
{
FUNCNAME("DOFVector::DoubleWell()");
Mesh* mesh = feSpace->getMesh();
int deg = 0;
dim = mesh->getDim();
if (!q) {
deg = 2 * feSpace->getBasisFcts()->getDegree();
q = Quadrature::provideQuadrature(dim, deg);
}
quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),
*q, INIT_PHI);
norm = 0.0;
traverseVector = const_cast<DOFVector<T>*>(this);
mesh->traverse(-1, Mesh::CALL_LEAF_EL|Mesh::FILL_COORDS|Mesh::FILL_DET, DoubleWell_fct);
return norm;
}
template<typename T>
int DOFVector<T>::DoubleWell_fct(ElInfo *elinfo)
{
double det = elinfo->getDet();
const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL,
quad_fast, NULL);
int numPoints = quad_fast->getNumPoints();
double normT = 0.0;
for (int iq = 0; iq < numPoints; iq++) {
normT += quad_fast->getWeight(iq)*sqr(uh_vec[iq])*sqr(1.0 - uh_vec[iq]);
}
norm += det * normT;
return 0;
}
}
AMDiS/src/Operator.cc
View file @ 3236301f
......@@ -251,7 +251,7 @@ namespace AMDiS {
Quadrature *quad1GrdPsi,
Quadrature *quad1GrdPhi,
Quadrature *quad0)
{
{
if (!assembler[rank])
if (optimized) {
assembler[rank] =
......@@ -280,6 +280,14 @@ namespace AMDiS {
{
vecAtQPs = subAssembler->getVectorAtQPs(vec, elInfo, quad);
}
void Vec2AtQP_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
vecAtQPs1 = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
vecAtQPs2 = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
}
void CoordsAtQP_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
......@@ -310,6 +318,14 @@ namespace AMDiS {
gradAtQPs = subAssembler->getGradientsAtQPs(vec, elInfo, quad);
}
void VecMatrixGradientAtQP_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
vecAtQPs = subAssembler->getVectorAtQPs(vec, elInfo, quad);
gradAtQPs = subAssembler->getGradientsAtQPs(vec, elInfo, quad);
}
void VecAndCoordsAtQP_SOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
......@@ -400,6 +416,15 @@ namespace AMDiS {
vecAtQPs2 = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
}
void Vec3AtQP_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
vecAtQPs1 = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
vecAtQPs2 = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
vecAtQPs3 = subAssembler->getVectorAtQPs(vec3, elInfo, quad);
}
void VecAndCoordsAtQP_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
......@@ -407,7 +432,16 @@ namespace AMDiS {
vecAtQPs = subAssembler->getVectorAtQPs(vec, elInfo, quad);
coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
}
void Vec2AndGradAtQP_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
{
vecAtQPs1 = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
vecAtQPs2 = subAssembler->getVectorAtQPs(vec2, elInfo, quad);
gradAtQPs = subAssembler->getGradientsAtQPs(vec1, elInfo, quad);
}
void FctGradientCoords_ZOT::initElement(const ElInfo* elInfo,
SubAssembler* subAssembler,
Quadrature *quad)
......@@ -476,7 +510,14 @@ namespace AMDiS {
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq]));
}
}
void Vec2AtQP_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < numPoints; iq++) {
l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
}
}
void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < numPoints; iq++) {
......@@ -506,6 +547,14 @@ namespace AMDiS {
}
}
void VecMatrixGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int numPoints,
DimMat<double> **LALt) const {
const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < numPoints; iq++) {
lalt(Lambda, (*f)(vecAtQPs[iq], gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
}
}
void VecAndCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt)
const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
for (int iq = 0; iq < numPoints; iq++) {
......@@ -644,6 +693,28 @@ namespace AMDiS {
}
}
void Vec3AtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);
}
}
void Vec3AtQP_ZOT::eval(int numPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double fac) const
{
for (int iq = 0; iq < numPoints; iq++) {
result[iq] +=
fac *
(*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) *
uhAtQP[iq];
}
}
void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
......@@ -704,6 +775,28 @@ namespace AMDiS {
}
}
void Vec2AndGradAtQP_ZOT::eval(int numPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double fac) const
{
for (int iq = 0; iq < numPoints; iq++) {
result[iq] +=
fac *
(*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]) *
uhAtQP[iq];
}
}
void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const
{
for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]);
}
}
void VecAndGradAtQP_ZOT::getC(const ElInfo *, int numPoints, double *C) const {
for (int iq = 0; iq < numPoints; iq++) {
C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
......@@ -959,6 +1052,39 @@ namespace AMDiS {
}
}
void Vec2AtQP_SOT::eval(int numPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double fac) const
{
int dow = Global::getGeo(WORLD);
if (D2UhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
double resultQP = 0.0;
for (int i = 0; i < dow; i++) {
resultQP += D2UhAtQP[iq][i][i];
}
result[iq] += fac * factor * resultQP;
}
}
}
void Vec2AtQP_SOT::weakEval(int numPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
{
if (grdUhAtQP) {
for (int iq = 0; iq < numPoints; iq++) {
double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
axpy(factor, grdUhAtQP[iq], result[iq]);
}
}
}
void CoordsAtQP_SOT::eval(int numPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
......@@ -1058,7 +1184,48 @@ namespace AMDiS {
}
}
void VecMatrixGradientAtQP_SOT::eval(int numPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double factor) const
{
int dow = Global::getGeo(WORLD);
for (int iq = 0; iq < numPoints; iq++) {
double resultQP = 0.0;
WorldMatrix<double> A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
if (D2UhAtQP) {
for (int i = 0; i < dow; i++) {
for (int j = 0; j < dow; j++) {
resultQP += A[i][j] * D2UhAtQP[iq][j][i];
}
}
}
if (grdUhAtQP) {
resultQP += (*divFct)(A) * grdUhAtQP[iq];
}
result[iq] += resultQP * factor;
}
}
void VecMatrixGradientAtQP_SOT::weakEval(int numPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const
{
if (grdUhAtQP) {
WorldMatrix<double> A;
for (int iq = 0; iq < numPoints; iq++) {
A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
result[iq] += A * grdUhAtQP[iq];
}
}
}
void VecAndCoordsAtQP_SOT::eval(int numPoints,
const double *uhAtQP,
......
AMDiS/src/Operator.h
View file @ 3236301f
......@@ -718,6 +718,77 @@ namespace AMDiS {
AbstractFunction<double, double> *f;
};
/**
* \ingroup Assembler
*
* \brief
* Laplace operator multiplied with a function evaluated at the quadrature
* points of a given DOFVector:
* \f$ f(v(\vec{x}), w(\vec{x})) \Delta u(\vec{x}) \f$
*/
class Vec2AtQP_SOT : public SecondOrderTerm {
public:
/** \brief
* Constructor.
*/
Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2,
BinaryAbstractFunction<double, double, double> *f_)
: SecondOrderTerm(f_->getDegree()),
vec1(dv1),
vec2(dv2),
f(f_)
{
setSymmetric(true);
};
/** \brief
* Implementation of \ref OperatorTerm::initElement().
*/
void initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
Quadrature *quad = NULL);
/** \brief
* Implements SecondOrderTerm::getLALt().
*/
void getLALt(const ElInfo *elInfo, int numPoints, DimMat<double> **LALt) const;
/** \brief
* Implenetation of SecondOrderTerm::eval().
*/
void eval(int numPoints,
const double *uhAtQP,
const WorldVector<double> *grdUhAtQP,
const WorldMatrix<double> *D2UhAtQP,
double *result,
double factor) const;
/** \brief
* Implenetation of SecondOrderTerm::weakEval().
*/
void weakEval(int numPoints,
const WorldVector<double> *grdUhAtQP,
WorldVector<double> *result) const;
protected:
/** \brief
* DOFVector to be evaluated at quadrature points.
*/
DOFVectorBase<double>* vec1;
DOFVectorBase<double>* vec2;
/** \brief
* Pointer to an array containing the DOFVector evaluated at quadrature
* points.
*/
double* vecAtQPs1;
double* vecAtQPs2;
/** \brief
* Function evaluated at \ref vecAtQPs.
*/
BinaryAbstractFunction<double, double, double> *f;
};