Skip to content

GitLab

A
amdis
Commit c2ebc435 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files
Options

Solved bug, when coupled problems have different basis functions degree.

parent b4ddbecb
Hide whitespace changes
Inline Side-by-side
Showing with 191 additions and 233 deletions
AMDiS/libtool
View file @ c2ebc435
......@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
# Libtool was configured on host NWRW03:
# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -328,7 +328,7 @@ link_all_deplibs=unknown
sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/"
# Run-time system search path for libraries
sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-3.0.1 /usr/lib/qt-3.3/lib "
sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
......@@ -6760,7 +6760,7 @@ build_old_libs=`case $build_libtool_libs in yes) $echo no;; *) $echo yes;; esac`
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
# Libtool was configured on host NWRW03:
# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -7041,7 +7041,7 @@ link_all_deplibs=unknown
sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/"
# Run-time system search path for libraries
sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-3.0.1 /usr/lib/qt-3.3/lib "
sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
......@@ -7065,7 +7065,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
# Libtool was configured on host NWRW03:
# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -7109,7 +7109,7 @@ LTCC="gcc"
LTCFLAGS="-g -O2"
# A language-specific compiler.
CC="f95"
CC="g77"
# Is the compiler the GNU C compiler?
with_gcc=yes
......@@ -7346,10 +7346,10 @@ variables_saved_for_relink="PATH LD_LIBRARY_PATH LD_RUN_PATH GCC_EXEC_PREFIX COM
link_all_deplibs=unknown
# Compile-time system search path for libraries
sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../i386-redhat-linux/4.1.2/ /usr/lib/gcc/i386-redhat-linux/4.1.2/../../../ /lib/i386-redhat-linux/4.1.2/ /lib/ /usr/lib/i386-redhat-linux/4.1.2/ /usr/lib/"
sys_lib_search_path_spec=" /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../../i386-redhat-linux/lib/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../i386-redhat-linux/3.4.6/ /usr/lib/gcc/i386-redhat-linux/3.4.6/../../../ /lib/i386-redhat-linux/3.4.6/ /lib/ /usr/lib/i386-redhat-linux/3.4.6/ /usr/lib/"
# Run-time system search path for libraries
sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-3.0.1 /usr/lib/qt-3.3/lib "
sys_lib_dlsearch_path_spec="/lib /usr/lib /usr/lib/octave-2.9.9 /usr/lib/qt-3.3/lib "
# Fix the shell variable $srcfile for the compiler.
fix_srcfile_path=""
......
AMDiS/src/Element.h
View file @ c2ebc435
......@@ -328,7 +328,6 @@ namespace AMDiS {
/** \} */
// ===== pure virtual methods =================================================
/** \name pure virtual methods
* \{
......@@ -339,12 +338,9 @@ namespace AMDiS {
* Used by Estimator for the jumps => same quadrature nodes from both sides!
*/
virtual const FixVec<int,WORLD>&
sortFaceIndices(int face, FixVec<int,WORLD> *vec) const = 0;
sortFaceIndices(int face, FixVec<int, WORLD> *vec) const = 0;
/** \brief
* Returns a copy of itself. Needed by Mesh to create Elements by a
* prototype.
*/
/// Returns a copy of itself. Needed by Mesh to create Elements by a prototype.
virtual Element *clone() = 0;
/** \brief
......
AMDiS/src/FiniteElemSpace.cc
View file @ c2ebc435
......@@ -10,12 +10,12 @@ namespace AMDiS {
std::vector<FiniteElemSpace*> FiniteElemSpace::feSpaces;
FiniteElemSpace::FiniteElemSpace(DOFAdmin* admin_,
const BasisFunction* bas_fcts_,
const BasisFunction* bas_fcts,
Mesh* aMesh,
const std::string& aString)
: name(aString),
admin(admin_),
basFcts(bas_fcts_),
basFcts(bas_fcts),
mesh(aMesh)
{
FUNCNAME("FiniteElemSpace::FiniteElemSpace()");
......@@ -30,17 +30,16 @@ namespace AMDiS {
for (int i = 0; i < mesh->getNumberOfDOFAdmin(); i++) {
admin_local = &(mesh->getDOFAdmin(i));
int j = 0;
for (; j <= mesh->getDim(); j++) {
for (; j <= mesh->getDim(); j++)
if (admin_local->getNumberOfDOFs(j) != (*ndof)[j])
break;
}
if (j > mesh->getDim())
break;
admin_local = NULL;
}
if (!admin_local) {
if (!admin_local)
admin_local = mesh->createDOFAdmin(name, *ndof);
}
admin = const_cast<DOFAdmin*>(admin_local);
}
......@@ -76,24 +75,13 @@ namespace AMDiS {
{
int numSpaces = static_cast<int>(feSpaces.size());
for (int i = 0; i < numSpaces; i++) {
if (admin) {
if (feSpaces[i]->admin == admin) {
return feSpaces[i];
}
} else {
if (feSpaces[i]->basFcts == basFcts && feSpaces[i]->mesh == mesh) {
return feSpaces[i];
}
}
}
for (int i = 0; i < numSpaces; i++)
if (feSpaces[i]->basFcts == basFcts &&
feSpaces[i]->mesh == mesh &&
(!admin || (admin && feSpaces[i]->admin == admin)))
return feSpaces[i];
if (admin) {
ERROR_EXIT("no fespace found for this admin\n");
return NULL;
} else {
return new FiniteElemSpace(NULL, basFcts, mesh, name_);
}
return new FiniteElemSpace(admin, basFcts, mesh, name_);
}
int FiniteElemSpace::calcMemoryUsage()
......
AMDiS/src/Lagrange.cc
View file @ c2ebc435
......@@ -1046,11 +1046,14 @@ namespace AMDiS {
int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(CENTER);
Element *el = list->getElement(0);
int node = drv->getFESpace()->getMesh()->getNode(CENTER);
DegreeOfFreedom dof0 = el->getDOF(node, n0); /* Parent center */
DegreeOfFreedom dof_new = el->getChild(0)->getDOF(node, n0); /* newest vertex is center */
// Parent center
DegreeOfFreedom dof0 = el->getDOF(node, n0);
// Newest vertex is center
DegreeOfFreedom dof_new = el->getChild(0)->getDOF(node, n0);
(*drv)[dof_new] = (*drv)[dof0];
dof_new = el->getChild(1)->getDOF(node, n0); /* newest vertex is center */
// Newest vertex is center
dof_new = el->getChild(1)->getDOF(node, n0);
(*drv)[dof_new] = (*drv)[dof0];
}
......
AMDiS/src/MacroReader.cc
View file @ c2ebc435
......@@ -158,9 +158,8 @@ namespace AMDiS {
associated = new VertexVector(mesh->getVertexAdmin(), "vertex vector");
mesh->periodicAssociations[boundaryType] = associated;
VertexVector::Iterator it(associated, ALL_DOFS);
for (it.reset2(); !it.end(); ++it) {
for (it.reset2(); !it.end(); ++it)
*it = it.getDOFIndex();
}
}
for (int j = 0; j < dim; j++) {
......@@ -213,11 +212,10 @@ namespace AMDiS {
}
}
for (int i = 0; i < mesh->getNumberOfVertices(); i++) {
if (periodicMap.getEntry(i) != -1) {
MSG("identification : vertex %d is now vertex %d\n", i, periodicMap.getEntry(i));
}
}
for (int i = 0; i < mesh->getNumberOfVertices(); i++)
if (periodicMap.getEntry(i) != -1)
MSG("identification : vertex %d is now vertex %d\n",
i, periodicMap.getEntry(i));
}
// =========================================================
......@@ -268,7 +266,6 @@ namespace AMDiS {
boundaryDOFs(mesh);
// initial boundary projections
//if(dim > 1) {
int numFaces = mesh->getGeo(FACE);
int dim = mesh->getDim();
mel = mesh->firstMacroElement();
......@@ -291,7 +288,6 @@ namespace AMDiS {
}
}
}
//}
macroInfo->fillBoundaryInfo(mesh);
......@@ -1409,46 +1405,35 @@ namespace AMDiS {
&(nei->projection[mesh->getGeo(FACE)+edge_no]);
if (mesh->getNumberOfDOFs(EDGE))
{
// if(periodic) {
// nei->element->setDOF(node+edge_no, mesh->getDOF(EDGE));
// } else {
nei->element->setDOF(node+edge_no,edge_dof);
// }
}
nei->element->setDOF(node+edge_no,edge_dof);
if (next_el[edge_no][0] != opp_v)
{
if(nei->getBoundary(next_el[edge_no][0])) {
lbound = newBound(nei->getBoundary(next_el[edge_no][0]), lbound);
Projection *neiProject = nei->getProjection(next_el[edge_no][0]);
if(!lproject)
if (next_el[edge_no][0] != opp_v) {
if(nei->getBoundary(next_el[edge_no][0])) {
lbound = newBound(nei->getBoundary(next_el[edge_no][0]), lbound);
Projection *neiProject = nei->getProjection(next_el[edge_no][0]);
if (!lproject) {
lproject = neiProject;
} else {
if (neiProject && (lproject->getID() < neiProject->getID()))
lproject = neiProject;
else {
if(neiProject && (lproject->getID() < neiProject->getID())) {
lproject = neiProject;
}
}
}
opp_v = nei->getOppVertex(next_el[edge_no][0]);
nei = nei->getNeighbour(next_el[edge_no][0]);
}
else
{
if(nei->getBoundary(next_el[edge_no][1])) {
lbound = newBound(nei->getBoundary(next_el[edge_no][1]), lbound);
Projection *neiProject = nei->getProjection(next_el[edge_no][1]);
if(!lproject)
opp_v = nei->getOppVertex(next_el[edge_no][0]);
nei = nei->getNeighbour(next_el[edge_no][0]);
} else {
if (nei->getBoundary(next_el[edge_no][1])) {
lbound = newBound(nei->getBoundary(next_el[edge_no][1]), lbound);
Projection *neiProject = nei->getProjection(next_el[edge_no][1]);
if (!lproject) {
lproject = neiProject;
} else {
if (neiProject && (lproject->getID() < neiProject->getID()))
lproject = neiProject;
else {
if(neiProject && (lproject->getID() < neiProject->getID())) {
lproject = neiProject;
}
}
}
opp_v = nei->getOppVertex(next_el[edge_no][1]);
nei = nei->getNeighbour(next_el[edge_no][1]);
}
opp_v = nei->getOppVertex(next_el[edge_no][1]);
nei = nei->getNeighbour(next_el[edge_no][1]);
}
}
if (!nei)
{
......@@ -1861,10 +1846,10 @@ namespace AMDiS {
{
FUNCNAME("MacroReader::checkMesh()");
int i, nused, nfree;
int i, nused, nfree;
DOFAdmin *localAdmin=mesh->admin[0];
Flag fill_flag;
int error_detected = 0;
Flag fill_flag;
int error_detected = 0;
MSG("checking mesh\n");
......
AMDiS/src/Mesh.cc
View file @ c2ebc435
......@@ -387,14 +387,15 @@ namespace AMDiS {
localAdmin->setMesh(this);
std::vector<DOFAdmin*>::iterator dai = std::find(admin.begin(),admin.end(),localAdmin);
std::vector<DOFAdmin*>::iterator dai =
std::find(admin.begin(),admin.end(),localAdmin);
TEST_EXIT(dai == admin.end())
("admin %s is already associated to mesh %s\n",
localAdmin->getName().c_str(), this->getName().c_str());
// if this will be required, see the untested code in revision < 224
TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work yet!\n");
// TEST_EXIT(!initialized)("Adding DOFAdmins to initilized meshes does not work yet!\n");
admin.push_back(localAdmin);
......@@ -435,9 +436,8 @@ namespace AMDiS {
}
node[CENTER] = nNodeEl;
if (nDOF[CENTER] > 0) {
if (nDOF[CENTER] > 0)
nNodeEl += 1;
}
}
void Mesh::dofCompress()
......@@ -1082,25 +1082,25 @@ namespace AMDiS {
// If there is no value file which should be written, we can delete
// the information of the macro file.
if (!valueFilename.length()) {
if (!valueFilename.length())
clearMacroFileInfo();
}
}
initialized = true;
}
bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
bool Mesh::associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2)
{
std::map<BoundaryType, VertexVector*>::iterator it;
std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
for (it = periodicAssociations.begin(); it != end; ++it) {
for (it = periodicAssociations.begin(); it != end; ++it)
if ((*(it->second))[dof1] == dof2)
return true;
}
return false;
}
bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2) {
bool Mesh::indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2)
{
std::vector<DegreeOfFreedom> associatedToDOF1;
std::map<BoundaryType, VertexVector*>::iterator it;
std::map<BoundaryType, VertexVector*>::iterator end = periodicAssociations.end();
......
AMDiS/src/ProblemVec.cc
View file @ c2ebc435
......@@ -72,17 +72,9 @@ namespace AMDiS {
WARNING("feSpaces already created\n");
} else {
if (initFlag.isSet(INIT_FE_SPACE) ||
(initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
(initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE)))
createFESpace(NULL);
if (!adoptFlag.isSet(INIT_MESH)) {
createFESpace(NULL);
} else {
TEST_EXIT(meshes.size() == 1)("I have to think about it!\n");
TEST_EXIT(meshes[0]->getNumberOfDOFAdmin() == 1)("I have to think about it!\n");
createFESpace(&(const_cast<DOFAdmin&>(meshes[0]->getDOFAdmin(0))));
}
}
if (adoptProblem &&
(adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
feSpaces = adoptProblem->getFESpaces();
......@@ -106,9 +98,9 @@ namespace AMDiS {
WARNING("no feSpace created\n");
// === create system ===
if (initFlag.isSet(INIT_SYSTEM)) {
if (initFlag.isSet(INIT_SYSTEM))
createMatricesAndVectors();
}
if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
solution = adoptProblem->getSolution();
rhs = adoptProblem->getRHS();
......@@ -189,9 +181,13 @@ namespace AMDiS {
&globalRefinements);
// Initialize the meshes if there is no serialization file.
for (int i = 0; i < static_cast<int>(meshes.size()); i++)
if (initFlag.isSet(INIT_MESH) && meshes[i] && !(meshes[i]->isInitialized()))
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("");
......@@ -233,9 +229,9 @@ namespace AMDiS {
sprintf(number, "%d", i);
int refSet = -1;
GET_PARAMETER(0, name + "->refinement set[" + number + "]", "%d", &refSet);
if (refSet < 0) {
if (refSet < 0)
refSet = 0;
}
if (meshForRefinementSet[refSet] == NULL) {
Mesh *newMesh = new Mesh(meshName, dim);
meshForRefinementSet[refSet] = newMesh;
......@@ -296,9 +292,9 @@ namespace AMDiS {
}
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
for (int j = 0; j < nComponents; j++)
traverseInfo.getMatrix(i, j).setFESpace(componentSpaces[i], componentSpaces[j]);
}
traverseInfo.getVector(i).setFESpace(componentSpaces[i]);
}
......
AMDiS/src/ResidualEstimator.cc
View file @ c2ebc435
......@@ -12,9 +12,11 @@ namespace AMDiS {
: Estimator(name, r),
C0(1.0),
C1(1.0),
C2(1.0),
C2(0.0),
C3(1.0)
{
FUNCNAME("ResidualEstimator::ResidualEstimator()");
GET_PARAMETER(0, name + "->C0", "%f", &C0);
GET_PARAMETER(0, name + "->C1", "%f", &C1);
GET_PARAMETER(0, name + "->C2", "%f", &C2);
......@@ -24,6 +26,8 @@ namespace AMDiS {
C1 = C1 > 1.e-25 ? sqr(C1) : 0.0;
C2 = C2 > 1.e-25 ? sqr(C2) : 0.0;
C3 = C3 > 1.e-25 ? sqr(C3) : 0.0;
TEST_EXIT(C2 == 0.0)("C2 is not used! Please remove it or set it to 0.0!\n");
}
void ResidualEstimator::init(double ts)
......@@ -46,22 +50,19 @@ namespace AMDiS {
basFcts[system] = uh[system]->getFESpace()->getBasisFcts();
degree = std::max(degree, basFcts[system]->getDegree());
}
degree *= 2;
quad = Quadrature::provideQuadrature(dim, degree);
nPoints = quad->getNumPoints();
Flag flag = INIT_PHI | INIT_GRD_PHI;
if (degree > 2) {
flag |= INIT_D2_PHI;
}
if (degree > 2)
flag |= INIT_D2_PHI;
for (int system = 0; system < nSystems; system++) {
for (int system = 0; system < nSystems; system++)
quadFast[system] = FastQuadrature::provideFastQuadrature(basFcts[system],
*quad,
flag);
}
flag);
uhEl = new double*[nSystems];
uhNeigh = new double*[nSystems];
......@@ -76,9 +77,7 @@ namespace AMDiS {
uhQP = timestep ? new double[nPoints] : NULL;
uhOldQP = timestep ? new double[nPoints] : NULL;
riq = new double[nPoints];
grdUh_qp = NULL;
D2uhqp = NULL;
......@@ -110,16 +109,16 @@ namespace AMDiS {
neighInfo = mesh->createNewElInfo();
// prepare date for computing jump residual
if (C1 && (dim > 1)) {
surfaceQuad_ = Quadrature::provideQuadrature(dim - 1, degree);
nPointsSurface_ = surfaceQuad_->getNumPoints();
grdUhEl_.resize(nPointsSurface_);
grdUhNeigh_.resize(nPointsSurface_);
jump_.resize(nPointsSurface_);
localJump_.resize(nPointsSurface_);
neighbours_ = Global::getGeo(NEIGH, dim);
lambdaNeigh_ = new DimVec<WorldVector<double> >(dim, NO_INIT);
lambda_ = new DimVec<double>(dim, NO_INIT);
if (C1 > 0.0 && dim > 1) {
surfaceQuad = Quadrature::provideQuadrature(dim - 1, degree);
nPointsSurface = surfaceQuad->getNumPoints();
grdUhEl.resize(nPointsSurface);
grdUhNeigh.resize(nPointsSurface);
jump.resize(nPointsSurface);
localJump.resize(nPointsSurface);
nNeighbours = Global::getGeo(NEIGH, dim);
lambdaNeigh = new DimVec<WorldVector<double> >(dim, NO_INIT);
lambda = new DimVec<double>(dim, NO_INIT);
}
}
......@@ -165,8 +164,8 @@ namespace AMDiS {
delete [] D2uhqp;
if (C1 && (dim > 1)) {
delete lambdaNeigh_;
delete lambda_;
delete lambdaNeigh;
delete lambda;
}
delete neighInfo;
......@@ -187,43 +186,40 @@ namespace AMDiS {
double est_el = el->getEstimation(row);
double h2 = h2_from_det(det, dim);
for (int iq = 0; iq < nPoints; iq++) {
for (int iq = 0; iq < nPoints; iq++)
riq[iq] = 0.0;
}
for (int system = 0; system < nSystems; system++) {
if (matrix[system] == NULL)
continue;
// init assemblers
// === init assemblers ===
for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
++it, ++itfac) {
if (*itfac == NULL || **itfac != 0.0) {
++it, ++itfac)
if (*itfac == NULL || **itfac != 0.0)
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
}
}
if (C0) {
if (C0 > 0.0)
for (it = const_cast<DOFVector<double>*>(fh[system])->getOperatorsBegin();
it != const_cast<DOFVector<double>*>(fh[system])->getOperatorsEnd();
++it) {
++it)
(*it)->getAssembler(omp_get_thread_num())->initElement(elInfo, NULL, quad);
}
}
if (timestep && uhOld[system]) {
TEST_EXIT_DBG(uhOld[system])("no uhOld\n");
uhOld[system]->getLocalVector(el, uhOldEl[system]);
// ===== time and element residuals
if (C0 || C3) {
// === Compute time error. ===
if (C0 > 0.0 || C3 > 0.0) {
uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
uhOld[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhOldQP);
if (C3 && uhOldQP && system == std::max(row, 0)) {
if (C3 > 0.0 && uhOldQP && system == std::max(row, 0)) {
val = 0.0;
for (int iq = 0; iq < nPoints; iq++) {
double tiq = (uhQP[iq] - uhOldQP[iq]);
......@@ -236,27 +232,32 @@ namespace AMDiS {
}
}
if (C0) {
// === Compute element residual. ===
if (C0 > 0.0) {
for (it = const_cast<DOFMatrix*>(matrix[system])->getOperatorsBegin(),
itfac = const_cast<DOFMatrix*>(matrix[system])->getOperatorEstFactorBegin();
it != const_cast<DOFMatrix*>(matrix[system])->getOperatorsEnd();
++it, ++itfac) {
if (*itfac == NULL || **itfac != 0.0) {
if ((uhQP == NULL) && (*it)->zeroOrderTerms()) {
if (uhQP == NULL && (*it)->zeroOrderTerms()) {
uhQP = new double[nPoints];
uh[system]->getVecAtQPs(elInfo, NULL, quadFast[system], uhQP);
}
if ((grdUh_qp == NULL) && ((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
if (grdUh_qp == NULL &&
((*it)->firstOrderTermsGrdPsi() || (*it)->firstOrderTermsGrdPhi())) {
grdUh_qp = new WorldVector<double>[nPoints];
uh[system]->getGrdAtQPs(elInfo, NULL, quadFast[system], grdUh_qp);
}
if ((D2uhqp == NULL) && (degree > 2) && (*it)->secondOrderTerms()) {
if (D2uhqp == NULL && degree > 2 && (*it)->secondOrderTerms()) {
D2uhqp = new WorldMatrix<double>[nPoints];
uh[system]->getD2AtQPs(elInfo, NULL, quadFast[system], D2uhqp);
}
}
}
// === Compute the element residual and store it in irq. ===
r(elInfo,
nPoints,
uhQP,
......@@ -284,17 +285,19 @@ namespace AMDiS {
est_el += val;
// ===== jump residuals
// === Compute jump residuals. ===
if (C1 && (dim > 1)) {
int dow = Global::getGeo(WORLD);
for (int face = 0; face < neighbours_; face++) {
for (int face = 0; face < nNeighbours; face++) {
Element *neigh = const_cast<Element*>(elInfo->getNeighbour(face));
if (neigh && neigh->getMark()) {
int oppV = elInfo->getOppVertex(face);
el->sortFaceIndices(face, &faceIndEl_);
neigh->sortFaceIndices(oppV, &faceIndNeigh_);
el->sortFaceIndices(face, &faceIndEl);