From 016ba45b02b4ca9beb57022aa69a71f9ed0b07ec Mon Sep 17 00:00:00 2001 From: Axel Larsson Date: Tue, 13 Feb 2024 12:05:56 -0800 Subject: [PATCH] Boilerplate for working BC --- src/interface_form.cpp | 759 +++++++++++++++++++------------------- src/interfaceinteg.cpp | 70 ++-- src/linelast_solver.cpp | 181 +++++++-- src/main_workflow.cpp | 251 +++++++++---- src/multiblock_solver.cpp | 4 +- test/dg_integ_mms.cpp | 2 +- 6 files changed, 737 insertions(+), 530 deletions(-) diff --git a/src/interface_form.cpp b/src/interface_form.cpp index d32395da..94b51fca 100644 --- a/src/interface_form.cpp +++ b/src/interface_form.cpp @@ -10,473 +10,482 @@ using namespace std; namespace mfem { -InterfaceForm::InterfaceForm( - Array &meshes_, Array &fes_, TopologyHandler *topol_) - : meshes(meshes_), fes(fes_), topol_handler(topol_), numSub(meshes_.Size()) -{ - assert(fes_.Size() == numSub); + InterfaceForm::InterfaceForm( + Array &meshes_, Array &fes_, TopologyHandler *topol_) + : meshes(meshes_), fes(fes_), topol_handler(topol_), numSub(meshes_.Size()) + { + assert(fes_.Size() == numSub); - block_offsets.SetSize(numSub + 1); - block_offsets = 0; - for (int i = 1; i < numSub + 1; i++) - block_offsets[i] = fes[i-1]->GetTrueVSize(); - block_offsets.PartialSum(); -} + block_offsets.SetSize(numSub + 1); + block_offsets = 0; + for (int i = 1; i < numSub + 1; i++) + block_offsets[i] = fes[i - 1]->GetTrueVSize(); + block_offsets.PartialSum(); + } -InterfaceForm::~InterfaceForm() -{ - DeletePointers(fnfi); -} + InterfaceForm::~InterfaceForm() + { + DeletePointers(fnfi); + } -void InterfaceForm::AssembleInterfaceMatrixes(Array2D &mats) const -{ - assert(mats.NumRows() == numSub); - assert(mats.NumCols() == numSub); - for (int i = 0; i < numSub; i++) - for (int j = 0; j < numSub; j++) assert(mats(i, j)); - - const PortInfo *pInfo; - Array midx(2); - Array2D mats_p(2,2); - Mesh *mesh1, *mesh2; - FiniteElementSpace *fes1, *fes2; - Array* interface_infos; - - for (int p = 0; p < topol_handler->GetNumPorts(); p++) + void InterfaceForm::AssembleInterfaceMatrixes(Array2D &mats) const { - pInfo = topol_handler->GetPortInfo(p); + assert(mats.NumRows() == numSub); + assert(mats.NumCols() == numSub); + for (int i = 0; i < numSub; i++) + for (int j = 0; j < numSub; j++) + assert(mats(i, j)); + + const PortInfo *pInfo; + Array midx(2); + Array2D mats_p(2, 2); + Mesh *mesh1, *mesh2; + FiniteElementSpace *fes1, *fes2; + Array *interface_infos; + + for (int p = 0; p < topol_handler->GetNumPorts(); p++) + { + pInfo = topol_handler->GetPortInfo(p); - midx[0] = pInfo->Mesh1; - midx[1] = pInfo->Mesh2; - - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) mats_p(i, j) = mats(midx[i], midx[j]); + midx[0] = pInfo->Mesh1; + midx[1] = pInfo->Mesh2; - mesh1 = meshes[midx[0]]; - mesh2 = meshes[midx[1]]; + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + mats_p(i, j) = mats(midx[i], midx[j]); - fes1 = fes[midx[0]]; - fes2 = fes[midx[1]]; + mesh1 = meshes[midx[0]]; + mesh2 = meshes[midx[1]]; - interface_infos = topol_handler->GetInterfaceInfos(p); - AssembleInterfaceMatrix(mesh1, mesh2, fes1, fes2, interface_infos, mats_p); - } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) -} + fes1 = fes[midx[0]]; + fes2 = fes[midx[1]]; -void InterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const -{ - x_tmp.Update(const_cast(x), block_offsets); - y_tmp.Update(y, block_offsets); - - Array midx(2); - Mesh *mesh1, *mesh2; - FiniteElementSpace *fes1, *fes2; + interface_infos = topol_handler->GetInterfaceInfos(p); + AssembleInterfaceMatrix(mesh1, mesh2, fes1, fes2, interface_infos, mats_p); + } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) + } - for (int p = 0; p < topol_handler->GetNumPorts(); p++) + void InterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const { - const PortInfo *pInfo = topol_handler->GetPortInfo(p); + x_tmp.Update(const_cast(x), block_offsets); + y_tmp.Update(y, block_offsets); - midx[0] = pInfo->Mesh1; - midx[1] = pInfo->Mesh2; + Array midx(2); + Mesh *mesh1, *mesh2; + FiniteElementSpace *fes1, *fes2; - mesh1 = meshes[midx[0]]; - mesh2 = meshes[midx[1]]; + for (int p = 0; p < topol_handler->GetNumPorts(); p++) + { + const PortInfo *pInfo = topol_handler->GetPortInfo(p); - fes1 = fes[midx[0]]; - fes2 = fes[midx[1]]; + midx[0] = pInfo->Mesh1; + midx[1] = pInfo->Mesh2; - Array* const interface_infos = topol_handler->GetInterfaceInfos(p); - AssembleInterfaceVector(mesh1, mesh2, fes1, fes2, interface_infos, - x_tmp.GetBlock(midx[0]), x_tmp.GetBlock(midx[1]), - y_tmp.GetBlock(midx[0]), y_tmp.GetBlock(midx[1])); - } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) + mesh1 = meshes[midx[0]]; + mesh2 = meshes[midx[1]]; - for (int i=0; i < y_tmp.NumBlocks(); ++i) - y_tmp.GetBlock(i).SyncAliasMemory(y); -} + fes1 = fes[midx[0]]; + fes2 = fes[midx[1]]; -void InterfaceForm::InterfaceGetGradient(const Vector &x, Array2D &mats) const -{ - assert(mats.NumRows() == numSub); - assert(mats.NumCols() == numSub); - for (int i = 0; i < numSub; i++) - for (int j = 0; j < numSub; j++) - assert(mats(i, j)); + Array *const interface_infos = topol_handler->GetInterfaceInfos(p); + AssembleInterfaceVector(mesh1, mesh2, fes1, fes2, interface_infos, + x_tmp.GetBlock(midx[0]), x_tmp.GetBlock(midx[1]), + y_tmp.GetBlock(midx[0]), y_tmp.GetBlock(midx[1])); + } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) - x_tmp.Update(const_cast(x), block_offsets); - - Array midx(2); - Array2D mats_p(2,2); - Mesh *mesh1, *mesh2; - FiniteElementSpace *fes1, *fes2; + for (int i = 0; i < y_tmp.NumBlocks(); ++i) + y_tmp.GetBlock(i).SyncAliasMemory(y); + } - for (int p = 0; p < topol_handler->GetNumPorts(); p++) + void InterfaceForm::InterfaceGetGradient(const Vector &x, Array2D &mats) const { - const PortInfo *pInfo = topol_handler->GetPortInfo(p); - - midx[0] = pInfo->Mesh1; - midx[1] = pInfo->Mesh2; + assert(mats.NumRows() == numSub); + assert(mats.NumCols() == numSub); + for (int i = 0; i < numSub; i++) + for (int j = 0; j < numSub; j++) + assert(mats(i, j)); - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) - mats_p(i, j) = mats(midx[i], midx[j]); + x_tmp.Update(const_cast(x), block_offsets); + + Array midx(2); + Array2D mats_p(2, 2); + Mesh *mesh1, *mesh2; + FiniteElementSpace *fes1, *fes2; - mesh1 = meshes[midx[0]]; - mesh2 = meshes[midx[1]]; + for (int p = 0; p < topol_handler->GetNumPorts(); p++) + { + const PortInfo *pInfo = topol_handler->GetPortInfo(p); - fes1 = fes[midx[0]]; - fes2 = fes[midx[1]]; + midx[0] = pInfo->Mesh1; + midx[1] = pInfo->Mesh2; - Array* const interface_infos = topol_handler->GetInterfaceInfos(p); - AssembleInterfaceGrad(mesh1, mesh2, fes1, fes2, interface_infos, - x_tmp.GetBlock(midx[0]), x_tmp.GetBlock(midx[1]), mats_p); - } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) -} + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + mats_p(i, j) = mats(midx[i], midx[j]); -void InterfaceForm::AssembleInterfaceMatrix( - Mesh *mesh1, Mesh *mesh2, FiniteElementSpace *fes1, FiniteElementSpace *fes2, - Array *interface_infos, Array2D &mats) const -{ - Array2D elemmats(2, 2); - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) elemmats(i, j) = new DenseMatrix; - FaceElementTransformations *tr1, *tr2; - const FiniteElement *fe1, *fe2; - Array *> vdofs(2); - vdofs[0] = new Array; - vdofs[1] = new Array; - - for (int bn = 0; bn < interface_infos->Size(); bn++) - { - InterfaceInfo *if_info = &((*interface_infos)[bn]); + mesh1 = meshes[midx[0]]; + mesh2 = meshes[midx[1]]; + + fes1 = fes[midx[0]]; + fes2 = fes[midx[1]]; - topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); + Array *const interface_infos = topol_handler->GetInterfaceInfos(p); + AssembleInterfaceGrad(mesh1, mesh2, fes1, fes2, interface_infos, + x_tmp.GetBlock(midx[0]), x_tmp.GetBlock(midx[1]), mats_p); + } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) + } - if ((tr1 != NULL) && (tr2 != NULL)) + void InterfaceForm::AssembleInterfaceMatrix( + Mesh *mesh1, Mesh *mesh2, FiniteElementSpace *fes1, FiniteElementSpace *fes2, + Array *interface_infos, Array2D &mats) const + { + Array2D elemmats(2, 2); + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + elemmats(i, j) = new DenseMatrix; + FaceElementTransformations *tr1, *tr2; + const FiniteElement *fe1, *fe2; + Array *> vdofs(2); + vdofs[0] = new Array; + vdofs[1] = new Array; + for (int bn = 0; bn < interface_infos->Size(); bn++) { - fes1->GetElementVDofs(tr1->Elem1No, *vdofs[0]); - fes2->GetElementVDofs(tr2->Elem1No, *vdofs[1]); - // Both domains will have the adjacent element as Elem1. - fe1 = fes1->GetFE(tr1->Elem1No); - fe2 = fes2->GetFE(tr2->Elem1No); + InterfaceInfo *if_info = &((*interface_infos)[bn]); - for (int itg = 0; itg < fnfi.Size(); itg++) - { - assert(fnfi[itg]); - fnfi[itg]->AssembleInterfaceMatrix(*fe1, *fe2, *tr1, *tr2, elemmats); + topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); - for (int i = 0; i < 2; i++) { - for (int j = 0; j < 2; j++) { - mats(i, j)->AddSubMatrix(*vdofs[i], *vdofs[j], *elemmats(i,j), skip_zeros); + if ((tr1 != NULL) && (tr2 != NULL)) + { + fes1->GetElementVDofs(tr1->Elem1No, *vdofs[0]); + fes2->GetElementVDofs(tr2->Elem1No, *vdofs[1]); + // Both domains will have the adjacent element as Elem1. + fe1 = fes1->GetFE(tr1->Elem1No); + fe2 = fes2->GetFE(tr2->Elem1No); + + for (int itg = 0; itg < fnfi.Size(); itg++) + { + assert(fnfi[itg]); + fnfi[itg]->AssembleInterfaceMatrix(*fe1, *fe2, *tr1, *tr2, elemmats); + for (int i = 0; i < 2; i++) + { + for (int j = 0; j < 2; j++) + { + mats(i, j)->AddSubMatrix(*vdofs[i], *vdofs[j], *elemmats(i, j), skip_zeros); + } } } - } - } // if ((tr1 != NULL) && (tr2 != NULL)) - } // for (int bn = 0; bn < interface_infos.Size(); bn++) + } // if ((tr1 != NULL) && (tr2 != NULL)) + } // for (int bn = 0; bn < interface_infos.Size(); bn++) - DeletePointers(elemmats); - DeletePointers(vdofs); -} - -void InterfaceForm::AssembleInterfaceMatrixAtPort( - const int p, Array &fes_comp, Array2D &mats_p) const -{ - const int num_ref_ports = topol_handler->GetNumRefPorts(); - - assert(topol_handler->GetType() == TopologyHandlerMode::COMPONENT); - assert((p >= 0) && (p < num_ref_ports)); - DeletePointers(mats_p); - mats_p.SetSize(2, 2); + DeletePointers(elemmats); + DeletePointers(vdofs); + } - int c1, c2; - topol_handler->GetComponentPair(p, c1, c2); - Mesh *comp1 = topol_handler->GetComponentMesh(c1); - Mesh *comp2 = topol_handler->GetComponentMesh(c2); + void InterfaceForm::AssembleInterfaceMatrixAtPort( + const int p, Array &fes_comp, Array2D &mats_p) const + { + const int num_ref_ports = topol_handler->GetNumRefPorts(); - // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. - // Need to use two copied instances. - Mesh mesh1(*comp1); - Mesh mesh2(*comp2); + assert(topol_handler->GetType() == TopologyHandlerMode::COMPONENT); + assert((p >= 0) && (p < num_ref_ports)); + DeletePointers(mats_p); + mats_p.SetSize(2, 2); - Array c_idx(2); - c_idx[0] = c1; - c_idx[1] = c2; + int c1, c2; + topol_handler->GetComponentPair(p, c1, c2); + Mesh *comp1 = topol_handler->GetComponentMesh(c1); + Mesh *comp2 = topol_handler->GetComponentMesh(c2); - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) - mats_p(i, j) = new SparseMatrix(fes_comp[c_idx[i]]->GetTrueVSize(), fes_comp[c_idx[j]]->GetTrueVSize()); + // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. + // Need to use two copied instances. + Mesh mesh1(*comp1); + Mesh mesh2(*comp2); - Array *if_infos = topol_handler->GetRefInterfaceInfos(p); + Array c_idx(2); + c_idx[0] = c1; + c_idx[1] = c2; - // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. - // Need to use two copied instances. - AssembleInterfaceMatrix(&mesh1, &mesh2, fes_comp[c1], fes_comp[c2], if_infos, mats_p); + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + mats_p(i, j) = new SparseMatrix(fes_comp[c_idx[i]]->GetTrueVSize(), fes_comp[c_idx[j]]->GetTrueVSize()); - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) mats_p(i, j)->Finalize(); -} + Array *if_infos = topol_handler->GetRefInterfaceInfos(p); -void InterfaceForm::AssembleInterfaceVector(Mesh *mesh1, Mesh *mesh2, - FiniteElementSpace *fes1, FiniteElementSpace *fes2, Array *interface_infos, - const Vector &x1, const Vector &x2, Vector &y1, Vector &y2) const -{ - assert(x1.Size() == fes1->GetTrueVSize()); - assert(x2.Size() == fes2->GetTrueVSize()); - assert(y1.Size() == x1.Size()); - assert(y2.Size() == x2.Size()); + // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. + // Need to use two copied instances. + AssembleInterfaceMatrix(&mesh1, &mesh2, fes_comp[c1], fes_comp[c2], if_infos, mats_p); - FaceElementTransformations *tr1, *tr2; - const FiniteElement *fe1, *fe2; - Array vdofs1, vdofs2; - Vector el_x1, el_x2, el_y1, el_y2; + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + mats_p(i, j)->Finalize(); + } - for (int bn = 0; bn < interface_infos->Size(); bn++) + void InterfaceForm::AssembleInterfaceVector(Mesh *mesh1, Mesh *mesh2, + FiniteElementSpace *fes1, FiniteElementSpace *fes2, Array *interface_infos, + const Vector &x1, const Vector &x2, Vector &y1, Vector &y2) const { - InterfaceInfo *if_info = &((*interface_infos)[bn]); + assert(x1.Size() == fes1->GetTrueVSize()); + assert(x2.Size() == fes2->GetTrueVSize()); + assert(y1.Size() == x1.Size()); + assert(y2.Size() == x2.Size()); - topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); + FaceElementTransformations *tr1, *tr2; + const FiniteElement *fe1, *fe2; + Array vdofs1, vdofs2; + Vector el_x1, el_x2, el_y1, el_y2; - if ((tr1 != NULL) && (tr2 != NULL)) + for (int bn = 0; bn < interface_infos->Size(); bn++) { - fes1->GetElementVDofs(tr1->Elem1No, vdofs1); - fes2->GetElementVDofs(tr2->Elem1No, vdofs2); - - x1.GetSubVector(vdofs1, el_x1); - x2.GetSubVector(vdofs2, el_x2); + InterfaceInfo *if_info = &((*interface_infos)[bn]); - // Both domains will have the adjacent element as Elem1. - fe1 = fes1->GetFE(tr1->Elem1No); - fe2 = fes2->GetFE(tr2->Elem1No); + topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); - for (int itg = 0; itg < fnfi.Size(); itg++) + if ((tr1 != NULL) && (tr2 != NULL)) { - assert(fnfi[itg]); + fes1->GetElementVDofs(tr1->Elem1No, vdofs1); + fes2->GetElementVDofs(tr2->Elem1No, vdofs2); - fnfi[itg]->AssembleInterfaceVector(*fe1, *fe2, *tr1, *tr2, el_x1, el_x2, el_y1, el_y2); + x1.GetSubVector(vdofs1, el_x1); + x2.GetSubVector(vdofs2, el_x2); - y1.AddElementVector(vdofs1, el_y1); - y2.AddElementVector(vdofs2, el_y2); - } - } // if ((tr1 != NULL) && (tr2 != NULL)) - } // for (int bn = 0; bn < interface_infos.Size(); bn++) -} + // Both domains will have the adjacent element as Elem1. + fe1 = fes1->GetFE(tr1->Elem1No); + fe2 = fes2->GetFE(tr2->Elem1No); -void InterfaceForm::AssembleInterfaceGrad(Mesh *mesh1, Mesh *mesh2, - FiniteElementSpace *fes1, FiniteElementSpace *fes2, Array *interface_infos, - const Vector &x1, const Vector &x2, Array2D &mats) const -{ - assert(x1.Size() == fes1->GetTrueVSize()); - assert(x2.Size() == fes2->GetTrueVSize()); - - Array2D elemmats(2, 2); - for (int i = 0; i < elemmats.NumRows(); i++) - for (int j = 0; j < elemmats.NumCols(); j++) elemmats(i, j) = new DenseMatrix; - FaceElementTransformations *tr1, *tr2; - const FiniteElement *fe1, *fe2; - Array vdofs1, vdofs2; - Vector el_x1, el_x2, el_y1, el_y2; - - for (int bn = 0; bn < interface_infos->Size(); bn++) - { - InterfaceInfo *if_info = &((*interface_infos)[bn]); + for (int itg = 0; itg < fnfi.Size(); itg++) + { + assert(fnfi[itg]); - topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); + fnfi[itg]->AssembleInterfaceVector(*fe1, *fe2, *tr1, *tr2, el_x1, el_x2, el_y1, el_y2); - if ((tr1 != NULL) && (tr2 != NULL)) - { - fes1->GetElementVDofs(tr1->Elem1No, vdofs1); - fes2->GetElementVDofs(tr2->Elem1No, vdofs2); + y1.AddElementVector(vdofs1, el_y1); + y2.AddElementVector(vdofs2, el_y2); + } + } // if ((tr1 != NULL) && (tr2 != NULL)) + } // for (int bn = 0; bn < interface_infos.Size(); bn++) + } - x1.GetSubVector(vdofs1, el_x1); - x2.GetSubVector(vdofs2, el_x2); + void InterfaceForm::AssembleInterfaceGrad(Mesh *mesh1, Mesh *mesh2, + FiniteElementSpace *fes1, FiniteElementSpace *fes2, Array *interface_infos, + const Vector &x1, const Vector &x2, Array2D &mats) const + { + assert(x1.Size() == fes1->GetTrueVSize()); + assert(x2.Size() == fes2->GetTrueVSize()); + + Array2D elemmats(2, 2); + for (int i = 0; i < elemmats.NumRows(); i++) + for (int j = 0; j < elemmats.NumCols(); j++) + elemmats(i, j) = new DenseMatrix; + FaceElementTransformations *tr1, *tr2; + const FiniteElement *fe1, *fe2; + Array vdofs1, vdofs2; + Vector el_x1, el_x2, el_y1, el_y2; + + for (int bn = 0; bn < interface_infos->Size(); bn++) + { + InterfaceInfo *if_info = &((*interface_infos)[bn]); - // Both domains will have the adjacent element as Elem1. - fe1 = fes1->GetFE(tr1->Elem1No); - fe2 = fes2->GetFE(tr2->Elem1No); + topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); - for (int itg = 0; itg < fnfi.Size(); itg++) + if ((tr1 != NULL) && (tr2 != NULL)) { - assert(fnfi[itg]); + fes1->GetElementVDofs(tr1->Elem1No, vdofs1); + fes2->GetElementVDofs(tr2->Elem1No, vdofs2); - fnfi[itg]->AssembleInterfaceGrad(*fe1, *fe2, *tr1, *tr2, el_x1, el_x2, elemmats); + x1.GetSubVector(vdofs1, el_x1); + x2.GetSubVector(vdofs2, el_x2); - mats(0, 0)->AddSubMatrix(vdofs1, vdofs1, *elemmats(0, 0), skip_zeros); - mats(0, 1)->AddSubMatrix(vdofs1, vdofs2, *elemmats(0, 1), skip_zeros); - mats(1, 0)->AddSubMatrix(vdofs2, vdofs1, *elemmats(1, 0), skip_zeros); - mats(1, 1)->AddSubMatrix(vdofs2, vdofs2, *elemmats(1, 1), skip_zeros); - } - } // if ((tr1 != NULL) && (tr2 != NULL)) - } // for (int bn = 0; bn < interface_infos.Size(); bn++) + // Both domains will have the adjacent element as Elem1. + fe1 = fes1->GetFE(tr1->Elem1No); + fe2 = fes2->GetFE(tr2->Elem1No); - DeletePointers(elemmats); -} + for (int itg = 0; itg < fnfi.Size(); itg++) + { + assert(fnfi[itg]); -/* - MixedInterfaceForm -*/ + fnfi[itg]->AssembleInterfaceGrad(*fe1, *fe2, *tr1, *tr2, el_x1, el_x2, elemmats); -MixedInterfaceForm::MixedInterfaceForm( - Array &meshes_, Array &trial_fes_, - Array &test_fes_, TopologyHandler *topol_) - : meshes(meshes_), trial_fes(trial_fes_), test_fes(test_fes_), topol_handler(topol_), numSub(meshes_.Size()) -{ - assert(trial_fes_.Size() == numSub); - assert(test_fes_.Size() == numSub); - - trial_block_offsets.SetSize(numSub + 1); - test_block_offsets.SetSize(numSub + 1); - trial_block_offsets = 0; - test_block_offsets = 0; - for (int i = 1; i < numSub + 1; i++) - { - trial_block_offsets[i] = trial_fes[i-1]->GetTrueVSize(); - test_block_offsets[i] = test_fes[i-1]->GetTrueVSize(); + mats(0, 0)->AddSubMatrix(vdofs1, vdofs1, *elemmats(0, 0), skip_zeros); + mats(0, 1)->AddSubMatrix(vdofs1, vdofs2, *elemmats(0, 1), skip_zeros); + mats(1, 0)->AddSubMatrix(vdofs2, vdofs1, *elemmats(1, 0), skip_zeros); + mats(1, 1)->AddSubMatrix(vdofs2, vdofs2, *elemmats(1, 1), skip_zeros); + } + } // if ((tr1 != NULL) && (tr2 != NULL)) + } // for (int bn = 0; bn < interface_infos.Size(); bn++) + + DeletePointers(elemmats); } - trial_block_offsets.PartialSum(); - test_block_offsets.PartialSum(); -} -MixedInterfaceForm::~MixedInterfaceForm() -{ - DeletePointers(fnfi); -} + /* + MixedInterfaceForm + */ -void MixedInterfaceForm::AssembleInterfaceMatrixes(Array2D &mats) const -{ - const PortInfo *pInfo; - Array midx(2); - Array2D b_mats_p(2,2); - Mesh *mesh1, *mesh2; - FiniteElementSpace *trial_fes1, *trial_fes2, *test_fes1, *test_fes2; - Array* interface_infos; - - for (int p = 0; p < topol_handler->GetNumPorts(); p++) + MixedInterfaceForm::MixedInterfaceForm( + Array &meshes_, Array &trial_fes_, + Array &test_fes_, TopologyHandler *topol_) + : meshes(meshes_), trial_fes(trial_fes_), test_fes(test_fes_), topol_handler(topol_), numSub(meshes_.Size()) { - pInfo = topol_handler->GetPortInfo(p); + assert(trial_fes_.Size() == numSub); + assert(test_fes_.Size() == numSub); + + trial_block_offsets.SetSize(numSub + 1); + test_block_offsets.SetSize(numSub + 1); + trial_block_offsets = 0; + test_block_offsets = 0; + for (int i = 1; i < numSub + 1; i++) + { + trial_block_offsets[i] = trial_fes[i - 1]->GetTrueVSize(); + test_block_offsets[i] = test_fes[i - 1]->GetTrueVSize(); + } + trial_block_offsets.PartialSum(); + test_block_offsets.PartialSum(); + } - midx[0] = pInfo->Mesh1; - midx[1] = pInfo->Mesh2; - - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) - b_mats_p(i, j) = mats(midx[i], midx[j]); + MixedInterfaceForm::~MixedInterfaceForm() + { + DeletePointers(fnfi); + } - mesh1 = meshes[midx[0]]; - mesh2 = meshes[midx[1]]; + void MixedInterfaceForm::AssembleInterfaceMatrixes(Array2D &mats) const + { + const PortInfo *pInfo; + Array midx(2); + Array2D b_mats_p(2, 2); + Mesh *mesh1, *mesh2; + FiniteElementSpace *trial_fes1, *trial_fes2, *test_fes1, *test_fes2; + Array *interface_infos; + + for (int p = 0; p < topol_handler->GetNumPorts(); p++) + { + pInfo = topol_handler->GetPortInfo(p); - trial_fes1 = trial_fes[midx[0]]; - trial_fes2 = trial_fes[midx[1]]; - test_fes1 = test_fes[midx[0]]; - test_fes2 = test_fes[midx[1]]; + midx[0] = pInfo->Mesh1; + midx[1] = pInfo->Mesh2; - interface_infos = topol_handler->GetInterfaceInfos(p); + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + b_mats_p(i, j) = mats(midx[i], midx[j]); - AssembleInterfaceMatrix( - mesh1, mesh2, trial_fes1, trial_fes2, test_fes1, test_fes2, interface_infos, b_mats_p); - } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) -} + mesh1 = meshes[midx[0]]; + mesh2 = meshes[midx[1]]; -void MixedInterfaceForm::AssembleInterfaceMatrix( - Mesh *mesh1, Mesh *mesh2, FiniteElementSpace *trial_fes1, FiniteElementSpace *trial_fes2, - FiniteElementSpace *test_fes1, FiniteElementSpace *test_fes2, Array *interface_infos, - Array2D &mats) const -{ - Array2D elemmats(2, 2); - for (int i = 0; i < elemmats.NumRows(); i++) - for (int j = 0; j < elemmats.NumCols(); j++) elemmats(i, j) = new DenseMatrix; + trial_fes1 = trial_fes[midx[0]]; + trial_fes2 = trial_fes[midx[1]]; + test_fes1 = test_fes[midx[0]]; + test_fes2 = test_fes[midx[1]]; - FaceElementTransformations *tr1, *tr2; - const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2; - Array *> test_vdofs(2), trial_vdofs(2); - trial_vdofs[0] = new Array; - trial_vdofs[1] = new Array; - test_vdofs[0] = new Array; - test_vdofs[1] = new Array; + interface_infos = topol_handler->GetInterfaceInfos(p); - InterfaceInfo *if_info; + AssembleInterfaceMatrix( + mesh1, mesh2, trial_fes1, trial_fes2, test_fes1, test_fes2, interface_infos, b_mats_p); + } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) + } - for (int bn = 0; bn < interface_infos->Size(); bn++) + void MixedInterfaceForm::AssembleInterfaceMatrix( + Mesh *mesh1, Mesh *mesh2, FiniteElementSpace *trial_fes1, FiniteElementSpace *trial_fes2, + FiniteElementSpace *test_fes1, FiniteElementSpace *test_fes2, Array *interface_infos, + Array2D &mats) const { - if_info = &((*interface_infos)[bn]); - - topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); - - if ((tr1 != NULL) && (tr2 != NULL)) + Array2D elemmats(2, 2); + for (int i = 0; i < elemmats.NumRows(); i++) + for (int j = 0; j < elemmats.NumCols(); j++) + elemmats(i, j) = new DenseMatrix; + + FaceElementTransformations *tr1, *tr2; + const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2; + Array *> test_vdofs(2), trial_vdofs(2); + trial_vdofs[0] = new Array; + trial_vdofs[1] = new Array; + test_vdofs[0] = new Array; + test_vdofs[1] = new Array; + + InterfaceInfo *if_info; + + for (int bn = 0; bn < interface_infos->Size(); bn++) { - trial_fes1->GetElementVDofs(tr1->Elem1No, *trial_vdofs[0]); - trial_fes2->GetElementVDofs(tr2->Elem1No, *trial_vdofs[1]); - test_fes1->GetElementVDofs(tr1->Elem1No, *test_vdofs[0]); - test_fes2->GetElementVDofs(tr2->Elem1No, *test_vdofs[1]); - // Both domains will have the adjacent element as Elem1. - trial_fe1 = trial_fes1->GetFE(tr1->Elem1No); - trial_fe2 = trial_fes2->GetFE(tr2->Elem1No); - test_fe1 = test_fes1->GetFE(tr1->Elem1No); - test_fe2 = test_fes2->GetFE(tr2->Elem1No); - - for (int itg = 0; itg < fnfi.Size(); itg++) - { - assert(fnfi[itg]); + if_info = &((*interface_infos)[bn]); - fnfi[itg]->AssembleInterfaceMatrix( - *trial_fe1, *trial_fe2, *test_fe1, *test_fe2, *tr1, *tr2, elemmats); + topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); - for (int i = 0; i < 2; i++) { - for (int j = 0; j < 2; j++) { - mats(i, j)->AddSubMatrix(*test_vdofs[i], *trial_vdofs[j], *elemmats(i,j), skip_zeros); + if ((tr1 != NULL) && (tr2 != NULL)) + { + trial_fes1->GetElementVDofs(tr1->Elem1No, *trial_vdofs[0]); + trial_fes2->GetElementVDofs(tr2->Elem1No, *trial_vdofs[1]); + test_fes1->GetElementVDofs(tr1->Elem1No, *test_vdofs[0]); + test_fes2->GetElementVDofs(tr2->Elem1No, *test_vdofs[1]); + // Both domains will have the adjacent element as Elem1. + trial_fe1 = trial_fes1->GetFE(tr1->Elem1No); + trial_fe2 = trial_fes2->GetFE(tr2->Elem1No); + test_fe1 = test_fes1->GetFE(tr1->Elem1No); + test_fe2 = test_fes2->GetFE(tr2->Elem1No); + + for (int itg = 0; itg < fnfi.Size(); itg++) + { + assert(fnfi[itg]); + + fnfi[itg]->AssembleInterfaceMatrix( + *trial_fe1, *trial_fe2, *test_fe1, *test_fe2, *tr1, *tr2, elemmats); + + for (int i = 0; i < 2; i++) + { + for (int j = 0; j < 2; j++) + { + mats(i, j)->AddSubMatrix(*test_vdofs[i], *trial_vdofs[j], *elemmats(i, j), skip_zeros); + } } - } - } // for (int itg = 0; itg < fnfi.Size(); itg++) - } // if ((tr1 != NULL) && (tr2 != NULL)) - } // for (int bn = 0; bn < interface_infos.Size(); bn++) + } // for (int itg = 0; itg < fnfi.Size(); itg++) + } // if ((tr1 != NULL) && (tr2 != NULL)) + } // for (int bn = 0; bn < interface_infos.Size(); bn++) - DeletePointers(test_vdofs); - DeletePointers(trial_vdofs); - DeletePointers(elemmats); -} + DeletePointers(test_vdofs); + DeletePointers(trial_vdofs); + DeletePointers(elemmats); + } -void MixedInterfaceForm::AssembleInterfaceMatrixAtPort( - const int p, Array &trial_fes_comp, - Array &test_fes_comp, Array2D &mats_p) const -{ - const int num_ref_ports = topol_handler->GetNumRefPorts(); + void MixedInterfaceForm::AssembleInterfaceMatrixAtPort( + const int p, Array &trial_fes_comp, + Array &test_fes_comp, Array2D &mats_p) const + { + const int num_ref_ports = topol_handler->GetNumRefPorts(); - assert(topol_handler->GetType() == TopologyHandlerMode::COMPONENT); - assert((p >= 0) && (p < num_ref_ports)); - DeletePointers(mats_p); - mats_p.SetSize(2, 2); + assert(topol_handler->GetType() == TopologyHandlerMode::COMPONENT); + assert((p >= 0) && (p < num_ref_ports)); + DeletePointers(mats_p); + mats_p.SetSize(2, 2); - const int num_comp = topol_handler->GetNumComponents(); - assert(trial_fes_comp.Size() == num_comp); - assert(test_fes_comp.Size() == num_comp); + const int num_comp = topol_handler->GetNumComponents(); + assert(trial_fes_comp.Size() == num_comp); + assert(test_fes_comp.Size() == num_comp); - int c1, c2; - topol_handler->GetComponentPair(p, c1, c2); - Mesh *comp1 = topol_handler->GetComponentMesh(c1); - Mesh *comp2 = topol_handler->GetComponentMesh(c2); + int c1, c2; + topol_handler->GetComponentPair(p, c1, c2); + Mesh *comp1 = topol_handler->GetComponentMesh(c1); + Mesh *comp2 = topol_handler->GetComponentMesh(c2); - // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. - // Need to use two copied instances. - Mesh mesh1(*comp1); - Mesh mesh2(*comp2); + // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. + // Need to use two copied instances. + Mesh mesh1(*comp1); + Mesh mesh2(*comp2); - Array c_idx(2); - c_idx[0] = c1; c_idx[1] = c2; + Array c_idx(2); + c_idx[0] = c1; + c_idx[1] = c2; - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) - mats_p(i, j) = new SparseMatrix(test_fes_comp[c_idx[i]]->GetTrueVSize(), trial_fes_comp[c_idx[j]]->GetTrueVSize()); + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + mats_p(i, j) = new SparseMatrix(test_fes_comp[c_idx[i]]->GetTrueVSize(), trial_fes_comp[c_idx[j]]->GetTrueVSize()); - Array* const if_infos = topol_handler->GetRefInterfaceInfos(p); + Array *const if_infos = topol_handler->GetRefInterfaceInfos(p); - // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. - // Need to use two copied instances. - AssembleInterfaceMatrix( - &mesh1, &mesh2, trial_fes_comp[c_idx[0]], trial_fes_comp[c_idx[1]], - test_fes_comp[c_idx[0]], test_fes_comp[c_idx[1]], if_infos, mats_p); + // NOTE: If comp1 == comp2, using comp1 and comp2 directly leads to an incorrect penalty matrix. + // Need to use two copied instances. + AssembleInterfaceMatrix( + &mesh1, &mesh2, trial_fes_comp[c_idx[0]], trial_fes_comp[c_idx[1]], + test_fes_comp[c_idx[0]], test_fes_comp[c_idx[1]], if_infos, mats_p); - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) - mats_p(i, j)->Finalize(); -} + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + mats_p(i, j)->Finalize(); + } } diff --git a/src/interfaceinteg.cpp b/src/interfaceinteg.cpp index 57b8b3d5..73bc4a72 100644 --- a/src/interfaceinteg.cpp +++ b/src/interfaceinteg.cpp @@ -1008,55 +1008,54 @@ namespace mfem FaceElementTransformations &Trans1, FaceElementTransformations &Trans2, Array2D &elmats) { -#ifdef MFEM_THREAD_SAFE - // For descriptions of these variables, see the class declaration. - Vector shape1, shape2; - DenseMatrix dshape1, dshape2; - DenseMatrix adjJ; - DenseMatrix dshape1_ps, dshape2_ps; - Vector nor1, nor2; - Vector nL1, nL2; - Vector nM1, nM2; - Vector dshape1_dnM, dshape2_dnM; - DenseMatrix jmat; -#endif + /* #ifdef MFEM_THREAD_SAFE + // For descriptions of these variables, see the class declaration. + Vector shape1, shape2; + DenseMatrix dshape1, dshape2; + DenseMatrix adjJ; + DenseMatrix dshape1_ps, dshape2_ps; + Vector nor1, nor2; + Vector nL1, nL2; + Vector nM1, nM2; + Vector dshape1_dnM, dshape2_dnM; + DenseMatrix jmat; + #endif */ bool boundary = false; const int dim = el1.GetDim(); const int ndofs1 = el1.GetDof(); const int ndofs2 = (boundary) ? 0 : el2.GetDof(); - const int nvdofs = dim * (ndofs1 + ndofs2); - + // const int nvdofs = dim * (ndofs1 + ndofs2); + const int nvdofs1 = dim * ndofs1; + const int nvdofs2 = dim * ndofs2; // Initially 'elmat' corresponds to the term: // < { sigma(u) . n }, [v] > = // < { (lambda div(u) I + mu (grad(u) + grad(u)^T)) . n }, [v] > // But eventually, it's going to be replaced by: // elmat := -elmat + alpha*elmat^T + jmat - elmats(0, 0)->SetSize(ndofs1, ndofs1); - elmats(0, 1)->SetSize(ndofs1, ndofs2); - elmats(1, 0)->SetSize(ndofs2, ndofs1); - elmats(1, 1)->SetSize(ndofs2, ndofs2); + elmats(0, 0)->SetSize(nvdofs1, nvdofs1); + elmats(0, 1)->SetSize(nvdofs1, nvdofs2); + elmats(1, 0)->SetSize(nvdofs2, nvdofs1); + elmats(1, 1)->SetSize(nvdofs2, nvdofs2); for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) *elmats(i, j) = 0.0; // elmat.SetSize(ndofs); // elmat = 0.0; - const bool kappa_is_nonzero = (kappa != 0.0); + jmats.SetSize(2, 2); if (kappa_is_nonzero) { - jmats.SetSize(2, 2); - jmats(0, 0) = new DenseMatrix(ndofs1, ndofs1); + jmats(0, 0) = new DenseMatrix(nvdofs1, nvdofs1); // only the lower-triangular part of jmat is assembled. jmats(0, 1) = NULL; - jmats(1, 0) = new DenseMatrix(ndofs2, ndofs1); - jmats(1, 1) = new DenseMatrix(ndofs2, ndofs2); + jmats(1, 0) = new DenseMatrix(nvdofs2, nvdofs1); + jmats(1, 1) = new DenseMatrix(nvdofs2, nvdofs2); for (int i = 0; i < 2; i++) for (int j = 0; j <= i; j++) *jmats(i, j) = 0.0; // jmat.SetSize(ndofs); // jmat = 0.; } - adjJ.SetSize(dim); shape1.SetSize(ndofs1); dshape1.SetSize(ndofs1, dim); @@ -1076,7 +1075,6 @@ namespace mfem nM2.SetSize(dim); dshape2_dnM.SetSize(ndofs2); } - const IntegrationRule *ir = IntRule; if (ir == NULL) { @@ -1085,7 +1083,6 @@ namespace mfem ir = &IntRules.Get(Trans1.GetGeometryType(), order); assert(Trans1.GetGeometryType() == Trans2.GetGeometryType()); } - for (int pind = 0; pind < ir->GetNPoints(); ++pind) { const IntegrationPoint &ip = ir->IntPoint(pind); @@ -1098,8 +1095,7 @@ namespace mfem // Note: eip1 and eip2 come from Element1 of Trans1 and Trans2 respectively. const IntegrationPoint &eip1 = Trans1.GetElement1IntPoint(); const IntegrationPoint &eip2 = Trans2.GetElement1IntPoint(); - - // computing outward normal vectors. + // computing outward normal vectors. if (dim == 1) { nor1(0) = 2 * eip1.x - 1.0; @@ -1166,7 +1162,7 @@ namespace mfem // (1,2) block AssembleBlock( - dim, ndofs1, ndofs2, 0, dim * ndofs1, jmatcoef, nL2, nM2, + dim, ndofs1, ndofs2, 0, dim * ndofs1, 0.0, nL2, nM2, shape1, shape2, dshape2_dnM, dshape2_ps, *elmats(0, 1), *jmats(0, 1)); // (2,1) block AssembleBlock( @@ -1256,21 +1252,15 @@ namespace mfem const Vector &col_dshape_dnM, const DenseMatrix &col_dshape, DenseMatrix &elmat, DenseMatrix &jmat) { // row_offset and col_offset are not needed for elmat. - for (int jm = 0, j = 0; jm < dim; ++jm) + for (int d = 0; d < dim; ++d) { + int j = d * col_ndofs; for (int jdof = 0; jdof < col_ndofs; ++jdof, ++j) { + int i = d * row_ndofs; const double t2 = col_dshape_dnM(jdof); - for (int im = 0, i = 0; im < dim; ++im) - { - const double t1 = col_dshape(jdof, jm) * col_nL(im); - const double t3 = col_dshape(jdof, im) * col_nM(jm); - const double tt = t1 + ((im == jm) ? t2 : 0.0) + t3; - for (int idof = 0; idof < row_ndofs; ++idof, ++i) - { - elmat(i, j) += row_shape(idof) * tt; - } - } + for (int idof = 0; idof < row_ndofs; ++idof, ++i) + elmat(i, j) += row_shape(idof) * t2; } } diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index c7fbcb3d..49104abf 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -59,6 +59,14 @@ LinElastSolver::~LinElastSolver() delete mumps; } +void LinElastSolver::SetupBCVariables() +{ + MultiBlockSolver::SetupBCVariables(); + + bdr_coeffs.SetSize(numBdr); + bdr_coeffs = NULL; +} + void LinElastSolver::SetupMaterialVariables() { int max_bdr_attr = -1; // A bit redundant... @@ -71,20 +79,22 @@ void LinElastSolver::SetupMaterialVariables() Vector lambda(max_bdr_attr); lambda = 1.0; // Set lambda = 1 for all element attributes. lambda(0) = 50.0; // Set lambda = 50 for element attribute 1. - PWConstCoefficient lambda_c(lambda); + // PWConstCoefficient lambda_c(lambda); Vector mu(max_bdr_attr); mu = 1.0; // Set mu = 1 for all element attributes. mu(0) = 50.0; // Set mu = 50 for element attribute 1. - PWConstCoefficient mu_c(mu); + // PWConstCoefficient mu_c(mu); lambda_cs.SetSize(numSub); mu_cs.SetSize(numSub); for (int m = 0; m < numSub; m++) { - lambda_cs[m] = &lambda_c; - mu_cs[m] = &mu_c; + // lambda_cs[m] = &lambda_c; + // mu_cs[m] = &mu_c; + lambda_cs[m] = new PWConstCoefficient(lambda); + mu_cs[m] = new PWConstCoefficient(mu); } } @@ -149,35 +159,65 @@ void LinElastSolver::BuildOperators() void LinElastSolver::BuildRHSOperators() { bs.SetSize(numSub); - for (int m = 0; m < numSub; m++) { bs[m] = new LinearForm(fes[m], RHS->GetBlock(m).GetData()); - bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator(*init_x, *(lambda_cs[m]), *(mu_cs[m]), alpha, kappa), *(bdr_markers[0])); // m ight be wrong TODO bdr_marker will be extended to include forces too, also, unsure if this is correct + } +} + +void LinElastSolver::SetupRHSBCOperators() +{ + assert(bs.Size() == numSub); + for (int m = 0; m < numSub; m++) + { + assert(bs[m]); + for (int b = 0; b < global_bdr_attributes.Size(); b++) + { + int idx = meshes[m]->bdr_attributes.Find(global_bdr_attributes[b]); + if (idx < 0) + continue; + if (!BCExistsOnBdr(b)) + continue; + + // bs[m]->AddBdrFaceIntegrator(new DGDirichletLFIntegrator(*bdr_coeffs[b], sigma, kappa), *bdr_markers[b]); + // bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator(*init_x, *lambda_cs[b], *mu_cs[b], alpha, kappa), *bdr_markers[b]); + bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator(*bdr_coeffs[b], *lambda_cs[b], *mu_cs[b], alpha, kappa), *bdr_markers[b]); + } } } void LinElastSolver::BuildDomainOperators() { // SanityCheckOnCoeffs(); - as.SetSize(numSub); for (int m = 0; m < numSub; m++) { as[m] = new BilinearForm(fes[m]); as[m]->AddDomainIntegrator(new ElasticityIntegrator(*(lambda_cs[m]), *(mu_cs[m]))); + if (full_dg) { as[m]->AddInteriorFaceIntegrator( new DGElasticityIntegrator(*(lambda_cs[m]), *(mu_cs[m]), alpha, kappa)); - as[m]->AddBdrFaceIntegrator( - new DGElasticityIntegrator(*(lambda_cs[m]), *(mu_cs[m]), alpha, kappa), *(bdr_markers[0])); + + for (int b = 0; b < global_bdr_attributes.Size(); b++) + { + int idx = meshes[m]->bdr_attributes.Find(global_bdr_attributes[b]); + if (idx < 0) + continue; + if (!BCExistsOnBdr(b)) + continue; + + as[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator(*bdr_coeffs[b], *lambda_cs[b], *mu_cs[b], alpha, kappa), *bdr_markers[b]); + } } + as[m]->Assemble(); + as[m]->Finalize(); } a_itf = new InterfaceForm(meshes, fes, topol_handler); - a_itf->AddIntefaceIntegrator(new InterfaceDGElasticityIntegrator(alpha, kappa)); + a_itf->AddIntefaceIntegrator(new InterfaceDGElasticityIntegrator(lambda_cs[0], mu_cs[0], alpha, kappa)); } void LinElastSolver::Assemble() @@ -189,13 +229,13 @@ void LinElastSolver::Assemble() void LinElastSolver::AssembleRHS() { // SanityCheckOnCoeffs(); - MFEM_ASSERT(bs.Size() == numSub, "LinearForm bs != numSub.\n"); - for (int m = 0; m < numSub; m++) { MFEM_ASSERT(bs[m], "LinearForm or BilinearForm pointer of a subdomain is not associated!\n"); + cout << bs[m]->GetFLFI()[0] << endl; bs[m]->Assemble(); + cout << "linear form norm: " << bs[m]->Norml2() << endl; } for (int m = 0; m < numSub; m++) @@ -206,15 +246,16 @@ void LinElastSolver::AssembleRHS() void LinElastSolver::AssembleOperator() { // SanityCheckOnCoeffs(); - MFEM_ASSERT(as.Size() == numSub, "BilinearForm bs != numSub.\n"); - for (int m = 0; m < numSub; m++) { MFEM_ASSERT(as[m], "LinearForm or BilinearForm pointer of a subdomain is not associated!\n"); as[m]->Assemble(); - } + as[m]->Finalize(); + double binorm = as[m]->SpMat().ToDenseMatrix()->FNorm(); + cout << "bilinear form norm: " << binorm << endl; + } mats.SetSize(numSub, numSub); for (int i = 0; i < numSub; i++) { @@ -231,7 +272,6 @@ void LinElastSolver::AssembleOperator() } } AssembleInterfaceMatrixes(); - for (int m = 0; m < numSub; m++) as[m]->Finalize(); @@ -282,8 +322,13 @@ bool LinElastSolver::Solve() double rtol = config.GetOption("solver/relative_tolerance", 1.e-15); double atol = config.GetOption("solver/absolute_tolerance", 1.e-15); int print_level = config.GetOption("solver/print_level", 0); - + for (size_t i = 0; i < U->NumBlocks(); i++) + { + cout << "Unorm " << i << ": " << U->GetBlock(i).Norml2() << endl; + } + cout << "RHSnorm: " << RHS->Norml2() << endl; // TODO: need to change when the actual parallelization is implemented. + cout << "direct_solve is: " << direct_solve << endl; if (direct_solve) { assert(mumps); @@ -351,21 +396,97 @@ bool LinElastSolver::Solve() return converged; } -void LinElastSolver::SetupBCVariables(){"LinElastSolver::SetupBCVariables is not implemented yet!\n";} -void LinElastSolver::AddBCFunction(std::function F, const int battr){"LinElastSolver::AddBCFunction is not implemented yet!\n";} -void LinElastSolver::AddBCFunction(const double &F, const int battr){"LinElastSolver::AddBCFunction is not implemented yet!\n";} -bool LinElastSolver::BCExistsOnBdr(const int &global_battr_idx){std::cout<<"LinElastSolver::BCExistsOnBdr is not implemented yet!\n"; return false;} -void LinElastSolver::SetupBCOperators(){"LinElastSolver::SetupBCOperators is not implemented yet!\n";} -void LinElastSolver::SetupRHSBCOperators(){"LinElastSolver::SetupRHSBCOperators is not implemented yet!\n";} -void LinElastSolver::SetupDomainBCOperators(){"LinElastSolver::SetupDomainBCOperators is not implemented yet!\n";} +void PrintVector(string filename, Vector &vec) +{ + std::ofstream outfile(filename); + double tol = 1e-7; + double val = 0.0; + for (size_t i = 0; i < vec.Size(); i++) + { + val = vec[i]; + if (abs(val) < tol) + { + val = 0.0; + } + + outfile << setprecision(8) << val << endl; + } + outfile.close(); + cout << "done printing vector" << endl; +} + +void PrintMatrix(string filename, DenseMatrix &mat) +{ + std::ofstream outfile(filename); + + double tol = 1e-7; + double val = 0.0; + + for (size_t i = 0; i < mat.Height(); i++) + { + for (size_t j = 0; j < mat.Width(); j++) + { + val = mat(i, j); + if (abs(val) < tol) + { + val = 0.0; + } + + outfile << setprecision(8) << val << " "; + } + outfile << endl; + } + outfile.close(); + cout << "done printing matrix" << endl; +} + +void PrintBlockVector(string filename, BlockVector &bvec) +{ + std::ofstream outfile(filename); + + for (size_t i = 0; i < bvec.GetBlock(0).Size(); i++) + { + for (size_t j = 0; j < bvec.NumBlocks(); j++) + { + outfile << setprecision(1) << bvec.GetBlock(j)[i] << " "; + } + outfile << endl; + } + outfile.close(); + cout << "done printing blockvector" << endl; +} + +void LinElastSolver::PrintOperators() +{ + PrintMatrix("scaleuprom_a.txt", *(as[0]->SpMat().ToDenseMatrix())); + PrintVector("scaleuprom_b.txt", *bs[0]); +} + +void LinElastSolver::SetupBCVariables() { "LinElastSolver::SetupBCVariables is not implemented yet!\n"; } +void LinElastSolver::AddBCFunction(std::function F, const int battr) { "LinElastSolver::AddBCFunction is not implemented yet!\n"; } +void LinElastSolver::AddBCFunction(const double &F, const int battr) { "LinElastSolver::AddBCFunction is not implemented yet!\n"; } +bool LinElastSolver::BCExistsOnBdr(const int &global_battr_idx) +{ + std::cout << "LinElastSolver::BCExistsOnBdr is not implemented yet!\n"; + return false; +} + +void LinElastSolver::SetupBCOperators() +{ + SetupRHSBCOperators(); + SetupDomainBCOperators(); +} + +void LinElastSolver::SetupRHSBCOperators() { "LinElastSolver::SetupRHSBCOperators is not implemented yet!\n"; } +void LinElastSolver::SetupDomainBCOperators() { "LinElastSolver::SetupDomainBCOperators is not implemented yet!\n"; } // Component-wise assembly -void LinElastSolver::BuildCompROMElement(Array &fes_comp){"LinElastSolver::BuildCompROMElement is not implemented yet!\n";} -void LinElastSolver::BuildBdrROMElement(Array &fes_comp){"LinElastSolver::BuildBdrROMElement is not implemented yet!\n";} -void LinElastSolver::BuildInterfaceROMElement(Array &fes_comp){"LinElastSolver::BuildInterfaceROMElement is not implemented yet!\n";} +void LinElastSolver::BuildCompROMElement(Array &fes_comp) { "LinElastSolver::BuildCompROMElement is not implemented yet!\n"; } +void LinElastSolver::BuildBdrROMElement(Array &fes_comp) { "LinElastSolver::BuildBdrROMElement is not implemented yet!\n"; } +void LinElastSolver::BuildInterfaceROMElement(Array &fes_comp) { "LinElastSolver::BuildInterfaceROMElement is not implemented yet!\n"; } -void LinElastSolver::ProjectOperatorOnReducedBasis(){"LinElastSolver::ProjectOperatorOnReducedBasis is not implemented yet!\n";} +void LinElastSolver::ProjectOperatorOnReducedBasis() { "LinElastSolver::ProjectOperatorOnReducedBasis is not implemented yet!\n"; } -void LinElastSolver::SanityCheckOnCoeffs(){"LinElastSolver::SanityCheckOnCoeffs is not implemented yet!\n";} +void LinElastSolver::SanityCheckOnCoeffs() { "LinElastSolver::SanityCheckOnCoeffs is not implemented yet!\n"; } -void LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem){"LinElastSolver::SetParameterizedProblem is not implemented yet!\n";} +void LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem) { "LinElastSolver::SetParameterizedProblem is not implemented yet!\n"; } diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index 433d93f5..14a7baa6 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -47,14 +47,26 @@ void RunExample() test.SaveVisualization(); } -MultiBlockSolver* InitSolver() +MultiBlockSolver *InitSolver() { std::string solver_type = config.GetRequiredOption("main/solver"); MultiBlockSolver *solver = NULL; - if (solver_type == "poisson") { solver = new PoissonSolver; } - else if (solver_type == "stokes") { solver = new StokesSolver; } - else if (solver_type == "steady-ns") { solver = new SteadyNSSolver; } - else if (solver_type == "linelast") { solver = new LinElastSolver; } + if (solver_type == "poisson") + { + solver = new PoissonSolver; + } + else if (solver_type == "stokes") + { + solver = new StokesSolver; + } + else if (solver_type == "steady-ns") + { + solver = new SteadyNSSolver; + } + else if (solver_type == "linelast") + { + solver = new LinElastSolver; + } else { printf("Unknown MultiBlockSolver %s!\n", solver_type.c_str()); @@ -64,9 +76,9 @@ MultiBlockSolver* InitSolver() return solver; } -SampleGenerator* InitSampleGenerator(MPI_Comm comm) +SampleGenerator *InitSampleGenerator(MPI_Comm comm) { - SampleGenerator* generator = NULL; + SampleGenerator *generator = NULL; std::string type = config.GetOption("sample_generation/type", "base"); @@ -94,8 +106,10 @@ std::vector GetGlobalBasisTagList(const TopologyHandlerMode &topol_ if (train_mode == TrainMode::INDIVIDUAL) { TopologyHandler *topol_handler; - if (topol_mode == TopologyHandlerMode::SUBMESH) topol_handler = new SubMeshTopologyHandler(); - else if (topol_mode == TopologyHandlerMode::COMPONENT) topol_handler = new ComponentTopologyHandler(); + if (topol_mode == TopologyHandlerMode::SUBMESH) + topol_handler = new SubMeshTopologyHandler(); + else if (topol_mode == TopologyHandlerMode::COMPONENT) + topol_handler = new ComponentTopologyHandler(); else mfem_error("GetGlobalBasisTagList - TopologyHandlerMode is not set!\n"); @@ -104,7 +118,7 @@ std::vector GetGlobalBasisTagList(const TopologyHandlerMode &topol_ delete topol_handler; } - else // if (train_mode == TrainMode::UNIVERSAL) + else // if (train_mode == TrainMode::UNIVERSAL) { if (topol_mode == TopologyHandlerMode::SUBMESH) { @@ -115,19 +129,22 @@ std::vector GetGlobalBasisTagList(const TopologyHandlerMode &topol_ YAML::Node component_dict = config.FindNode("mesh/component-wise/components"); assert(component_dict); for (int p = 0; p < component_dict.size(); p++) - component_list.push_back(config.GetRequiredOptionFromDict("name", component_dict[p])); + component_list.push_back(config.GetRequiredOptionFromDict("name", component_dict[p])); } else mfem_error("GetGlobalBasisTagList - TopologyHandlerMode is not set!\n"); - } // if (train_mode == TrainMode::UNIVERSAL) + } // if (train_mode == TrainMode::UNIVERSAL) std::vector var_list(0); if (separate_variable_basis) { std::string solver_type = config.GetRequiredOption("main/solver"); - if (solver_type == "poisson") var_list = PoissonSolver::GetVariableNames(); - else if (solver_type == "stokes") var_list = StokesSolver::GetVariableNames(); - else if (solver_type == "steady-ns") var_list = SteadyNSSolver::GetVariableNames(); + if (solver_type == "poisson") + var_list = PoissonSolver::GetVariableNames(); + else if (solver_type == "stokes") + var_list = StokesSolver::GetVariableNames(); + else if (solver_type == "steady-ns") + var_list = SteadyNSSolver::GetVariableNames(); else { printf("Unknown MultiBlockSolver %s!\n", solver_type.c_str()); @@ -160,7 +177,8 @@ void GenerateSamples(MPI_Comm comm) int s = 0; while (s < sample_generator->GetTotalSampleSize()) { - if (!sample_generator->IsMyJob(s)) continue; + if (!sample_generator->IsMyJob(s)) + continue; // NOTE: this will change config.dict_ sample_generator->SetSampleParams(s); @@ -194,7 +212,7 @@ void GenerateSamples(MPI_Comm comm) continue; } } - + test->SaveSolution(sol_file); test->SaveVisualization(); @@ -249,7 +267,7 @@ void TrainROM(MPI_Comm comm) { // Find if additional inputs are specified for basis_tags[p]. YAML::Node basis_tag_input = config.LookUpFromDict("name", basis_tags[p], basis_list); - + // If basis_tags[p] has additional inputs, parse them. if (basis_tag_input) { @@ -258,14 +276,14 @@ void TrainROM(MPI_Comm comm) // parse the sample snapshot file list. file_list = config.GetOptionFromDict>( - "snapshot_files", std::vector(0), basis_tag_input); + "snapshot_files", std::vector(0), basis_tag_input); YAML::Node snapshot_format = config.FindNodeFromDict("snapshot_format", basis_tag_input); if (snapshot_format) { FilenameParam snapshot_param("", snapshot_format); snapshot_param.ParseFilenames(file_list); } - } // if (basis_tag_input) + } // if (basis_tag_input) else num_basis = num_basis_default; } @@ -284,7 +302,7 @@ void TrainROM(MPI_Comm comm) } sample_generator->FormReducedBasis(basis_prefix, basis_tags[p], file_list, num_basis); - } // for (int p = 0; p < basis_tags.size(); p++) + } // for (int p = 0; p < basis_tags.size(); p++) delete sample_generator; @@ -301,10 +319,17 @@ void AuxiliaryTrainROM(MPI_Comm comm) { ParameterizedProblem *problem = InitParameterizedProblem(); StokesSolver *solver = NULL; - if (solver_type == "stokes") { solver = new StokesSolver; } - else if (solver_type == "steady-ns") { solver = new SteadyNSSolver; } + if (solver_type == "stokes") + { + solver = new StokesSolver; + } + else if (solver_type == "steady-ns") + { + solver = new SteadyNSSolver; + } - if (!solver->UseRom()) mfem_error("ROM must be enabled for supremizer enrichment!\n"); + if (!solver->UseRom()) + mfem_error("ROM must be enabled for supremizer enrichment!\n"); solver->InitVariables(); // This time needs to be ROMHandler, in order not to run StokesSolver::LoadSupremizer. @@ -325,7 +350,8 @@ void BuildROM(MPI_Comm comm) MultiBlockSolver *test = NULL; test = InitSolver(); - if (!test->UseRom()) mfem_error("ROM must be enabled for BuildROM!\n"); + if (!test->UseRom()) + mfem_error("ROM must be enabled for BuildROM!\n"); test->InitVariables(); // test->InitVisualization(); @@ -336,9 +362,9 @@ void BuildROM(MPI_Comm comm) // TODO: there are skippable operations depending on rom/fom mode. test->BuildOperators(); test->SetupBCOperators(); - + test->LoadReducedBasis(); - + TopologyHandlerMode topol_mode = test->GetTopologyMode(); ROMBuildingLevel save_operator = test->GetROMHandler()->GetBuildingLevel(); @@ -348,29 +374,29 @@ void BuildROM(MPI_Comm comm) switch (save_operator) { - case ROMBuildingLevel::COMPONENT: - { - if (topol_mode == TopologyHandlerMode::SUBMESH) - mfem_error("Submesh does not support component rom building level!\n"); + case ROMBuildingLevel::COMPONENT: + { + if (topol_mode == TopologyHandlerMode::SUBMESH) + mfem_error("Submesh does not support component rom building level!\n"); - test->AllocateROMElements(); - test->BuildROMElements(); - std::string filename = test->GetROMHandler()->GetOperatorPrefix() + ".h5"; - test->SaveROMElements(filename); - break; - } - case ROMBuildingLevel::GLOBAL: - { - test->ProjectOperatorOnReducedBasis(); - test->SaveROMOperator(); - break; - } - case ROMBuildingLevel::NONE: - { - printf("BuildROM - ROM building level is set to none. No ROM is saved.\n"); - break; - } - } // switch (save_operator) + test->AllocateROMElements(); + test->BuildROMElements(); + std::string filename = test->GetROMHandler()->GetOperatorPrefix() + ".h5"; + test->SaveROMElements(filename); + break; + } + case ROMBuildingLevel::GLOBAL: + { + test->ProjectOperatorOnReducedBasis(); + test->SaveROMOperator(); + break; + } + case ROMBuildingLevel::NONE: + { + printf("BuildROM - ROM building level is set to none. No ROM is saved.\n"); + break; + } + } // switch (save_operator) test->SaveBasisVisualization(); @@ -384,7 +410,7 @@ double SingleRun(MPI_Comm comm, const std::string output_file) { RandomSampleGenerator *generator = new RandomSampleGenerator(comm); generator->SetParamSpaceSizes(); - int idx = UniformRandom(0, generator->GetTotalSampleSize()-1); + int idx = UniformRandom(0, generator->GetTotalSampleSize() - 1); // NOTE: this will change config.dict_ generator->SetSampleParams(idx); delete generator; @@ -408,8 +434,10 @@ double SingleRun(MPI_Comm comm, const std::string output_file) const int num_var = test->GetNumVar(); Vector rom_assemble(1), rom_solve(1), fom_assemble(1), fom_solve(1), error(num_var); - rom_assemble = -1.0; rom_solve = -1.0; - fom_assemble = -1.0; fom_solve = -1.0; + rom_assemble = -1.0; + rom_solve = -1.0; + fom_assemble = -1.0; + fom_solve = -1.0; error = -1.0; ROMHandlerBase *rom = NULL; @@ -435,46 +463,46 @@ double SingleRun(MPI_Comm comm, const std::string output_file) switch (save_operator) { - case ROMBuildingLevel::COMPONENT: - { - if (topol_mode == TopologyHandlerMode::SUBMESH) - mfem_error("Submesh does not support component rom building level!\n"); - - printf("Loading component operator file.. "); - test->AllocateROMElements(); - std::string filename = rom->GetOperatorPrefix() + ".h5"; - test->LoadROMElements(filename); - test->AssembleROM(); - break; - } - case ROMBuildingLevel::GLOBAL: - { - printf("Loading global operator file.. "); - test->LoadROMOperatorFromFile(); - break; - } - case ROMBuildingLevel::NONE: - { - printf("Building operator file all the way from FOM.. "); - test->BuildDomainOperators(); - test->SetupDomainBCOperators(); - test->AssembleOperator(); - test->ProjectOperatorOnReducedBasis(); - break; - } + case ROMBuildingLevel::COMPONENT: + { + if (topol_mode == TopologyHandlerMode::SUBMESH) + mfem_error("Submesh does not support component rom building level!\n"); + + printf("Loading component operator file.. "); + test->AllocateROMElements(); + std::string filename = rom->GetOperatorPrefix() + ".h5"; + test->LoadROMElements(filename); + test->AssembleROM(); + break; + } + case ROMBuildingLevel::GLOBAL: + { + printf("Loading global operator file.. "); + test->LoadROMOperatorFromFile(); + break; + } + case ROMBuildingLevel::NONE: + { + printf("Building operator file all the way from FOM.. "); + test->BuildDomainOperators(); + test->SetupDomainBCOperators(); + test->AssembleOperator(); + test->ProjectOperatorOnReducedBasis(); + break; + } } printf("Done!\n"); printf("Projecting RHS to ROM.. "); test->ProjectRHSOnReducedBasis(); printf("Done!\n"); - } // if (test->UseRom()) + } // if (test->UseRom()) else { test->BuildDomainOperators(); test->SetupDomainBCOperators(); test->AssembleOperator(); - } // not if (test->UseRom()) + } // not if (test->UseRom()) solveTimer.Stop(); printf("%s-assemble time: %f seconds.\n", solveType.c_str(), solveTimer.RealTime()); @@ -573,7 +601,7 @@ double SingleRun(MPI_Comm comm, const std::string output_file) // Save solution and visualization. test->SaveSolution(); test->SaveVisualization(); - + delete test; delete problem; @@ -581,3 +609,62 @@ double SingleRun(MPI_Comm comm, const std::string output_file) return error.Max(); } +void TEMPRunAndCompare() +{ + LinElastSolver test; + test.InitVariables(); + const std::string visual_path = "test_"+ test.GetVisualizationPrefix(); + std::string sol_file = "test_"+ test.GetSolutionFilePrefix(); + sol_file += ".h5"; + + cout<<"sol_file is: "< mats; mats.SetSize(topol_data.numSub, topol_data.numSub);