1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580
|
#ifndef FE_ENGINE_H
#define FE_ENGINE_H
#include <vector>
#include <map>
#include <set>
#include <utility>
#include <string>
#include "ATC_TypeDefs.h"
#include "Array.h"
#include "Array2D.h"
#include "FE_Mesh.h"
#include "PhysicsModel.h"
#include "OutputManager.h"
#include "MeshReader.h"
#include "mpi.h"
namespace ATC {
class ATC_Method;
class FE_Element;
class XT_Function;
class KernelFunction;
/**
* @class FE_Engine
* @brief Base class for computing and assembling mass matrix
* and rhs vectors
*/
class FE_Engine{
public:
/** constructor/s */
FE_Engine(MPI_Comm comm);
/** destructor */
~FE_Engine();
/** initialize */
void initialize();
MPI_Comm communicator() {return communicator_;}
void partition_mesh();
void departition_mesh();
bool is_partitioned() const { return feMesh_->is_partitioned(); }
int map_elem_to_myElem(int elemID) const
{ return feMesh_->map_elem_to_myElem(elemID); }
int map_myElem_to_elem(int myElemID) const
{ return feMesh_->map_myElem_to_elem(myElemID); }
// note: it is misleading to declare the following const
// because it touches the nIPsPer* data members, which
// are now declared mutable. Why? Well, set_quadrature
// has to be called from a const function, and all the
// matrices dependent on nIPsPer* are declared mutable
// as well (and have been). I think this is because a
// const engine needs to be able to deal with various
// quadratures and update its data members directly, which
// are really convenience-copies of data members that
// are more pertinent to other classes (FE_Interpolate,
// for the most part) that it uses temporarily for space/
// time speedups while doing it's computations.
//
// I approve of this usage of mutable, but the const/
// non-const member function declaring in this class is
// really all wrong to begin with.
/** set quadrature scheme, resize matrices if necessary as per
* initialize() */
void set_quadrature(FeIntQuadrature quadType, bool temp=true) const;
/** parser/modifier */
bool modify(int narg, char **arg);
/** finish up */
void finish();
/** print out the "global connectivity" of all elements */
void print_mesh() const;
//----------------------------------------------------------------
/** \name output */
//----------------------------------------------------------------
/*@{*/
/** these assume the caller is handling the parallel collection */
void initialize_output(int rank, std::string outputPrefix, std::set<int> otypes);
/** write geometry */
void write_geometry(void);
/** write data: data is arrayed over _unique_ nodes
and then mapped by the engine */
void write_data(double time, FIELDS &soln, OUTPUT_LIST *data=nullptr);
void write_data(double time, OUTPUT_LIST *data);
void write_restart_file(std::string fileName, RESTART_LIST *data)
{ outputManager_.write_restart_file(fileName,data); }
void read_restart_file(std::string fileName, RESTART_LIST *data)
{ outputManager_.read_restart_file(fileName,data); }
void delete_elements(const std::set<int> &elementList);
void cut_mesh(const std::set<PAIR> &cutFaces, const std::set<int> &edgeNodes);
void add_global(const std::string name, const double value)
{ outputManager_.add_global(name,value); }
void add_field_names(const std::string field, const std::vector<std::string> & names)
{ outputManager_.add_field_names(field,names); }
void reset_globals() { outputManager_.reset_globals(); }
/** pass through to access output manager */
OutputManager *output_manager() { return &outputManager_; }
/*@}*/
//----------------------------------------------------------------
/** \name assembled matrices and vectors */
//----------------------------------------------------------------
/*@{*/
DENS_VEC interpolate_field(const DENS_VEC & x, const FIELD & f) const;
/** interpolate fields */
void interpolate_fields(const int ielem,
const FIELDS &fields,
AliasArray<int> &conn,
DENS_MAT &N,
DIAG_MAT &weights,
std::map<FieldName,DENS_MAT> &fieldsAtIPs) const;
/** interpolate fields & gradients */
void interpolate_fields(const int ielem,
const FIELDS &fields,
AliasArray<int> &conn,
DENS_MAT &N,
DENS_MAT_VEC &dN,
DIAG_MAT &weights,
FIELD_MATS &fieldsAtIPs,
GRAD_FIELD_MATS &grad_fieldsAtIPs) const;
/** compute a dimensionless stiffness matrix */
void stiffness_matrix(SPAR_MAT &matrix) const;
/** compute tangent matrix for a pair of fields - native quadrature */
void compute_tangent_matrix(
const RHS_MASK &rhsMask,
const std::pair<FieldName,FieldName> row_col,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<int> &elementMaterials,
SPAR_MAT &tangent,
const DenseMatrix<bool> *elementMask=nullptr) const;
/** compute tangent matrix for a pair of fields - given quadrature */
void compute_tangent_matrix(const RHS_MASK &rhsMask,
const std::pair<FieldName,FieldName> row_col,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<std::set<int> > &pointMaterialGroups,
const DIAG_MAT &weights,
const SPAR_MAT &N,
const SPAR_MAT_VEC &dN,
SPAR_MAT &tangent,
const DenseMatrix<bool> *elementMask=nullptr) const;
/** compute a consistent mass matrix for a field */
void compute_mass_matrix(
const Array<FieldName> &mask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<int> &elementMaterials,
CON_MASS_MATS &mass_matrix,
const DenseMatrix<bool> *elementMask=nullptr) const;
/** compute a dimensionless mass matrix */
void compute_mass_matrix(SPAR_MAT &mass_matrix) const;
/** computes a dimensionless mass matrix for the given-quadrature */
void compute_mass_matrix(const DIAG_MAT &weights,
const SPAR_MAT &N,
SPAR_MAT &mass_matrix) const;
/** compute a single dimensionless mass matrix */
void compute_lumped_mass_matrix(DIAG_MAT &lumped_mass_matrix) const;
/** compute lumped mass matrix = diag (\int \rho N_I dV) */
void compute_lumped_mass_matrix(
const Array<FieldName> &mask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<int> &elementMaterials,
MASS_MATS &mass_matrix,
const DenseMatrix<bool> *elementMask=nullptr) const;
/** compute dimensional lumped mass matrix using given quadrature */
void compute_lumped_mass_matrix(
const Array<FieldName> &mask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<std::set<int> > &pointMaterialGroups,
const DIAG_MAT &weights,
const SPAR_MAT &N,
MASS_MATS &mass_matrix) const;
/** compute an approximation to a finite difference gradient from mesh */
void compute_gradient_matrix(SPAR_MAT_VEC &grad_matrix) const;
/** compute energy */
void compute_energy(const Array<FieldName> &mask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<int> &elementMaterials,
FIELD_MATS &energy,
const DenseMatrix<bool> *elementMask=nullptr,
const IntegrationDomainType domain=FULL_DOMAIN) const;
/** compute residual or RHS of the dynamic weak eqn */
void compute_rhs_vector(
const RHS_MASK &rhsMask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<int> &elementMaterials,
FIELDS &rhs,
bool freeOnly=false,
const DenseMatrix<bool> *elementMask=nullptr) const;
/** compute RHS for given quadrature */
void compute_rhs_vector(const RHS_MASK &rhsMask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<std::set<int> > &pointMaterialGroups,
const DIAG_MAT &weights,
const SPAR_MAT &N,
const SPAR_MAT_VEC &dN,
FIELDS &rhs) const;
/** compute pointwise source for given quadrature */
void compute_source(const Array2D<bool> &rhsMask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<std::set<int> > &pointMaterialGroups,
const DIAG_MAT &weights,
const SPAR_MAT &N,
const SPAR_MAT_VEC &dN,
FIELD_MATS &sources) const;
/** compute flux in domain i.e. N^T B_integrand */
void compute_flux(const RHS_MASK &rhsMask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<int> &elementMaterials,
GRAD_FIELD_MATS &flux,
const DenseMatrix<bool> *elementMask=nullptr) const;
/** compute the flux on the MD/FE boundary */
void compute_boundary_flux(const RHS_MASK &rhsMask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<int> &elementMaterials,
const std::set<PAIR> &faceSet,
FIELDS &rhs) const;
/** compute the flux on using an L2 interpolation of the flux */
void compute_boundary_flux(const RHS_MASK &rhsMask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<int> &elementMaterials,
const Array<std::set<int> > &pointMaterialGroups,
const DIAG_MAT &weights,
const SPAR_MAT &N,
const SPAR_MAT_VEC &dN,
const DIAG_MAT &flux_mask,
FIELDS &rhs,
const DenseMatrix<bool> *elementMask=nullptr,
const std::set<int> *nodeSet=nullptr) const;
/** compute prescribed flux given an array of functions of x & t */
void add_fluxes(const Array<bool> &fieldMask,
const double time,
const SURFACE_SOURCE &sourceFunctions,
FIELDS &nodalSources) const;
void compute_fluxes(const Array<bool> &fieldMask,
const double time,
const SURFACE_SOURCE &sourceFunctions,
FIELDS &nodalSources) const
{
SURFACE_SOURCE::const_iterator src_iter;
for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++) {
_fieldName_ = src_iter->first;
if (!fieldMask((int)_fieldName_)) continue;
if (nodalSources[_fieldName_].nRows()==0) {
nodalSources[_fieldName_].reset(nNodesUnique_,1);
}
}
add_fluxes(fieldMask, time, sourceFunctions, nodalSources);
}
/** compute robin flux given an array of functions of u, x & t */
void add_robin_fluxes(const Array2D<bool> &rhsMask,
const FIELDS &fields,
const double time,
const ROBIN_SURFACE_SOURCE &sourceFunctions,
FIELDS &nodalSources) const;
void add_robin_tangent(const Array2D<bool> &rhsMask,
const FIELDS &fields,
const double time,
const ROBIN_SURFACE_SOURCE &sourceFunctions,
SPAR_MAT &tangent) const;
/** compute open flux given a face set */
void add_open_fluxes(const Array2D<bool> &rhsMask,
const FIELDS &fields,
const OPEN_SURFACE &openFaces,
FIELDS &nodalSources,
const FieldName velocity = ELECTRON_VELOCITY) const;
void add_open_tangent(const Array2D<bool> &rhsMask,
const FIELDS &fields,
const OPEN_SURFACE &openFaces,
SPAR_MAT &tangent,
const FieldName velocity = ELECTRON_VELOCITY) const;
/** compute nodal vector of volume based sources */
void add_sources(const Array<bool> &fieldMask,
const double time,
const VOLUME_SOURCE &sourceFunctions,
FIELDS &nodalSources) const;
void add_sources(const Array<bool> &fieldMask,
const double time,
const FIELDS &sources,
FIELDS &nodalSources) const;
/** compute surface flux of a nodal field */
void field_surface_flux(const DENS_MAT &field,
const std::set<PAIR> &faceSet,
DENS_MAT &values,
const bool contour=false,
const int axis=2) const;
/** integrate a nodal field over an element set */
DENS_VEC integrate(const DENS_MAT &field, const ESET & eset) const;
/** integrate a nodal field over an face set */
DENS_VEC integrate(const DENS_MAT & /* field */, const FSET & /* fset */) const
{ throw ATC_Error(FILELINE,"unimplemented function"); }
/*@}*/
//----------------------------------------------------------------
/** \name shape functions */
//----------------------------------------------------------------
/*@{*/
/** evaluate shape function at a list of points in R^3 */
void evaluate_shape_functions(const MATRIX &coords,
SPAR_MAT &N) const;
/** evaluate shape function & derivatives at a list of points in R^3 */
void evaluate_shape_functions(const MATRIX &coords,
SPAR_MAT &N,
SPAR_MAT_VEC &dN) const;
/** evaluate shape function at a list of points in R^3 */
void evaluate_shape_functions(const MATRIX &coords,
const INT_ARRAY &pointToEltMap,
SPAR_MAT &N) const;
/** evaluate shape function & derivatives at a list of points in R^3 */
void evaluate_shape_functions(const MATRIX &coords,
const INT_ARRAY &pointToEltMap,
SPAR_MAT &N,
SPAR_MAT_VEC &dN) const;
/** evaluate shape derivatives at a list of points in R^3 */
void evaluate_shape_function_derivatives(const MATRIX &coords,
const INT_ARRAY &pointToEltMap,
SPAR_MAT_VEC &dN) const;
void shape_functions(const VECTOR &x,
DENS_VEC &shp,
Array<int> &node_list) const
{ feMesh_->shape_functions(x,shp,node_list); }
void shape_functions(const VECTOR & x,
DENS_VEC& shp,
DENS_MAT& dshp,
Array<int> &node_list) const
{ feMesh_->shape_functions(x,shp,dshp,node_list); }
void shape_functions(const VECTOR &x,
const int eltId,
DENS_VEC& shp,
Array<int> &node_list) const
{ feMesh_->shape_functions(x,eltId,shp,node_list); }
void shape_functions(const VECTOR &x,
DENS_VEC& shp,
Array<int> &node_list,
int &eltId) const
{ feMesh_->shape_functions(x,shp,node_list,eltId); }
void shape_functions(const VECTOR &x,
const int eltId,
DENS_VEC &shp,
DENS_MAT &dshp,
Array<int> &node_list) const
{ feMesh_->shape_functions(x,eltId,shp,dshp,node_list); }
/*@}*/
//----------------------------------------------------------------
/** \name kernel functions */
//----------------------------------------------------------------
/** evaluate kernel function */
void evaluate_kernel_functions(const MATRIX &pt_coords,
SPAR_MAT &N) const;
/** kernel matrix bandwidth */
int kernel_matrix_bandwidth(const MATRIX &pt_coords) const;
//----------------------------------------------------------------
/** \name nodeset */
//----------------------------------------------------------------
/** pass through */
void create_nodeset(const std::string &name, const std::set<int> &nodeset)
{ feMesh_->create_nodeset(name,nodeset); }
//----------------------------------------------------------------
/** \name accessors */
//----------------------------------------------------------------
/*@{*/
/** even though these are pass-throughs there is a necessary
* translation */
/** return number of unique nodes */
int num_nodes() const { return feMesh_->num_nodes_unique(); }
/** return number of total nodes */
int nNodesTotal() const { return feMesh_->num_nodes(); }
/** return number of elements */
int num_elements() const { return feMesh_->num_elements(); }
int my_num_elements() const { return feMesh_->my_num_elements(); }
/** return number of nodes per element */
int num_nodes_per_element() const { return feMesh_->num_nodes_per_element(); }
/** return element connectivity */
void element_connectivity(const int eltID,
Array<int> & nodes) const
{ feMesh_->element_connectivity_unique(eltID, nodes); }
/** return face connectivity */
void face_connectivity(const PAIR &faceID,
Array<int> &nodes) const
{ feMesh_->face_connectivity_unique(faceID, nodes); }
/** in lieu of pass-throughs const accessors ... */
/** return const ptr to mesh */
const FE_Mesh* fe_mesh() const { return feMesh_; }
/** return number of spatial dimensions */
int nsd() const { return feMesh_->num_spatial_dimensions(); }
/** return if the FE mesh has been created */
int has_mesh() const { return feMesh_!=nullptr; }
/** get nodal coordinates for a given element */
void element_coordinates(const int eltIdx, DENS_MAT &coords)
{ feMesh_->element_coordinates(eltIdx,coords); }
/** get nodal coordinates for a given element */
void element_field(const int eltIdx, const DENS_MAT field,
DENS_MAT &local_field)
{ feMesh_->element_field(eltIdx, field, local_field); }
/** access list of elements to be deleted */
const std::set<int> &null_elements(void) const
{ return nullElements_; }
/** access to the amended nodal coordinate values */
const DENS_MAT &nodal_coordinates(void) const
{ return (*feMesh_->coordinates()); }
/** map global node numbering to unique node numbering for
* amended mesh */
int map_global_to_unique(const int global_id) const
{ return (*feMesh_->node_map())(global_id); }
int number_of_global_nodes(void) const { return nNodes_; }
/*@}*/
/** set kernel */
void set_kernel(KernelFunction* ptr);
KernelFunction *kernel(int /* i */) { return kernelFunction_; }
KernelFunction *kernel() { return kernelFunction_; }
private:
//----------------------------------------------------------------
/** mesh setup commands (called from modify) */
//----------------------------------------------------------------
/*@{*/
MPI_Comm communicator_;
/** finite element mesh */
FE_Mesh *feMesh_;
/** auxiliary kernel function */
KernelFunction *kernelFunction_;
/** initialized flag */
bool initialized_;
void parse_partitions(int & argIdx, int narg, char ** arg,
int idof, Array<double> & dx ) const;
void print_partitions(double xmin, double xmax, Array<double> &dx) const;
/** create a uniform, structured mesh */
void create_mesh(Array<double> &dx,
Array<double> &dy,
Array<double> &dz,
const char *regionName,
Array<bool> periodic);
void create_mesh(int nx, int ny, int nz,
const char *regionName,
Array<bool> periodic);
/** read an unstructured mesh from a file */
void read_mesh(std::string meshFile, Array<bool> & periodicity);
/*@}*/
/** data that can be used for a subset of original mesh */
std::set<int> nullElements_;
/** faces upon which nodes are duplicated */
std::set<PAIR> cutFaces_;
std::set<int> cutEdge_;
/** workspace */
int nNodesPerElement_;
int nSD_;
int nElems_;
int nNodes_; /** number of global nodes */
int nNodesUnique_; /** number of unique nodes */
mutable int nIPsPerElement_;
mutable int nIPsPerFace_;
mutable FeIntQuadrature quadrature_;
mutable FIELDS::const_iterator _fieldItr_;
mutable FieldName _fieldName_;
/** sized arrays */
mutable DIAG_MAT _weights_;
mutable DENS_MAT _N_, _Nw_;
mutable DENS_MAT_VEC _dN_, _dNw_;
mutable DIAG_MAT _fweights_;
mutable DENS_MAT _fN_;
mutable DENS_MAT_VEC _fdN_, _nN_;
/** unsized arrays */
mutable DENS_MAT _Nmat_;
mutable FIELD_MATS _fieldsAtIPs_;
mutable GRAD_FIELD_MATS _gradFieldsAtIPs_;
mutable DENS_MAT _Nfluxes_;
mutable AliasArray<int> _conn_;
mutable DENS_MAT_VEC _Bfluxes_;
/** output object */
OutputManager outputManager_;
};
}; // end namespace ATC
#endif
|