#include namespace ngsolve { using namespace ngsolve; PDE :: PDE (const MeshAccess & ama) : ma (ama) { levelsolved = -1; } PDE :: ~PDE() { for (int i = 0; i < coefficients.Size(); i++) delete coefficients[i]; for (int i = 0; i < spaces.Size(); i++) delete spaces[i]; for (int i = 0; i < gridfunctions.Size(); i++) delete gridfunctions[i]; for (int i = 0; i < bilinearforms.Size(); i++) delete bilinearforms[i]; for (int i = 0; i < linearforms.Size(); i++) delete linearforms[i]; for (int i = 0; i < preconditioners.Size(); i++) delete preconditioners[i]; for (int i = 0; i < numprocs.Size(); i++) delete numprocs[i]; Ng_ClearSolutionData (); } void PDE :: SavePDE (const string & filename) { ; } void PDE :: SaveSolution (const string & filename) { // ofstream outfile(filename.c_str()); ofstream outfile(filename.c_str(), ios_base::binary); for (int i = 0; i < gridfunctions.Size(); i++) gridfunctions[i]->GetVector().Save (outfile); } /// void PDE :: LoadSolution (const string & filename) { #ifdef NETGEN_ELTRANS int geometryorder = 1; if (constants.Used ("geometryorder")) geometryorder = int (constants["geometryorder"]); Ng_HighOrder (geometryorder); #endif // ifstream infile(filename.c_str()); ifstream infile(filename.c_str(), ios_base::binary); for (int i = 0; i < spaces.Size(); i++) spaces[i]->Update(); for (int i = 0; i < gridfunctions.Size(); i++) { gridfunctions[i]->Update(); gridfunctions[i]->GetVector().Load (infile); } } void PDE :: PrintReport (ostream & ost) { ost << endl << "PDE Description:" << endl; for (int i = 0; i < constants.Size(); i++) ost << "constant " << constants.GetName(i) << " = " << constants[i] << endl; for (int i = 0; i < variables.Size(); i++) ost << "variable " << variables.GetName(i) << " = " << variables[i] << endl; ost << endl; ost << "Coefficients:" << endl << "-------------" << endl; for (int i = 0; i < coefficients.Size(); i++) { ost << "coefficient " << coefficients.GetName(i) << ":" << endl; coefficients[i]->PrintReport (ost); } ost << endl << "Spaces:" << endl << "-------" << endl; for (int i = 0; i < spaces.Size(); i++) { ost << "space " << spaces.GetName(i) << ":" << endl; spaces[i]->PrintReport (ost); } ost << endl << "Bilinear-forms:" << endl << "---------------" << endl; for (int i = 0; i < bilinearforms.Size(); i++) { ost << "bilinear-form " << bilinearforms.GetName(i) << ":" << endl; bilinearforms[i]->PrintReport (ost); } ost << endl << "Linear-forms:" << endl << "-------------" << endl; for (int i = 0; i < linearforms.Size(); i++) { ost << "linear-form " << linearforms.GetName(i) << ":" << endl; linearforms[i]->PrintReport (ost); } ost << endl << "Grid-functions:" << endl << "---------------" << endl; for (int i = 0; i < gridfunctions.Size(); i++) { ost << "grid-function " << gridfunctions.GetName(i) << ":" << endl; gridfunctions[i]->PrintReport (ost); } ost << endl << "Preconditioners:" << endl << "----------------" << endl; for (int i = 0; i < preconditioners.Size(); i++) { ost << "preconditioner " << preconditioners.GetName(i) << ":" << endl; preconditioners[i]->PrintReport (ost); } ost << endl << "Numprocs:" << endl << "---------" << endl; for (int i = 0; i < numprocs.Size(); i++) { ost << "numproc " << numprocs.GetName(i) << ":" << endl; numprocs[i]->PrintReport (ost); } } void PDE :: PrintMemoryUsage (ostream & ost) { ARRAY memuse; for (int i = 0; i < spaces.Size(); i++) spaces[i]->MemoryUsage (memuse); for (int i = 0; i < bilinearforms.Size(); i++) bilinearforms[i]->MemoryUsage (memuse); for (int i = 0; i < linearforms.Size(); i++) linearforms[i]->MemoryUsage (memuse); for (int i = 0; i < gridfunctions.Size(); i++) gridfunctions[i]->MemoryUsage (memuse); for (int i = 0; i < preconditioners.Size(); i++) preconditioners[i]->MemoryUsage (memuse); int sumbytes = 0, sumblocks = 0; for (int i = 0; i < memuse.Size(); i++) { ost << memuse[i]->Name() << ": " << memuse[i]->NBytes() << " bytes in " << memuse[i]->NBlocks() << " blocks." << endl; sumbytes += memuse[i]->NBytes(); sumblocks += memuse[i]->NBlocks(); } cout << "total bytes " << sumbytes << " in " << sumblocks << " blocks." << endl; } bool PDE :: ConstantUsed (const string & aname) const { return constants.Used(aname); } double PDE :: GetConstant (const string & name, bool opt) const { if (constants.Used(name)) return constants[name]; if (opt) return 0; stringstream str; str << "Constant '" << name << "' not defined\n"; throw Exception (str.str()); } bool PDE :: VariableUsed (const string & aname) const { return variables.Used(aname); } double & PDE :: GetVariable (const string & name, bool opt) { if (variables.Used(name)) return variables[name]; static double dummy; if (opt) return dummy; stringstream str; str << "Variable '" << name << "' not defined\n"; throw Exception (str.str()); } CoefficientFunction * PDE :: GetCoefficientFunction (const string & name, bool opt) { if (coefficients.Used(name)) return coefficients[name]; if (opt) return 0; stringstream str; str << "CoefficientFunction '" << name << "' not defined\n"; throw Exception (str.str()); } FESpace * PDE :: GetFESpace (const string & name, bool opt) { if (spaces.Used(name)) return spaces[name]; if (opt) return 0; stringstream str; str << "FESpace '" << name << "' not defined\n"; throw Exception (str.str()); } GridFunction * PDE :: GetGridFunction (const string & name, bool opt) { if (gridfunctions.Used(name)) return gridfunctions[name]; if (opt) return 0; stringstream str; str << "GridFunction '" << name << "' not defined\n"; throw Exception (str.str()); } BilinearForm * PDE :: GetBilinearForm (const string & name, bool opt) { if (bilinearforms.Used(name)) return bilinearforms[name]; if (opt) return 0; stringstream str; str << "Bilinear-form '" << name << "' not defined\n"; throw Exception (str.str()); } /* BEMBilinearForm * PDE :: GetBEMBilinearForm (const string & name, bool opt) { if (bembilinearforms.Used(name)) return bembilinearforms[name]; if (opt) return 0; stringstream str; str << "BEM Bilinear-form '" << name << "' not defined\n"; throw Exception (str.str()); } */ LinearForm * PDE :: GetLinearForm (const string & name, bool opt) { if (linearforms.Used(name)) return linearforms[name]; if (opt) return 0; stringstream str; str << "Linear-form '" << name << "' not defined\n"; throw Exception (str.str()); } Preconditioner * PDE :: GetPreconditioner (const string & name, bool opt) { if (preconditioners.Used(name)) return preconditioners[name]; if (opt) return 0; stringstream str; str << "Preconditioner '" << name << "' not defined\n"; throw Exception (str.str()); } NumProc * PDE :: GetNumProc (const string & name, bool opt) { if (numprocs.Used(name)) return numprocs[name]; if (opt) return 0; stringstream str; str << "Numproc '" << name << "' not defined\n"; throw Exception (str.str()); } const CoefficientFunction * PDE :: GetCoefficientFunction (const string & name, bool opt) const { if (coefficients.Used(name)) return coefficients[name]; if (opt) return 0; stringstream str; str << "CoefficientFunction '" << name << "' not defined\n"; throw Exception (str.str()); } const FESpace * PDE :: GetFESpace (const string & name, bool opt) const { if (spaces.Used(name)) return spaces[name]; if (opt) return 0; stringstream str; str << "FESpace '" << name << "' not defined\n"; throw Exception (str.str()); } const GridFunction * PDE :: GetGridFunction (const string & name, bool opt) const { if (gridfunctions.Used(name)) return gridfunctions[name]; if (opt) return 0; stringstream str; str << "Grid-function '" << name << "' not defined\n"; throw Exception (str.str()); } const BilinearForm * PDE :: GetBilinearForm (const string & name, bool opt) const { if (bilinearforms.Used(name)) return bilinearforms[name]; if (opt) return 0; stringstream str; str << "Bilinear-form '" << name << "' not defined\n"; throw Exception (str.str()); } /* const BEMBilinearForm * PDE :: GetBEMBilinearForm (const string & name, bool opt) const { if (bembilinearforms.Used(name)) return bembilinearforms[name]; if (opt) return 0; stringstream str; str << "BEM Bilinear-form '" << name << "' not defined\n"; throw Exception (str.str()); } */ const LinearForm * PDE :: GetLinearForm (const string & name, bool opt) const { if (linearforms.Used(name)) return linearforms[name]; if (opt) return 0; stringstream str; str << "Linear-form '" << name << "' not defined\n"; throw Exception (str.str()); } const Preconditioner * PDE :: GetPreconditioner (const string & name, bool opt) const { if (preconditioners.Used(name)) return preconditioners[name]; if (opt) return 0; stringstream str; str << "Preconditioner '" << name << "' not defined\n"; throw Exception (str.str()); } const NumProc * PDE :: GetNumProc (const string & name, bool opt) const { if (numprocs.Used(name)) return numprocs[name]; if (opt) return 0; stringstream str; str << "Numproc '" << name << "' not defined\n"; throw Exception (str.str()); } void PDE :: SolveBVP () { long int heapsize = 1000000; if (constants.Used ("heapsize")) heapsize = (long int) (constants["heapsize"]); LocalHeap lh(heapsize); clock_t starttime, endtime; starttime = clock(); if (levelsolved >= 0) { if (constants.Used ("refinehp")) Ng_Refine(NG_REFINE_HP); else if (constants.Used ("refinep")) Ng_Refine(NG_REFINE_P); else Ng_Refine(NG_REFINE_H); } int secondorder = 0; if (constants.Used ("secondorder")) secondorder = int (constants["secondorder"]); if (constants.Used ("hpref") && levelsolved == -1) Ng_HPRefinement (int (constants["hpref"])); int geometryorder = 0; #ifdef NETGEN_ELTRANS geometryorder = 1; secondorder = 0; if (constants.Used ("geometryorder")) geometryorder = int (constants["geometryorder"]); #endif if (secondorder) Ng_SecondOrder(); if (geometryorder) Ng_HighOrder (geometryorder); cout << "Solve at level " << ma.GetNLevels()-1 << ", NE = " << ma.GetNE() << ", NP = " << ma.GetNP() << endl; for (int i = 0; i < spaces.Size(); i++) { try { cout << "Update " << spaces[i] -> GetClassName() << " " << spaces[i] -> GetName () << flush; spaces[i] -> Update(lh); lh.CleanUp(); if (spaces[i]->GetDimension() == 1) cout << ", ndof = " << spaces[i] -> GetNDof() << endl; else cout << ", ndof = " << spaces[i] -> GetDimension() << " x " << spaces[i] -> GetNDof() << endl; } catch (exception & e) { throw Exception (e.what() + string ("\nthrown by update space ") + string (spaces.GetName(i))); } catch (Exception & e) { e.Append (string ("\nthrown by update space ") + string (spaces.GetName(i))); throw; } } for (int i = 0; i < gridfunctions.Size(); i++) { try { cout << "Update gridfunction " << gridfunctions[i]->GetName() << endl; gridfunctions[i]->Update(); } catch (exception & e) { throw Exception (e.what() + string ("\nthrown by update grid-function ") + string (gridfunctions.GetName(i))); } catch (Exception & e) { e.Append (string ("\nthrown by update grid-function ") + string (gridfunctions.GetName(i))); throw; } } for (int i = 0; i < linearforms.Size(); i++) { try { cout << "Update linear-form " << linearforms[i]->GetName() << endl; linearforms[i]->Assemble(lh); lh.CleanUp(); } catch (exception & e) { throw Exception (e.what() + string ("\nthrown by update linear-form ") + string (linearforms.GetName(i))); } catch (Exception & e) { e.Append (string ("\nthrown by update linear-form ") + string (linearforms.GetName(i))); throw; } } for (int i = 0; i < bilinearforms.Size(); i++) { try { cout << "Update bilinear-form " << bilinearforms[i]->GetName() << endl; bilinearforms[i]->Assemble(lh); lh.CleanUp(); } catch (exception & e) { throw Exception (e.what() + string ("\nthrown by update bilinear-form ") + string (bilinearforms.GetName(i))); } catch (Exception & e) { e.Append (string ("\nthrown by update bilinear-form ") + string (bilinearforms.GetName(i))); throw; } } /* for (int i = 0; i < bembilinearforms.Size(); i++) { try { cout << "Update BEM bilinear-form " << bembilinearforms.GetName(i) << endl; bembilinearforms[i]->Assemble(); } catch (exception & e) { throw Exception (e.what() + string ("\nthrown by update bem bilinear-form ") + string (bembilinearforms.GetName(i))); } catch (Exception & e) { e.Append (string ("thrown by update bem bilinear-form ") + string (bembilinearforms.GetName(i)) + string ("\n")); throw; } } */ for (int i = 0; i < preconditioners.Size(); i++) { try { cout << "Update " << preconditioners[i]->ClassName() << " " << preconditioners.GetName(i) << endl; preconditioners[i]->Update(); // preconditioners[i]->Test(); } catch (exception & e) { throw Exception (e.what() + string ("\nthrown by update preconditioner ") + string (preconditioners.GetName(i))); } catch (Exception & e) { e.Append (string ("\nthrown by update preconditioner ") + string (preconditioners.GetName(i))); throw; } } for (int i = 0; i < numprocs.Size(); i++) { try { cout << "Call numproc " << numprocs[i]->GetClassName() << " " << numprocs[i]->GetName() << endl; numprocs[i]->Do(lh); lh.CleanUp(); } catch (exception & e) { throw Exception (e.what() + string ("\nthrown by update numproc ") + string (numprocs.GetName(i))); } catch (Exception & e) { e.Append (string ("\nthrown by update numproc ") + string (numprocs.GetName(i))); throw; } } // set solution data for (int i = 0; i < gridfunctions.Size(); i++) gridfunctions[i]->Visualize(gridfunctions.GetName(i)); Ng_Redraw(); levelsolved++; cout << "Problem Solved" << endl; endtime = clock(); cout << "Total Time = " << double(endtime - starttime)/CLOCKS_PER_SEC << endl << endl; } /* CoefficientFunction * PDE :: GetCoefficientFunction (const string & name) { return coefficients.Get (name); } FESpace * PDE :: GetFESpace (const string & name) { return spaces.Get (name); } GridFunction * PDE :: GetGridFunction (const string & name) { return gridfunctions.Get (name); } BilinearForm * PDE :: GetBilinearForm (const string & name) { return bilinearforms.Get (name); } LinearForm * PDE :: GetLinearForm (const string & name) { return linearforms.Get (name); } Preconditioner * PDE :: GetPreconditioner (const string & name) { if (preconditioners.Used(name)) return preconditioners.Get (name); else return NULL; } NumProc * PDE :: GetNumProc (const string & name) { return numprocs.Get (name); } const CoefficientFunction * PDE :: GetCoefficientFunction (const string & name) const { return coefficients.Get (name); } const FESpace * PDE :: GetFESpace (const string & name) const { return spaces.Get (name); } const GridFunction * PDE :: GetGridFunction (const string & name) const { return gridfunctions.Get (name); } const BilinearForm * PDE :: GetBilinearForm (const string & name) const { return bilinearforms.Get (name); } const LinearForm * PDE :: GetLinearForm (const string & name) const { return linearforms.Get (name); } const Preconditioner * PDE :: GetPreconditioner (const string & name) const { if (preconditioners.Used(name)) return preconditioners.Get (name); else return NULL; } const NumProc * PDE :: GetNumProc (const string & name) const { return numprocs.Get (name); } */ void PDE :: AddConstant (const string & name, double val) { cout << "add constant " << name << " = " << val << endl; constants.Set (name.c_str(), val); } void PDE :: AddVariable (const string & name, double val) { cout << "add variable " << name << " = " << val << endl; variables.Set (name.c_str(), val); } void PDE :: AddCoefficientFunction (const string & name, CoefficientFunction* fun) { cout << "add coefficient-function, name = " << name << endl; coefficients.Set (name.c_str(), fun); } FESpace * PDE :: AddFESpace (const string & name, Flags & flags) { cout << "add fespace " << name << flush; FESpace * space = 0; if (flags.GetDefineFlag ("vec")) flags.SetFlag ("dim", ma.GetDimension()); if (flags.GetDefineFlag ("tensor")) flags.SetFlag ("dim", sqr (ma.GetDimension())); if (flags.GetDefineFlag ("symtensor")) flags.SetFlag ("dim", ma.GetDimension()*(ma.GetDimension()+1) / 2); int order = int(flags.GetNumFlag ("order", -1)); for (int i = 0; i < GetFESpaceClasses().GetFESpaces().Size(); i++) { if (flags.GetDefineFlag (GetFESpaceClasses().GetFESpaces()[i]->name)) space = GetFESpaceClasses().GetFESpaces()[i]->creator (ma, flags); } if (space) ; else if (flags.GetDefineFlag ("compound")) { const ARRAY & spacenames = flags.GetStringListFlag ("spaces"); ARRAY spaces (spacenames.Size()); for (int i = 0; i < spaces.Size(); i++) spaces[i] = GetFESpace (spacenames[i]); space = new CompoundFESpace (GetMeshAccess(), spaces); } else { if (order <= 1) space = new NodalFESpace (GetMeshAccess(), flags); else space = new H1HighOrderFESpace (GetMeshAccess(), flags); } if (flags.GetDefineFlag ("bem")) space->SetBEM (1); if (flags.NumListFlagDefined ("dirichletboundaries")) { BitArray dirbnds(ma.GetNBoundaries()); dirbnds.Clear(); const ARRAY & array = flags.GetNumListFlag ("dirichletboundaries"); for (int i = 0; i < array.Size(); i++) dirbnds.Set (int(array[i])-1); space->SetDirichletBoundaries (dirbnds); } if (flags.NumListFlagDefined ("domains")) { BitArray definedon(ma.GetNDomains()); definedon.Clear(); const ARRAY & domains = flags.GetNumListFlag ("domains"); for (int i = 0; i < domains.Size(); i++) definedon.Set (int(domains[i])-1); space->SetDefinedOn (definedon); } if (flags.NumListFlagDefined ("boundaries")) { BitArray definedon(ma.GetNBoundaries()); definedon.Clear(); const ARRAY & boundaries = flags.GetNumListFlag ("boundaries"); for (int i = 0; i < boundaries.Size(); i++) definedon.Set (int(boundaries[i])-1); space->SetDefinedOnBoundary (definedon); } space->SetName (name); spaces.Set (name, space); cout << ", type = " << space->GetClassName() << endl; return space; } GridFunction * PDE :: AddGridFunction (const string & name, Flags & flags) { cout << "add grid-function " << name << endl; string spacename = flags.GetStringFlag ("fespace", ""); if (!spaces.Used (spacename)) { throw Exception (string ("Gridfuncton '") + name + "' uses undefined space '" + spacename + "'"); } const FESpace * space = GetFESpace(spacename); GridFunction * gf = CreateGridFunction (space, name, flags); gridfunctions.Set (name, gf); return gf; } BilinearForm * PDE :: AddBilinearForm (const string & name, Flags & flags) { cout << "add bilinear-form " << name << endl; string spacename = flags.GetStringFlag ("fespace", ""); if (!spaces.Used (spacename)) { cerr << "space " << spacename << " not defined " << endl; return 0; } const FESpace * space = spaces[spacename]; const FESpace * space2 = NULL; if (flags.StringFlagDefined ("fespace2")) space2 = spaces[flags.GetStringFlag ("fespace2", "")]; if (!space2) { bilinearforms.Set (name, CreateBilinearForm (space, name, flags)); } // else // bilinearforms.Set (name, new T_BilinearForm (*space, *space2, name)); if (flags.StringFlagDefined ("linearform")) bilinearforms[name] -> SetLinearForm (GetLinearForm (flags.GetStringFlag ("linearform", 0))); return bilinearforms[name]; } void PDE :: AddBEMElement (const string & name, Flags & flags) { cout << "add BEM-Element " << name << ", flags = " << flags << endl; // BilinearForm * blf = GetBilinearForm (flags.GetStringFlag ("bilinearform", NULL)); FESpace * spaced = spaces[flags.GetStringFlag ("fespaced", NULL)]; FESpace * spacen = spaces[flags.GetStringFlag ("fespacen", NULL)]; if (flags.GetDefineFlag ("helmholtz")) { spaced -> specialelements.Append (new Helmholtz_BEMElement (*spaced, *spacen, flags.GetNumFlag ("k", 1.0))); } else spaced -> specialelements.Append (new S_LaplaceBEMElement (*spaced, *spacen)); } LinearForm * PDE :: AddLinearForm (const string & name, Flags & flags) { cout << "add linear-form " << name << endl; // << ", flags = " << endl; // flags.PrintFlags (cout); const char * spacename = flags.GetStringFlag ("fespace", ""); if (!spaces.Used (spacename)) { cerr << "space " << spacename << " not defined " << endl; return 0; } const FESpace * space = spaces[spacename]; /* LinearForm * lf; CreateVecObject2 (lf, T_LinearForm, space->GetDimension(), space->IsComplex(), *space, name); linearforms.Set (name, lf); */ linearforms.Set (name, CreateLinearForm (space, name, flags)); return linearforms[name]; } Preconditioner * PDE :: AddPreconditioner (const string & name, Flags & flags) { cout << "add preconditioner " << name << flush; // flags.PrintFlags (cout); Preconditioner * pre = NULL; const char * type = flags.GetStringFlag ("type", NULL); if (strcmp (type, "multigrid") == 0) pre = new MGPreconditioner (this, flags); else if (strcmp (type, "direct") == 0) pre = new DirectPreconditioner (this, flags); else if (strcmp (type, "local") == 0) pre = new LocalPreconditioner (this, flags); // else if (strcmp (type, "vefc") == 0) // pre = new VEFC_Preconditioner (this, flags); else if (strcmp (type, "twolevel") == 0) pre = new TwoLevelPreconditioner (this, flags); else if (strcmp (type, "complex") == 0) pre = new ComplexPreconditioner (this, flags); else if (strcmp (type, "chebychev") == 0) pre = new ChebychevPreconditioner (this, flags); // else if (strcmp (type, "constrained") == 0) // pre = new ConstrainedPreconditioner (this, flags); else if (strcmp (type, "amg") == 0) pre = new CommutingAMGPreconditioner (this, flags); else if (strcmp (type, "dndd") == 0) pre = new DNDDPreconditioner (this, flags); // changed 08/19/2003, Bachinger else if (strcmp (type, "nonsymmetric") == 0) pre = new NonsymmetricPreconditioner (this, flags); else for (int i = 0; i < GetPreconditionerClasses().GetPreconditioners().Size(); i++) { if (flags.GetDefineFlag (GetPreconditionerClasses().GetPreconditioners()[i]->name)) pre = GetPreconditionerClasses().GetPreconditioners()[i]->creator (*this, flags); } if (!pre) throw Exception ("Unknown preconditioner type " + string(type)); preconditioners.Set (name, pre); cout << ", type = " << pre->ClassName() << endl; return pre; } void PDE :: AddNumProc (const string & name, NumProc * np) { cout << "add numproc " << name << ", type = " << np->GetClassName() << endl; np->SetName (name); numprocs.Set (name, np); } void PDE :: AddBilinearFormIntegrator (const string & name, BilinearFormIntegrator * part) { BilinearForm * form = GetBilinearForm (name); if (form && part) { form->AddIntegrator (part); cout << "integrator " << part->Name() << endl; } else { cerr << "Bilinearform = " << form << ", part = " << part << endl; } } void PDE :: AddLinearFormIntegrator (const string & name, LinearFormIntegrator * part) { LinearForm * form = GetLinearForm (name); if (form && part) { form->AddIntegrator (part); cout << "integrator " << part->Name() << endl; } else { cerr << "Linearform = " << form << ", part = " << part << endl; } } }