Commit f7866e7c authored by Dave Moxey's avatar Dave Moxey

Merge branch 'master' into 'fix/linearise_mapping'

# Conflicts:
#   utilities/NekMesh/ProcessModules/ProcessLinear.cpp
parents d3a92bd9 200fd1d0
......@@ -31,6 +31,9 @@ v4.4.0
file (!678)
- Extend ExtractDataToCoeffs to support interpolation between basis types for
quads and hexahedra (!682)
- Enabled MUMPS support in PETSc if a Fortran compiler was found and added 3D
support to the Helmholtz smoother used e.g. in FieldConverts C0Projection
module (!714)
**ADRSolver:**
- Add a projection equation system for C^0 projections (!675)
......@@ -48,6 +51,7 @@ v4.4.0
**FieldConvert:**
- Allow equi-spaced output for 1D and 2DH1D fields (!613)
- Update quality metric to include scaled Jacobian output (!695)
- Allow multiple XML files to be specified in InterpField module (!705)
**NekMesh:**
- Modify curve module to allow for spline input (!628)
......@@ -68,6 +72,9 @@ v4.4.0
- Add flag to `insertsurface` process for non-conforming geometries (!700)
- Bug fix to get two meshgen regression tests working (!700)
- Remove libANN in deference to boost::geometry (!703)
- 2D to 3D mesh extrusion module (!715)
- Add a mesh extract option to the linearise module to visualise the result
(!712)
**FieldConvert:**
- Move all modules to a new library, FieldUtils, to support post-processing
......@@ -82,6 +89,7 @@ v4.3.5
- Fix bug in DG with hybrid meshes (!694)
- Fix issue with parallel output (!699)
- Fix performance issue with iterative full solver (!693)
- Enforced precision on history point output (!706)
**Documentation**
- Update build instructions in user guide for Windows (!692)
......
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.7)
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
SET(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build,
options are: None(CMAKE_CXX_FLAGS or CMAKE_C_FLAGS used) Debug Release
RelWithDebInfo MinSizeRel.")
PROJECT(Nektar++)
PROJECT(Nektar++ C CXX)
INCLUDE(CheckLanguage)
CHECK_LANGUAGE(Fortran)
IF(CMAKE_Fortran_COMPILER)
ENABLE_LANGUAGE(Fortran)
ELSE()
MESSAGE(STATUS "No Fortran support")
ENDIF()
# Helps organize projects in IDEs.
SET_PROPERTY(GLOBAL PROPERTY USE_FOLDERS ON)
......
......@@ -30,16 +30,30 @@ IF (NEKTAR_USE_PETSC)
SET(PETSC_C_COMPILER "${CMAKE_C_COMPILER}")
SET(PETSC_CXX_COMPILER "${CMAKE_CXX_COMPILER}")
SET(PETSC_Fortran_COMPILER "${CMAKE_Fortran_COMPILER}")
IF (NEKTAR_USE_MPI)
IF (NOT MPI_BUILTIN)
SET(PETSC_C_COMPILER "${MPI_C_COMPILER}")
SET(PETSC_CXX_COMPILER "${MPI_CXX_COMPILER}")
SET(PETSC_Fortran_COMPILER "${MPI_Fortran_COMPILER}")
ENDIF (NOT MPI_BUILTIN)
ELSE (NEKTAR_USE_MPI)
SET(PETSC_NO_MPI "--with-mpi=0")
ENDIF (NEKTAR_USE_MPI)
IF(CMAKE_Fortran_COMPILER)
IF(NEKTAR_USE_MPI AND NOT MPI_Fortran_COMPILER)
MESSAGE(ERROR "MPI_Fortran_COMPILER not set")
ENDIF()
# we use a MUMPS build in ordering here, in the future it might make
# sense to hook it up with metis/scotch since this MIGHT be faster
SET(PETSC_MUMPS --download-scalapack --download-mumps)
ELSE()
MESSAGE(WARNING "No Fortran support. Building PETSc without MUMPS support")
SET(PETSC_Fortran_COMPILER "0")
ENDIF()
EXTERNALPROJECT_ADD(
petsc-3.7.2
PREFIX ${TPSRC}
......@@ -52,17 +66,21 @@ IF (NEKTAR_USE_PETSC)
URL http://www.nektar.info/thirdparty/petsc-lite-3.7.2.tar.gz
URL_MD5 "26c2ff8eaaa9e49aea063f839f5daa7e"
CONFIGURE_COMMAND
OMPI_FC=${CMAKE_Fortran_COMPILER}
OMPI_CC=${CMAKE_C_COMPILER}
OMPI_CXX=${CMAKE_CXX_COMPILER}
${PYTHON_EXECUTABLE} ./configure
./configure
--with-fc=${PETSC_Fortran_COMPILER}
--with-cc=${PETSC_C_COMPILER}
--with-cxx=${PETSC_CXX_COMPILER}
--with-shared-libraries=1
--with-pic=1
--with-x=0
--with-ssl=0
--prefix=${TPDIST}
--with-petsc-arch=c-opt
--with-fc=0
${PETSC_MUMPS}
${PETSC_NO_MPI}
BUILD_COMMAND MAKEFLAGS= make)
......
......@@ -632,6 +632,20 @@ match the interior domain:
This module should be added to the module chain for any complex if the user
suspects there may be a mesh issue. The module will print a warning if there is.
\subsection{2D mesh extrusion}
This module allows a 2D mesh, quads, triangles or both, to be extruded in the
$z$ direction to make a simple 3D mesh made of prisms and hexahedra. It is also
capable of extruding the high-order curvature within the 2D mesh. The module
requires two parameters:
\begin{lstlisting}[style=BashInputStyle]
NekMesh -m extrude:layers=n:length=l 2D.xml 3D.xml
\end{lstlisting}
length which determines how long the $z$ extrusion will be and layers, the
number of elements in the $z$ direction.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
......
......@@ -90,7 +90,7 @@ void ProcessInterpField::Process(po::variables_map &vm)
std::vector<std::string> files;
// set up session file for from field
files.push_back(m_config["fromxml"].as<string>());
ParseUtils::GenerateOrderedStringVector(m_config["fromxml"].as<string>().c_str(), files);
m_fromField->m_session =
LibUtilities::SessionReader::CreateInstance(0, 0, files);
......
......@@ -445,8 +445,8 @@ void ProcessInterpPoints::Process(po::variables_map &vm)
FieldSharedPtr fromField = boost::shared_ptr<Field>(new Field());
std::vector<std::string> files;
ParseUtils::GenerateOrderedStringVector(m_config["fromxml"].as<string>().c_str(), files);
// set up session file for from field
files.push_back(m_config["fromxml"].as<string>());
fromField->m_session =
LibUtilities::SessionReader::CreateInstance(0, 0, files);
......
......@@ -986,7 +986,7 @@ namespace Nektar
bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
}
}
m_locToGloMap->UniversalAssemble(wsp);
m_locToGloMap->UniversalAssemble(gamma);
// Add weak boundary conditions to forcing
Vmath::Vadd(contNcoeffs, wsp, 1, gamma, 1, wsp, 1);
......
......@@ -631,7 +631,83 @@ namespace Nektar
GeneralMatrixOp_IterPerExp(gkey,inarray,outarray);
}
}
/**
* First compute the inner product of forcing function with respect to
* base, and then solve the system with the linear advection operator.
* @param velocity Array of advection velocities in physical space
* @param inarray Forcing function.
* @param outarray Result.
* @param lambda reaction coefficient
* @param coeffstate State of Coefficients, Local or Global
* @param dirForcing Dirichlet Forcing.
*/
// could combine this with HelmholtzCG.
void ContField3D::v_LinearAdvectionDiffusionReactionSolve(
const Array<OneD, Array<OneD, NekDouble> > &velocity,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
const NekDouble lambda,
CoeffState coeffstate,
const Array<OneD, const NekDouble> &dirForcing)
{
// Inner product of forcing
int contNcoeffs = m_locToGloMap->GetNumGlobalCoeffs();
Array<OneD, NekDouble> wsp(contNcoeffs);
IProductWRTBase(inarray, wsp, eGlobal);
// Note -1.0 term necessary to invert forcing function to
// be consistent with matrix definition
Vmath::Neg(contNcoeffs, wsp, 1);
// Forcing function with weak boundary conditions
int i, j;
int bndcnt = 0;
Array<OneD, NekDouble> gamma(contNcoeffs, 0.0);
for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
{
if (m_bndConditions[i]->GetBoundaryConditionType() !=
SpatialDomains::eDirichlet)
{
for (j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); j++)
{
gamma[m_locToGloMap->GetBndCondCoeffsToGlobalCoeffsMap(
bndcnt++)] += (m_bndCondExpansions[i]->GetCoeffs())[j];
}
}
else
{
bndcnt += m_bndCondExpansions[i]->GetNcoeffs();
}
}
m_locToGloMap->UniversalAssemble(gamma);
// Add weak boundary conditions to forcing
Vmath::Vadd(contNcoeffs, wsp, 1, gamma, 1, wsp, 1);
// Solve the system
StdRegions::ConstFactorMap factors;
factors[StdRegions::eFactorLambda] = lambda;
StdRegions::VarCoeffMap varcoeffs;
varcoeffs[StdRegions::eVarCoeffVelX] = velocity[0];
varcoeffs[StdRegions::eVarCoeffVelY] = velocity[1];
varcoeffs[StdRegions::eVarCoeffVelZ] = velocity[2];
GlobalLinSysKey key(StdRegions::eLinearAdvectionDiffusionReaction,
m_locToGloMap,
factors,
varcoeffs);
if (coeffstate == eGlobal)
{
GlobalSolve(key, wsp, outarray, dirForcing);
}
else
{
Array<OneD, NekDouble> tmp(contNcoeffs, 0.0);
GlobalSolve(key, wsp, tmp, dirForcing);
GlobalToLocal(tmp, outarray);
}
}
int ContField3D::GetGlobalMatrixNnz(const GlobalMatrixKey &gkey)
{
ASSERTL1(gkey.LocToGloMapIsDefined(),
......
......@@ -186,6 +186,16 @@ namespace Nektar
const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate);
// Solve the linear advection problem assuming that m_coeffs
// vector contains an intial estimate for solution
MULTI_REGIONS_EXPORT virtual void v_LinearAdvectionDiffusionReactionSolve(
const Array<OneD, Array<OneD, NekDouble> > &velocity,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
const NekDouble lambda,
CoeffState coeffstate = eLocal,
const Array<OneD, const NekDouble> &dirForcing = NullNekDouble1DArray);
virtual void v_ClearGlobalLinSysManager(void);
......
......@@ -81,6 +81,17 @@ namespace Nektar
PetscInitialized(&isInitialized);
if (!isInitialized)
{
#ifdef NEKTAR_USE_MPI
std::string commType =
m_expList.lock()->GetSession()->GetComm()->GetType();
if (commType.find("MPI") != std::string::npos)
{
LibUtilities::CommMpiSharedPtr comm =
boost::static_pointer_cast<LibUtilities::CommMpi>(
m_expList.lock()->GetSession()->GetComm());
PETSC_COMM_WORLD = comm->GetComm();
}
#endif
PetscInitializeNoArguments();
}
......@@ -159,6 +170,13 @@ namespace Nektar
// Do system solve
KSPSolve(m_ksp, m_b, m_x);
KSPConvergedReason reason;
KSPGetConvergedReason(m_ksp, &reason);
ASSERTL0(reason > 0,
"PETSc solver diverged, reason is: " +
std::string(KSPConvergedReasons[reason]));
// Scatter results to local vector
VecScatterBegin(m_ctx, m_x, m_locVec,
INSERT_VALUES, SCATTER_FORWARD);
......
......@@ -402,9 +402,9 @@ void FilterHistoryPoints::v_Initialise(
gloCoord[2]);
m_outputStream << "# " << boost::format("%6.0f") % i;
m_outputStream << " " << boost::format("%25e") % gloCoord[0];
m_outputStream << " " << boost::format("%25e") % gloCoord[1];
m_outputStream << " " << boost::format("%25e") % gloCoord[2];
m_outputStream << " " << boost::format("%25.19e") % gloCoord[0];
m_outputStream << " " << boost::format("%25.19e") % gloCoord[1];
m_outputStream << " " << boost::format("%25.19e") % gloCoord[2];
m_outputStream << endl;
}
......@@ -575,10 +575,10 @@ void FilterHistoryPoints::v_Update(const Array<OneD, const MultiRegions::ExpList
// Write data values point by point
for (k = 0; k < m_historyPoints.size(); ++k)
{
m_outputStream << boost::format("%25e") % time;
m_outputStream << boost::format("%25.19e") % time;
for (int j = 0; j < numFields; ++j)
{
m_outputStream << " " << boost::format("%25e") % data[k*numFields+j];
m_outputStream << " " << boost::format("%25.19e") % data[k*numFields+j];
}
m_outputStream << endl;
}
......
......@@ -897,7 +897,7 @@ namespace Nektar
v_BwdTrans(inarray,tmp);
VarCoeffType varcoefftypes[] = {eVarCoeffVelX, eVarCoeffVelY};
VarCoeffType varcoefftypes[] = {eVarCoeffVelX, eVarCoeffVelY, eVarCoeffVelZ};
//calculate u dx + v dy + ..
Vmath::Zero(totpts,tmp_adv,1);
......
......@@ -207,7 +207,8 @@ namespace Nektar
eVarCoeffD02,
eVarCoeffD12,
eVarCoeffVelX,
eVarCoeffVelY
eVarCoeffVelY,
eVarCoeffVelZ
};
const char* const VarCoeffTypeMap[] = {
......
......@@ -27,6 +27,7 @@ SET(NekMeshHeaders
ProcessModules/ProcessTetSplit.h
ProcessModules/ProcessOptiExtract.h
ProcessModules/ProcessInsertSurface.h
ProcessModules/ProcessExtrude.h
)
SET(NekMeshSources
......@@ -59,6 +60,7 @@ SET(NekMeshSources
ProcessModules/ProcessTetSplit.cpp
ProcessModules/ProcessOptiExtract.cpp
ProcessModules/ProcessInsertSurface.cpp
ProcessModules/ProcessExtrude.cpp
)
IF (NEKTAR_USE_CCM)
......
///////////////////////////////////////////////////////////////////////////////
//
// File: ProcessExtrude.cpp
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Extrude a two-dimensional mesh to a three-dimensional mesh.
//
////////////////////////////////////////////////////////////////////////////////
#include <NekMeshUtils/MeshElements/Element.h>
#include "ProcessExtrude.h"
using namespace std;
using namespace Nektar::NekMeshUtils;
namespace Nektar
{
namespace Utilities
{
ModuleKey ProcessExtrude::className =
GetModuleFactory().RegisterCreatorFunction(
ModuleKey(eProcessModule, "extrude"), ProcessExtrude::create);
ProcessExtrude::ProcessExtrude(MeshSharedPtr m) : ProcessModule(m)
{
m_config["layers"] =
ConfigOption(false, "5", "Number of layers to extrude");
m_config["length"] = ConfigOption(false, "1.0", "Length of extrusion");
}
ProcessExtrude::~ProcessExtrude()
{
}
void ProcessExtrude::Process()
{
int nLayers = m_config["layers"].as<int>();
NekDouble length = m_config["length"].as<NekDouble>();
NekDouble dz = length / nLayers;
ASSERTL0(m_mesh->m_spaceDim == 2,
"Extrude should only be called for a two dimensional mesh");
// Increment space and expansion dimensions.
m_mesh->m_spaceDim++;
m_mesh->m_expDim++;
// Grab a copy of the existing two-dimensional elements.
vector<ElementSharedPtr> el = m_mesh->m_element[2];
// Reset mesh.
for (int d = 0; d <= 3; ++d)
{
m_mesh->m_element[d].clear();
}
NodeSet nodes = m_mesh->m_vertexSet;
map<int, NodeSharedPtr> id2node;
NodeSet::iterator it;
for (it = nodes.begin(); it != nodes.end(); ++it)
{
id2node[(*it)->m_id] = *it;
}
// Create vertices for subsequent layers.
for (int i = 1; i < nLayers + 1; ++i)
{
for (it = nodes.begin(); it != nodes.end(); ++it)
{
NodeSharedPtr n = *it;
NodeSharedPtr newNode(
new Node(i * nodes.size() + n->m_id, n->m_x, n->m_y, i * dz));
m_mesh->m_vertexSet.insert(newNode);
id2node[i * nodes.size() + n->m_id] = newNode;
}
}
EdgeSet es = m_mesh->m_edgeSet; // copy edges for curvature
for (int j = 0; j < nLayers; ++j)
{
for (int i = 0; i < el.size(); ++i)
{
vector<NodeSharedPtr> verts = el[i]->GetVertexList();
if (verts.size() == 4)
{
vector<NodeSharedPtr> nodeList(8);
nodeList[0] = id2node[verts[0]->m_id + j * nodes.size()];
nodeList[1] = id2node[verts[1]->m_id + j * nodes.size()];
nodeList[2] = id2node[verts[1]->m_id + (j + 1) * nodes.size()];
nodeList[3] = id2node[verts[0]->m_id + (j + 1) * nodes.size()];
nodeList[4] = id2node[verts[3]->m_id + j * nodes.size()];
nodeList[5] = id2node[verts[2]->m_id + j * nodes.size()];
nodeList[6] = id2node[verts[2]->m_id + (j + 1) * nodes.size()];
nodeList[7] = id2node[verts[3]->m_id + (j + 1) * nodes.size()];
vector<int> tags(1);
tags[0] = 0;
ElmtConfig conf(LibUtilities::eHexahedron, 1, false, false);
ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::eHexahedron, conf, nodeList, tags);
m_mesh->m_element[3].push_back(E);
}
else
{
vector<NodeSharedPtr> nodeList(6);
nodeList[0] = id2node[verts[0]->m_id + (j + 1) * nodes.size()];
nodeList[1] = id2node[verts[1]->m_id + (j + 1) * nodes.size()];
nodeList[2] = id2node[verts[1]->m_id + j * nodes.size()];
nodeList[3] = id2node[verts[0]->m_id + j * nodes.size()];
nodeList[4] = id2node[verts[2]->m_id + (j + 1) * nodes.size()];
nodeList[5] = id2node[verts[2]->m_id + j * nodes.size()];
vector<int> tags(1);
tags[0] = 1;
ElmtConfig conf(LibUtilities::ePrism, 1, false, false);
ElementSharedPtr E = GetElementFactory().CreateInstance(
LibUtilities::ePrism, conf, nodeList, tags);
m_mesh->m_element[3].push_back(E);
}
}
}
ProcessEdges();
ProcessFaces();
ProcessElements();
ProcessComposites();
EdgeSet::iterator eit;
for (eit = es.begin(); eit != es.end(); eit++)
{
if ((*eit)->m_edgeNodes.size() > 0)
{
for (int j = 0; j < nLayers + 1; ++j)
{
vector<NodeSharedPtr> ns((*eit)->m_edgeNodes.size());
for (int i = 0; i < ns.size(); i++)
{
NodeSharedPtr n = (*eit)->m_edgeNodes[i];
ns[i] = boost::shared_ptr<Node>(
new Node(0, n->m_x, n->m_y, j * dz));
}
EdgeSharedPtr e = boost::shared_ptr<Edge>(
new Edge(id2node[(*eit)->m_n1->m_id + j * nodes.size()],
id2node[(*eit)->m_n2->m_id + j * nodes.size()]));
EdgeSet::iterator f = m_mesh->m_edgeSet.find(e);
ASSERTL0(f != m_mesh->m_edgeSet.end(), "could not find edge");
if ((*f)->m_n1 == e->m_n1)
{
(*f)->m_edgeNodes = ns;
}
else
{
reverse(ns.begin(), ns.end());
(*f)->m_edgeNodes = ns;
}
}
}
}
}
}
}
////////////////////////////////////////////////////////////////////////////////
//
// File: ProcessExtrude.h
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: extrude a 2d mesh.
//
////////////////////////////////////////////////////////////////////////////////
#ifndef UTILITIES_NEKMESH_PROCESSEXTRUDE
#define UTILITIES_NEKMESH_PROCESSEXTRUDE
#include "../Module.h"
namespace Nektar
{
namespace Utilities
{
/**
* @brief This processing module extrudes a 2d mesh in the z direction
*/
class ProcessExtrude : public ProcessModule
{
public:
/// Creates an instance of this class
static boost::shared_ptr<Module> create(MeshSharedPtr m)
{
return MemoryManager<ProcessExtrude>::AllocateSharedPtr(m);
}
static ModuleKey className;
ProcessExtrude(MeshSharedPtr m);
virtual ~ProcessExtrude();
virtual void Process();
};
}
}
#endif
......@@ -55,6 +55,8 @@ ProcessLinear::ProcessLinear(MeshSharedPtr m) : ProcessModule(m)
ConfigOption(false, "0", "remove curve nodes if element is invalid.");
m_config["prismonly"] =
ConfigOption(true, "0", "only acts on prims");
m_config["extract"] =
ConfigOption(false, "", "dump a mesh of the extracted elements");
}
ProcessLinear::~ProcessLinear()
......@@ -138,6 +140,8 @@ void ProcessLinear::Process()
vector<NodeSharedPtr> zeroNodes;
boost::unordered_set<int> clearedEdges, clearedFaces, clearedElmts;
vector<ElementSharedPtr> dumpEls;
// Iterate over list of elements of expansion dimension.
while(el.size() > 0)
{
......@@ -153,6 +157,7 @@ void ProcessLinear::Process()
if (Invalid(el[i],thr))//(!gfac->IsValid())
{
dumpEls.push_back(el[i]);
clearedElmts.insert(el[i]->GetId());;
el[i]->SetVolumeNodes(zeroNodes);
......@@ -216,6 +221,26 @@ void ProcessLinear::Process()
<< " elements (" << clearedEdges.size() << " edges, "
<< clearedFaces.size() << " faces)" << endl;
}
if(m_config["extract"].beenSet)
{
MeshSharedPtr dmp = boost::shared_ptr<Mesh>(new Mesh());
dmp->m_expDim = 3;
dmp->m_spaceDim = 3;
dmp->m_nummode = 2;
dmp->m_element[3] = dumpEls;
ModuleSharedPtr mod = GetModuleFactory().CreateInstance(
ModuleKey(eOutputModule, "xml"), dmp);
mod->RegisterConfig("outfile", m_config["extract"].as<string>().c_str());
mod->ProcessVertices();
mod->ProcessEdges();
mod->ProcessFaces();
mod->ProcessElements();
mod->ProcessComposites();
mod->Process();
}
}
}
......
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