Commit f8646256 authored by Thomas Witkowski's avatar Thomas Witkowski

* Merge of PardisoSolver and parallel assembler

parent d56091a2
......@@ -127,6 +127,8 @@ $(SOURCE_DIR)/Element.h $(SOURCE_DIR)/ElementConnection.h \
$(SOURCE_DIR)/ElInfo.h $(SOURCE_DIR)/ElInfo1d.h $(SOURCE_DIR)/ElInfo2d.h $(SOURCE_DIR)/ElInfo3d.h \
$(SOURCE_DIR)/Error.h $(SOURCE_DIR)/Error.hh \
$(SOURCE_DIR)/Estimator.h $(SOURCE_DIR)/Estimator.cc \
$(SOURCE_DIR)/ResidualEstimator.h $(SOURCE_DIR)/ResidualEstimator.cc \
$(SOURCE_DIR)/ResidualParallelEstimator.h $(SOURCE_DIR)/ResidualParallelEstimator.cc \
$(SOURCE_DIR)/FiniteElemSpace.h \
$(SOURCE_DIR)/FixVec.h $(SOURCE_DIR)/FixVec.hh \
$(SOURCE_DIR)/FixVecConvert.h \
......@@ -134,7 +136,7 @@ $(SOURCE_DIR)/Flag.h \
$(SOURCE_DIR)/Global.h \
$(SOURCE_DIR)/GMResSolver.h $(SOURCE_DIR)/GMResSolver.hh \
$(SOURCE_DIR)/GMResSolver2.h $(SOURCE_DIR)/GMResSolver2.hh \
$(SOURCE_DIR)/TFQMR.h \ $(SOURCE_DIR)/TFQMR.hh \
$(SOURCE_DIR)/TFQMR.h $(SOURCE_DIR)/TFQMR.hh \
$(SOURCE_DIR)/VecSymSolver.h $(SOURCE_DIR)/VecSymSolver.hh \
$(SOURCE_DIR)/UmfPackSolver.h $(SOURCE_DIR)/UmfPackSolver.hh \
$(SOURCE_DIR)/PardisoSolver.h $(SOURCE_DIR)/PardisoSolver.hh $(SOURCE_DIR)/PardisoSolver.cc \
......@@ -152,7 +154,9 @@ $(SOURCE_DIR)/Preconditioner.h \
$(SOURCE_DIR)/Quadrature.h \
$(SOURCE_DIR)/RCNeighbourList.h \
$(SOURCE_DIRe)/RefinementManager.h \
$(SOURCE_DIR)/RefinementManager1d.h $(SOURCE_DIR)/RefinementManager2d.h $(SOURCE_DIR)/RefinementManager3d.h \
$(SOURCE_DIR)/RefinementManager1d.h \
$(SOURCE_DIR)/RefinementManager2d.h \
$(SOURCE_DIR)/RefinementManager3d.h \
$(SOURCE_DIR)/TecPlotWriter.h $(SOURCE_DIR)/TecPlotWriter.hh \
$(SOURCE_DIR)/Tetrahedron.h \
$(SOURCE_DIR)/Traverse.h \
......@@ -177,7 +181,9 @@ $(SOURCE_DIR)/demangle.cc \
$(SOURCE_DIR)/DOFAdmin.cc \
$(SOURCE_DIR)/DOFMatrix.cc \
$(SOURCE_DIR)/Element.cc \
$(SOURCE_DIR)/ElInfo1d.cc $(SOURCE_DIR)/ElInfo2d.cc $(SOURCE_DIR)/ElInfo3d.cc \
$(SOURCE_DIR)/ElInfo1d.cc \
$(SOURCE_DIR)/ElInfo2d.cc \
$(SOURCE_DIR)/ElInfo3d.cc \
$(SOURCE_DIR)/FiniteElemSpace.cc \
$(SOURCE_DIR)/FixVec.cc \
$(SOURCE_DIR)/Flag.cc \
......@@ -196,6 +202,7 @@ $(SOURCE_DIR)/RefinementManager2d.cc \
$(SOURCE_DIR)/RefinementManager3d.cc \
$(SOURCE_DIR)/Tetrahedron.cc \
$(SOURCE_DIR)/Traverse.cc \
$(SOURCE_DIR)/TraverseParallel.h $(SOURCE_DIR)/TraverseParallel.cc \
$(SOURCE_DIR)/Triangle.cc \
$(SOURCE_DIR)/TecPlotWriter.cc \
$(SOURCE_DIR)/ValueWriter.cc \
......@@ -211,11 +218,8 @@ $(SOURCE_DIR)/OpenMP.h
COMPOSITE_SOURCE_DIR = ../compositeFEM/src
if AMDIS_DEBUG
libcompositeFEM_la_CXXFLAGS = -g -O0 -Wall -DDEBUG=1 -I$(SOURCE_DIR)
else
libcompositeFEM_la_CXXFLAGS = -O2 -Wall -DDEBUG=0 -I$(SOURCE_DIR)
endif
libcompositeFEM_la_CXXFLAGS = $(libamdis_la_CXXFLAGS)
libcompositeFEM_la_SOURCES = $(COMPOSITE_SOURCE_DIR)/CFE_Integration.h \
$(COMPOSITE_SOURCE_DIR)/CFE_Integration.cc \
......
This diff is collapsed.
......@@ -2,6 +2,7 @@
#include "ElementMatrix.h"
#include "ElementVector.h"
#include "OpenMP.h"
#include "SubElementAssembler.h"
#include "SubElInfo.h"
......@@ -116,9 +117,11 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
subElementAssembler->getNCol());
subPolMat2->set(0.0);
if(!assembler) {
assembler = NEW StandardAssembler(this, NULL, NULL, NULL, NULL,
rowFESpace, colFESpace);
int myRank = omp_get_thread_num();
if (!assembler[myRank]) {
assembler[myRank] =
NEW StandardAssembler(this, NULL, NULL, NULL, NULL, rowFESpace, colFESpace);
}
if (elLS->getLevelSetDomain() ==
......@@ -129,7 +132,7 @@ CompositeFEMOperator::getElementMatrix(const ElInfo *elInfo,
elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
}
assembler->calculateElementMatrix(elInfo, elMat, 1.0);
assembler[myRank]->calculateElementMatrix(elInfo, elMat, 1.0);
subElementAssembler->getSubPolytopeMatrix(subPolytope,
subElementAssembler,
elInfo,
......@@ -259,20 +262,21 @@ CompositeFEMOperator::getElementVector(const ElInfo *elInfo,
subPolVec2 = NEW ElementVector(subElementAssembler->getNRow());
subPolVec2->set(0.0);
if(!assembler) {
assembler = NEW StandardAssembler(this, NULL, NULL, NULL, NULL,
rowFESpace, colFESpace);
int myRank = omp_get_thread_num();
if (!assembler[myRank]) {
assembler[myRank] =
NEW StandardAssembler(this, NULL, NULL, NULL, NULL, rowFESpace, colFESpace);
}
if (elLS->getLevelSetDomain() ==
ElementLevelSet::LEVEL_SET_INTERIOR) {
elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_EXTERIOR);
}
else {
} else {
elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
}
assembler->calculateElementVector(elInfo, elVec, 1.0);
assembler[myRank]->calculateElementVector(elInfo, elVec, 1.0);
subElementAssembler->getSubPolytopeVector(subPolytope,
subElementAssembler,
elInfo,
......
......@@ -52,7 +52,7 @@ WARN_LOGFILE =
# configuration options related to the input files
#---------------------------------------------------------------------------
INPUT = ../src/
FILE_PATTERNS = *.h *.hh *.cc
FILE_PATTERNS = *.h
RECURSIVE = NO
EXCLUDE =
EXCLUDE_PATTERNS =
......
......@@ -44,7 +44,7 @@ available_tags=" CXX F77"
# ### BEGIN LIBTOOL CONFIG
# Libtool was configured on host NWRW09:
# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -7266,7 +7266,7 @@ disable_libs=static
# End:
# ### BEGIN LIBTOOL TAG CONFIG: CXX
# Libtool was configured on host NWRW09:
# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......@@ -7574,7 +7574,7 @@ include_expsyms=""
# ### BEGIN LIBTOOL TAG CONFIG: F77
# Libtool was configured on host NWRW09:
# Libtool was configured on host NWRW15:
# Shell to use when invoking shell scripts.
SHELL="/bin/sh"
......
......@@ -63,7 +63,7 @@ namespace AMDiS {
coarsenAllowed(0),
refinementAllowed(1),
refineBisections(1),
coarseBisections(1)
coarseBisections(1)
{
double totalTol = 1.0;
double relSpaceErr = 1.0;
......
This diff is collapsed.
......@@ -118,7 +118,7 @@ namespace AMDiS {
* Returns \ref terms
*/
inline ::std::vector<OperatorTerm*> *getTerms() {
return &terms;
return &terms[omp_get_thread_num()];
};
/** \brief
......@@ -273,7 +273,7 @@ namespace AMDiS {
/** \brief
* List of all terms with a contribution to this SubAssembler
*/
::std::vector<OperatorTerm*> terms;
::std::vector< ::std::vector<OperatorTerm*> > terms;
/** \brief
*
......
......@@ -69,12 +69,13 @@ namespace AMDiS {
const WorldVector<double>& BasisFunction::evalGrdUh(const DimVec<double>& lambda,
const DimVec<WorldVector<double> >& grd_lambda,
const double *uh_loc, WorldVector<double>* grd_uh) const
const double *uh_loc,
WorldVector<double>* grd_uh) const
{
static WorldVector<double> grd;
static WorldVector<double> grd;
DimVec<double> grd_b(dim, DEFAULT_VALUE, 0.);
WorldVector<double> *val;
DimVec<double> grd_b(dim, DEFAULT_VALUE, 0.0);
WorldVector<double> *val;
DimVec<double> grd1(dim, DEFAULT_VALUE, 0.);
val = grd_uh ? grd_uh : &grd;
......@@ -83,7 +84,8 @@ namespace AMDiS {
grd1[j] = 0.0;
for (int i = 0; i < nBasFcts; i++) {
grd_b = (*(*grdPhi)[i])(lambda);
(*(*grdPhi)[i])(lambda, grd_b);
for (int j = 0; j < dim + 1; j++)
grd1[j] += uh_loc[i] * grd_b[j];
}
......@@ -96,7 +98,7 @@ namespace AMDiS {
(*val)[i] += grd_lambda[j][i] * grd1[j];
}
return((*val));
return ((*val));
}
const WorldMatrix<double>& BasisFunction::evalD2Uh(const DimVec<double>& lambda,
......@@ -104,14 +106,12 @@ namespace AMDiS {
const double *uh_loc, WorldMatrix<double>* D2_uh) const
{
static WorldMatrix<double> D2(DEFAULT_VALUE, 0.);
DimMat<double> D2_b(dim, NO_INIT);
WorldMatrix<double> *val;
DimMat<double> D2_b(dim, DEFAULT_VALUE, 0.0);
DimMat<double> D2_tmp(dim, DEFAULT_VALUE, 0.0);
val = D2_uh ? D2_uh : &D2;
WorldMatrix<double> *val = D2_uh ? D2_uh : &D2;
for (int i = 0; i < nBasFcts; i++) {
D2_b = (*(*d2Phi)[i])(lambda);
(*(*d2Phi)[i])(lambda, D2_b);
for (int k = 0; k < dim + 1; k++)
for (int l = 0; l < dim + 1; l++)
D2_tmp[k][l] += uh_loc[i] * D2_b[k][l];
......@@ -122,12 +122,12 @@ namespace AMDiS {
for (int i = 0; i < dow; i++)
for (int j = 0; j < dow; j++) {
(*val)[i][j] = 0.0;
for (int k = 0; k < dim+1; k++)
for (int l = 0; l < dim+1; l++)
for (int k = 0; k < dim + 1; k++)
for (int l = 0; l < dim + 1; l++)
(*val)[i][j] += grd_lambda[k][i] * grd_lambda[l][j] * D2_tmp[k][l];
}
return((*val));
return ((*val));
}
......
......@@ -38,7 +38,6 @@ namespace AMDiS {
// ============================================================================
class DOFAdmin;
//template<typename T> class DOFMatrix;
class Element;
class ElInfo;
class RCNeighbourList;
......@@ -54,17 +53,64 @@ namespace AMDiS {
template <typename T, GeoIndex d> class FixVec;
template <typename T> class VectorOfFixVecs;
// ============================================================================
// ===== function interfaces===================================================
// ============================================================================
/** \brief
* Function interface for evaluating basis functions.
*/
class BasFctType
{
public:
BasFctType() {};
virtual ~BasFctType() {};
virtual double operator()(const DimVec<double>&) const = 0;
};
/** \brief
* Function interface for evaluating gradients of basis functions.
*/
class GrdBasFctType
{
public:
GrdBasFctType() {};
virtual ~GrdBasFctType() {};
virtual void operator()(const DimVec<double>&,
DimVec<double>&) const = 0;
};
/** \brief
* Function interface for evaluating second derivative of basis functions.
*/
class D2BasFctType
{
public:
D2BasFctType() {};
virtual ~D2BasFctType() {};
virtual void operator()(const DimVec<double>&,
DimMat<double>&) const = 0;
};
// ============================================================================
// ===== typedefs =============================================================
// ============================================================================
typedef AbstractFunction<double, DimVec<double> > BasFctType;
typedef AbstractFunction<DimVec<double>, DimVec<double> > GrdBasFctType;
typedef AbstractFunction<DimMat<double>, DimVec<double> > D2BasFctType;
typedef BasFctType *BFptr;
typedef GrdBasFctType *GBFptr;
typedef D2BasFctType *DBFptr;
// ============================================================================
// ===== class BasisFunction ==================================================
// ============================================================================
......
......@@ -27,28 +27,27 @@ namespace AMDiS {
{
// ===== fill local conditions ==============================================
const FiniteElemSpace *feSpace = vec->getFESpace();
const BoundaryType *localBound = NULL;
const DegreeOfFreedom *dofIndices = NULL;
const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber();
DOFAdmin *admin = feSpace->getAdmin();
Vector<DegreeOfFreedom> dofIndices;
const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber();
DOFAdmin *admin = feSpace->getAdmin();
::std::map<BoundaryType, BoundaryCondition*>::iterator it;
if (localBCs.size() > 0) {
// get boundaries of all DOFs
localBound = basisFcts->getBound(elInfo, NULL);
BoundaryType *localBound = GET_MEMORY(BoundaryType, nBasFcts);
basisFcts->getBound(elInfo, localBound);
// get dof indices
dofIndices = basisFcts->getLocalIndices(elInfo->getElement(),
admin, NULL);
basisFcts->getLocalIndicesVec(elInfo->getElement(), admin, &dofIndices);
// apply non dirichlet boundary conditions
for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) {
if (!(*it).second->isDirichlet()) {
(*it).second->fillBoundaryCondition(vec, elInfo, dofIndices, localBound, nBasFcts);
(*it).second->fillBoundaryCondition(vec, elInfo, &dofIndices[0], localBound, nBasFcts);
}
}
}
......@@ -57,10 +56,12 @@ namespace AMDiS {
for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) {
if ((*it).second->isDirichlet()) {
(*it).second->fillBoundaryCondition(vec, elInfo, dofIndices, localBound, nBasFcts);
(*it).second->fillBoundaryCondition(vec, elInfo, &dofIndices[0], localBound, nBasFcts);
}
}
}
FREE_MEMORY(localBound, BoundaryType, nBasFcts);
}
}
......@@ -69,38 +70,37 @@ namespace AMDiS {
{
// ===== fill local conditions ==============================================
const FiniteElemSpace *feSpace = mat->getRowFESpace();
const BoundaryType *localBound = NULL;
const DegreeOfFreedom *dofIndices = NULL;
const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber();
DOFAdmin *admin = feSpace->getAdmin();
Vector<DegreeOfFreedom> dofIndices;
const BasisFunction *basisFcts = feSpace->getBasisFcts();
int nBasFcts = basisFcts->getNumber();
DOFAdmin *admin = feSpace->getAdmin();
::std::map<BoundaryType, BoundaryCondition*>::iterator it;
if (localBCs.size() > 0) {
// get boundaries of all DOFs
localBound = basisFcts->getBound(elInfo, NULL);
const BoundaryType *localBound = basisFcts->getBound(elInfo, NULL);
// get dof indices
dofIndices = basisFcts->getLocalIndices(elInfo->getElement(), admin, NULL);
basisFcts->getLocalIndicesVec(elInfo->getElement(), admin, &dofIndices);
// apply non dirichlet boundary conditions
for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) {
if (!(*it).second->isDirichlet()) {
(*it).second->fillBoundaryCondition(mat, elInfo, dofIndices, localBound, nBasFcts);
(*it).second->fillBoundaryCondition(mat, elInfo, &dofIndices[0], localBound, nBasFcts);
}
}
}
// apply dirichlet boundary conditions
for (it = localBCs.begin(); it != localBCs.end(); ++it) {
if ((*it).second) {
if ((*it).second->isDirichlet()) {
(*it).second->fillBoundaryCondition(mat, elInfo, dofIndices, localBound, nBasFcts);
(*it).second->fillBoundaryCondition(mat, elInfo, &dofIndices[0], localBound, nBasFcts);
}
}
}
}
}
}
......
......@@ -98,17 +98,6 @@ namespace AMDiS {
mesh = aMesh;
// if(mesh->getDim() == 3) {
// ::std::map<BoundaryType, PeriodicBC*> &periodicBCs =
// mesh->getPeriodicBCMap();
// ::std::map<BoundaryType, PeriodicBC*>::iterator it;
// ::std::map<BoundaryType, PeriodicBC*>::iterator end = periodicBCs.end();
// for(it = periodicBCs.begin(); it != end; ++it) {
// it->second->updateAssociatedDOFs();
// }
// }
n_elements = mesh->getNumberOfLeaves();
spreadCoarsenMark();
......
......@@ -8,6 +8,8 @@
#include "MatrixVector.h"
#include "SystemVector.h"
#include "Estimator.h"
#include "ResidualEstimator.h"
#include "ResidualParallelEstimator.h"
#include "RecoveryEstimator.h"
#include "LeafData.h"
#include "SurfaceRegion_ED.h"
......@@ -212,6 +214,9 @@ namespace AMDiS {
creator = NEW ResidualEstimator::Creator;
addCreator("residual", creator);
creator = NEW ResidualParallelEstimator::Creator;
addCreator("residualparallel", creator);
creator = NEW RecoveryEstimator::Creator;
addCreator("recovery", creator);
}
......
......@@ -23,7 +23,6 @@ namespace AMDiS {
{
rowFESpace = NULL;
colFESpace = NULL;
elementMatrix = NULL;
}
DOFMatrix::DOFMatrix(const FiniteElemSpace* rowFESpace_,
......@@ -32,7 +31,6 @@ namespace AMDiS {
: rowFESpace(rowFESpace_),
colFESpace(colFESpace_),
name(name_),
elementMatrix(NULL),
coupleMatrix(false)
{
TEST_EXIT(rowFESpace)("no rowFESpace\n");
......@@ -50,22 +48,18 @@ namespace AMDiS {
DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
: name(rhs.name+"copy")
{
*this=rhs;
elementMatrix = NULL;
*this = rhs;
if (rowFESpace && rowFESpace->getAdmin())
(const_cast<DOFAdmin*>( rowFESpace->getAdmin()))->addDOFIndexed(this);
}
DOFMatrix::~DOFMatrix()
{
FUNCNAME("DOFMatrix::~DOFMatrix");
FUNCNAME("DOFMatrix::~DOFMatrix()");
if(rowFESpace && rowFESpace->getAdmin()) {
(const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
}
if(elementMatrix != NULL) {
DELETE elementMatrix;
}
//DELETE boundaryManager;
}
void DOFMatrix::print() const
......@@ -245,23 +239,45 @@ namespace AMDiS {
for (int i = 0; i < n_row; i++) { // for all rows of element matrix
row = elMat.rowIndices[i];
// if (omp_get_thread_num() == 0) {
// ::std::cout << "bound[i] = " << bound[i] << ::std::endl;
// }
BoundaryCondition *condition =
bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;
if (condition && condition->isDirichlet()) {
// if (omp_get_thread_num() == 0) {
// ::std::cout << "IN BC" << ::std::endl;
// }
MatrixRow *matrixRow = &(matrix[row]);
if (coupleMatrix) {
// if (omp_get_thread_num() == 0) {
// ::std::cout << "RESIZE 0" << ::std::endl;
// }
matrixRow->resize(0);
} else {
// if (omp_get_thread_num() == 0) {
// ::std::cout << "RESIZE 1 with col = " << row << ::std::endl;
// }
matrixRow->resize(1);
((*matrixRow)[0]).col = row;
((*matrixRow)[0]).entry = 1.0;
}
} else {
} else {
// if (omp_get_thread_num() == 0) {
// ::std::cout << "OUT BC" << ::std::endl;
// }
for (int j = 0; j < n_col; j++) { // for all columns
col = elMat.colIndices[j];
entry = elMat[i][j];
addSparseDOFEntry(sign, row, col, entry, add);
// if (omp_get_thread_num() == 0) {
// ::std::cout << "addSparseDOFEntry(" << sign << "," << row << "," << col << "," << entry << "," << add << ")" << ::std::endl;
// }
}
}
}
......@@ -308,11 +324,12 @@ namespace AMDiS {
int jcol, double entry,
bool add)
{
FUNCNAME("add_sparse_dof_entry");
FUNCNAME("DOFMatrix::addSparseDOFEntry()");
MatrixRow *row = &(matrix[irow]);
if(add && !entry) return NULL;
if (add && !entry)
return NULL;
double *result = NULL;
......@@ -328,6 +345,9 @@ namespace AMDiS {
MatEntry newEntry = {irow, 0.0};
row->push_back(newEntry);
rowSize = 1;
// if (omp_get_thread_num() == 0) {
// ::std::cout << "PB: " << irow << " " << 0.0 << ::std::endl;
// }
}
// search jcol
......@@ -355,16 +375,25 @@ namespace AMDiS {
if(!add) (*row)[i].entry = 0.0;
(*row)[i].entry += sign * entry;
result = &((*row)[i].entry);
if (omp_get_thread_num() == 0) {
// ::std::cout << "ADD: " << i << " " << sign * entry << ::std::endl;
}
} else {
if(freeCol == -1) {
MatEntry newEntry = {jcol, sign * entry};
row->push_back(newEntry);
result = &((*row)[row->size() - 1].entry);
if (omp_get_thread_num() == 0) {
// ::std::cout << "PB: " << jcol << " " << sign * entry << ::std::endl;
}
} else {
(*row)[freeCol].col = jcol;
if(!add) (*row)[freeCol].entry = 0.0;
(*row)[freeCol].entry += sign * entry;
result = &((*row)[freeCol].entry);
if (omp_get_thread_num() == 0) {
// ::std::cout << "ADD: " << jcol << " " << sign * entry << ::std::endl;
}
}
}
......@@ -446,19 +475,20 @@ namespace AMDiS {
}
ElementMatrix *DOFMatrix::assemble(double factor, ElInfo *elInfo,
const BoundaryType *bound, Operator *op)
void DOFMatrix::assemble(double factor, ElInfo *elInfo,
const BoundaryType *bound, Operator *op)
{
FUNCNAME("DOFMatrix::assemble()");
if (!op && operators.size() == 0) {
return NULL;
return;
}
Operator *operat = op ? op : operators[0];
elementMatrix =
operat->getAssembler()->initElementMatrix(elementMatrix, elInfo);
// Do no combine the next two lines into one!
ElementMatrix *elementMatrix = NULL;
elementMatrix = operat->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
if (op) {
op->getElementMatrix(elInfo, elementMatrix);
......