diff --git a/.gitignore b/.gitignore index 0b5ccc801a6f7627529bc44c66358512f3bc8845..c863b8d3279b49d51fd0f92421ea97a139f76acd 100644 --- a/.gitignore +++ b/.gitignore @@ -12,3 +12,5 @@ Source/Solver/solver_fcbiilu2/fcbiilu2.cpp Source/Solver/solver_k3biilu2/k3d.cpp Source/Solver/solver_k3biilu2/k3d.h + +.vscode/* diff --git a/CMakeLists.txt b/CMakeLists.txt index ce755deb6a17304a814babe412009fef75705320..311f5592a4458d6657685d879f95ca0d2e1f96a6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -265,6 +265,9 @@ endif(COMPILE_TESTS) set(INMOST_INSTALL_HEADERS Source/Headers/inmost.h Source/Headers/inmost_autodiff.h + Source/Headers/inmost_residual.h + Source/Headers/inmost_model.h + Source/Headers/inmost_operator.h Source/Headers/inmost_common.h Source/Headers/inmost_data.h Source/Headers/inmost_dense.h @@ -276,6 +279,7 @@ set(INMOST_INSTALL_HEADERS Source/Headers/inmost.h Source/Headers/inmost_sparse.h Source/Headers/inmost_xml.h Source/Headers/inmost_variable.h + Source/Headers/inmost_block_variable.h Source/Headers/container.hpp) @@ -295,6 +299,9 @@ set_property(TARGET inmost PROPERTY PUBLIC_HEADER "${PROJECT_BINARY_DIR}/inmost_options.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_autodiff.h" + "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_residual.h" + "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_operator.h" + "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_model.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_common.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_data.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_dense.h" @@ -305,6 +312,7 @@ set_property(TARGET inmost PROPERTY PUBLIC_HEADER "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_solver.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_sparse.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_variable.h" + "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_block_variable.h" "${PROJECT_SOURCE_DIR}/Source/Headers/inmost_xml.h" "${PROJECT_SOURCE_DIR}/Source/Headers/container.hpp") diff --git a/Examples/DrawMatrix/main.cpp b/Examples/DrawMatrix/main.cpp index 6a630ea85ad76121a2ae93bbc3fb045b19668740..1a2066471794a15c67aafe7ed1322678a2875858 100644 --- a/Examples/DrawMatrix/main.cpp +++ b/Examples/DrawMatrix/main.cpp @@ -420,7 +420,7 @@ void draw() glBegin(GL_QUADS); for(Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it) for(Sparse::Row::iterator jt = it->Begin(); jt != it->End(); ++jt) - if( jt->first != it - m->Begin() ) + if( jt->first != it - m->Begin() && jt->second) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(jt->first));//, sqrt((jt->second-min)/(max-min))); DrawEntry(jt->first, m->Size() - (it - m->Begin()));//, sqrt((jt->second-min)/(max-min))); glEnd(); @@ -451,17 +451,17 @@ void draw() } glEnd(); - zoom += 1; + //zoom += 1; glColor3f(1.0, 0, 0); glBegin(GL_QUADS); for (Sparse::Matrix::iterator it = m->Begin(); it != m->End(); ++it) { int ind = it - m->Begin(); - if (fabs((*it)[ind]) < 1e-9) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind)); + if (fabs((*it)[ind]) < 1e-13) //DrawEntry(ord->position((it - m->Begin())), m->Size() - ord->position(ind)); DrawEntry((it - m->Begin()), m->Size() - ind); } glEnd(); - zoom -= 1; + //zoom -= 1; if (CommonInput != NULL) { diff --git a/Examples/GridTools/CMakeLists.txt b/Examples/GridTools/CMakeLists.txt index a53dc37584d1cebdcc12adf5b0000e51dab765de..960f5d4c3141017da6f301a0802bef45e9101f4c 100644 --- a/Examples/GridTools/CMakeLists.txt +++ b/Examples/GridTools/CMakeLists.txt @@ -9,6 +9,10 @@ add_executable(SliceFunc slice_func.cpp) add_executable(SplitNonplanar split_nonplanar.cpp) add_executable(CollapseDegenerate collapse_degenerate.cpp) add_executable(Bnd2Stl bnd2stl.cpp) +add_executable(Sector sector.cpp) +add_executable(Sew sew.cpp) +add_executable(Scale scale.cpp) +add_executable(Move move.cpp) target_link_libraries(FixFaults inmost) if(USE_MPI) @@ -106,4 +110,44 @@ if(USE_MPI) endif(USE_MPI) install(TARGETS Bnd2Stl EXPORT inmost-targets RUNTIME DESTINATION bin) +target_link_libraries(Sector inmost) +if(USE_MPI) + message("linking Sector with MPI") + target_link_libraries(Sector ${MPI_LIBRARIES}) + if(MPI_LINK_FLAGS) + set_target_properties(Sector PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") + endif() +endif(USE_MPI) +install(TARGETS Sector EXPORT inmost-targets RUNTIME DESTINATION bin) + +target_link_libraries(Sew inmost) +if(USE_MPI) + message("linking Sew with MPI") + target_link_libraries(Sew ${MPI_LIBRARIES}) + if(MPI_LINK_FLAGS) + set_target_properties(Sew PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") + endif() +endif(USE_MPI) +install(TARGETS Sew EXPORT inmost-targets RUNTIME DESTINATION bin) + +target_link_libraries(Scale inmost) +if(USE_MPI) + message("linking Scale with MPI") + target_link_libraries(Scale ${MPI_LIBRARIES}) + if(MPI_LINK_FLAGS) + set_target_properties(Scale PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") + endif() +endif(USE_MPI) +install(TARGETS Scale EXPORT inmost-targets RUNTIME DESTINATION bin) + + +target_link_libraries(Move inmost) +if(USE_MPI) + message("linking Move with MPI") + target_link_libraries(Move ${MPI_LIBRARIES}) + if(MPI_LINK_FLAGS) + set_target_properties(Move PROPERTIES LINK_FLAGS "${MPI_LINK_FLAGS}") + endif() +endif(USE_MPI) +install(TARGETS Move EXPORT inmost-targets RUNTIME DESTINATION bin) diff --git a/Examples/GridTools/move.cpp b/Examples/GridTools/move.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1baea0ace51176942415f72acf11e837a9766e2d --- /dev/null +++ b/Examples/GridTools/move.cpp @@ -0,0 +1,73 @@ +#include +#include "inmost.h" + + +using namespace INMOST; +typedef Storage::real real; +typedef Storage::real_array real_array; + +int main(int argc, char ** argv) +{ + if( argc < 2 ) + { + printf("Usage: %s input_mesh [move_x=1] [move_y=1] [move_z=1] [output_mesh]\n",argv[0]); + return -1; + } + + Mesh m; + m.Load(argv[1]); + double mx = argc > 2 ? atof(argv[2]) : 1; + double my = argc > 3 ? atof(argv[3]) : 1; + double mz = argc > 4 ? atof(argv[4]) : 1; + std::string save_file = argc > 5 ? std::string(argv[5]) : "out.vtk"; + + double xmax = -1.0e20, xmin = 1.0e20; + double ymax = -1.0e20, ymin = 1.0e20; + double zmax = -1.0e20, zmin = 1.0e20; + + for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n) + { + double x = n->Coords()[0], y = n->Coords()[1], z = n->Coords()[2]; + if( x > xmax ) xmax = x; + if( x < xmin ) xmin = x; + if( y > ymax ) ymax = y; + if( y < ymin ) ymin = y; + if( z > zmax ) zmax = z; + if( z < zmin ) zmin = z; + } + + std::cout << "x: " << xmin << ":" << xmax << std::endl; + std::cout << "y: " << ymin << ":" << ymax << std::endl; + std::cout << "z: " << zmin << ":" << zmax << std::endl; + + for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n) + { + n->Coords()[0] += mx; + n->Coords()[1] += my; + n->Coords()[2] += mz; + } + + xmax = -1.0e20, xmin = 1.0e20; + ymax = -1.0e20, ymin = 1.0e20; + zmax = -1.0e20, zmin = 1.0e20; + + for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n) + { + double x = n->Coords()[0], y = n->Coords()[1], z = n->Coords()[2]; + if( x > xmax ) xmax = x; + if( x < xmin ) xmin = x; + if( y > ymax ) ymax = y; + if( y < ymin ) ymin = y; + if( z > zmax ) zmax = z; + if( z < zmin ) zmin = z; + } + + std::cout << "x: " << xmin << ":" << xmax << std::endl; + std::cout << "y: " << ymin << ":" << ymax << std::endl; + std::cout << "z: " << zmin << ":" << zmax << std::endl; + + std::cout << "Save to " << save_file << std::endl; + m.Save(save_file); + + return 0; +} diff --git a/Examples/GridTools/scale.cpp b/Examples/GridTools/scale.cpp new file mode 100644 index 0000000000000000000000000000000000000000..40c46e421b388c320061b1370c55b362040d19e2 --- /dev/null +++ b/Examples/GridTools/scale.cpp @@ -0,0 +1,73 @@ +#include +#include "inmost.h" + + +using namespace INMOST; +typedef Storage::real real; +typedef Storage::real_array real_array; + +int main(int argc, char ** argv) +{ + if( argc < 2 ) + { + printf("Usage: %s input_mesh [scale_x=1] [scale_y=1] [scale_z=1] [output_mesh]\n",argv[0]); + return -1; + } + + Mesh m; + m.Load(argv[1]); + double sx = argc > 2 ? atof(argv[2]) : 1; + double sy = argc > 3 ? atof(argv[3]) : 1; + double sz = argc > 4 ? atof(argv[4]) : 1; + std::string save_file = argc > 5 ? std::string(argv[5]) : "out.vtk"; + + double xmax = -1.0e20, xmin = 1.0e20; + double ymax = -1.0e20, ymin = 1.0e20; + double zmax = -1.0e20, zmin = 1.0e20; + + for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n) + { + double x = n->Coords()[0], y = n->Coords()[1], z = n->Coords()[2]; + if( x > xmax ) xmax = x; + if( x < xmin ) xmin = x; + if( y > ymax ) ymax = y; + if( y < ymin ) ymin = y; + if( z > zmax ) zmax = z; + if( z < zmin ) zmin = z; + } + + std::cout << "x: " << xmin << ":" << xmax << std::endl; + std::cout << "y: " << ymin << ":" << ymax << std::endl; + std::cout << "z: " << zmin << ":" << zmax << std::endl; + + for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n) + { + n->Coords()[0] *= sx; + n->Coords()[1] *= sy; + n->Coords()[2] *= sz; + } + + xmax = -1.0e20, xmin = 1.0e20; + ymax = -1.0e20, ymin = 1.0e20; + zmax = -1.0e20, zmin = 1.0e20; + + for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n) + { + double x = n->Coords()[0], y = n->Coords()[1], z = n->Coords()[2]; + if( x > xmax ) xmax = x; + if( x < xmin ) xmin = x; + if( y > ymax ) ymax = y; + if( y < ymin ) ymin = y; + if( z > zmax ) zmax = z; + if( z < zmin ) zmin = z; + } + + std::cout << "x: " << xmin << ":" << xmax << std::endl; + std::cout << "y: " << ymin << ":" << ymax << std::endl; + std::cout << "z: " << zmin << ":" << zmax << std::endl; + + std::cout << "Save to " << save_file << std::endl; + m.Save(save_file); + + return 0; +} diff --git a/Examples/GridTools/sector.cpp b/Examples/GridTools/sector.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b77b581c5ec4be0b4834a6ec74b617fad6bcfbf7 --- /dev/null +++ b/Examples/GridTools/sector.cpp @@ -0,0 +1,81 @@ +#include +#include "inmost.h" + + +using namespace INMOST; +typedef Storage::real real; +typedef Storage::real_array real_array; + +int main(int argc, char ** argv) +{ + if( argc < 2 ) + { + printf("Usage: %s input_mesh [rotation_angle=0 degrees] [output_mesh]\n",argv[0]); + return -1; + } + + Mesh m; + m.Load(argv[1]); + + const double pi = 3.1415926536; + double theta = 0, ct, st; + if( argc > 2 ) theta = atof(argv[2])/180.0*pi; + ct = cos(theta); + st = sin(theta); + + double xmin = 1.0e20, xmax = -1.0e20; + double ymin = 1.0e20, ymax = -1.0e20; + + for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n) + { + double x = n->Coords()[0]; + double y = n->Coords()[1]; + if( x > xmax ) xmax = x; + if( x < xmin ) xmin = x; + if( y > ymax ) ymax = y; + if( y < ymin ) ymin = y; + } + + std::cout << "x: " << xmin << ":" << xmax << std::endl; + std::cout << "y: " << ymin << ":" << ymax << std::endl; + + for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n) + { + double x = n->Coords()[0], mx = x; + double y = n->Coords()[1], my = y; + double alpha = pi/4*(2*y-1); + mx += 1 - x*(1-cos(alpha)) + x*cos(2*alpha)*0.2; + my += x*sin(alpha); + n->Coords()[0] = (mx-0.5)*ct - (my-0.5)*st + 0.5; + n->Coords()[1] = (mx-0.5)*st + (my-0.5)*ct + 0.5; + } + + xmin = 1.0e20, xmax = -1.0e20; + ymin = 1.0e20, ymax = -1.0e20; + + for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n) + { + double x = n->Coords()[0]; + double y = n->Coords()[1]; + if( x > xmax ) xmax = x; + if( x < xmin ) xmin = x; + if( y > ymax ) ymax = y; + if( y < ymin ) ymin = y; + } + + std::cout << "x: " << xmin << ":" << xmax << std::endl; + std::cout << "y: " << ymin << ":" << ymax << std::endl; + + if( argc > 3 ) + { + std::cout << "Save to " << argv[3] << std::endl; + m.Save(argv[3]); + } + else + { + std::cout << "Save to out.vtk" << std::endl; + m.Save("out.vtk"); + } + + return 0; +} diff --git a/Examples/GridTools/sew.cpp b/Examples/GridTools/sew.cpp new file mode 100644 index 0000000000000000000000000000000000000000..def4811a9b81d7277c277e98620b4a7fc19ac428 --- /dev/null +++ b/Examples/GridTools/sew.cpp @@ -0,0 +1,44 @@ +#include +#include "inmost.h" + + +using namespace INMOST; +typedef Storage::real real; +typedef Storage::real_array real_array; + +int main(int argc, char ** argv) +{ + if( argc < 2 ) + { + printf("Usage: %s input_mesh1 [input_mesh2] ... [output_mesh]\n",argv[0]); + return -1; + } + + Mesh m; + for(int k = 1 ; k < argc; ++k) + m.Load(argv[k]); + + double xmax = -1.0e20, xmin = 1.0e20; + double ymax = -1.0e20, ymin = 1.0e20; + double zmax = -1.0e20, zmin = 1.0e20; + + for(Mesh::iteratorNode n = m.BeginNode(); n != m.EndNode(); ++n) + { + double x = n->Coords()[0], y = n->Coords()[1], z = n->Coords()[2]; + if( x > xmax ) xmax = x; + if( x < xmin ) xmin = x; + if( y > ymax ) ymax = y; + if( y < ymin ) ymin = y; + if( z > zmax ) zmax = z; + if( z < zmin ) zmin = z; + } + + std::cout << "x: " << xmin << ":" << xmax << std::endl; + std::cout << "y: " << ymin << ":" << ymax << std::endl; + std::cout << "z: " << zmin << ":" << zmax << std::endl; + + std::cout << "Save to out.pmf" << std::endl; + m.Save("out.pmf"); + + return 0; +} diff --git a/Examples/GridTools/slice_func.cpp b/Examples/GridTools/slice_func.cpp index 0dcbc3636b7dc62b052b4a64d97d8fd83da23819..68dd4d83e13436b583f874d9311a0009296bcad3 100644 --- a/Examples/GridTools/slice_func.cpp +++ b/Examples/GridTools/slice_func.cpp @@ -34,6 +34,10 @@ double func(double x, double y, double z, int n) if (y < Ly) return y - Ly; return fmin( fmin(x-Lx,Rx-x), fmin(y-Ly, Ry-y) ); } + else if( n == 6 ) + { + return std::min(std::min(sqrt((x-0.48)*(x-0.48) + (z-1)*(z-1))-0.25,sqrt((y-0.53)*(y-0.53) + (z-3)*(z-3))-0.24),sqrt((x-0.6)*(x-0.6) + (z-5)*(z-5))-0.3); + } return 1; } @@ -1057,7 +1061,7 @@ int main(int argc, char ** argv) nvv++; } vv /= nvv; - if( v < vv*0.15 ) + if( v < vv*0.01 ) { if( false ) { @@ -1117,7 +1121,7 @@ int main(int argc, char ** argv) } for(Mesh::iteratorCell it = m.BeginCell(); it != m.EndCell(); ++it) - if( material[*it] == 0 ) + if( material[*it] == 0 || it->Integer(collapse) ) it->Delete(); for(ElementType etype = FACE; etype >= NODE; etype = PrevElementType(etype) ) @@ -1130,4 +1134,4 @@ int main(int argc, char ** argv) m.ReleaseMarker(mrk,EDGE|FACE|NODE); m.Save(grid_out); return 0; -} \ No newline at end of file +} diff --git a/Examples/OldDrawGrid/CMakeLists.txt b/Examples/OldDrawGrid/CMakeLists.txt index 9ca8254e451a4420bde7a42a8f00139ed91bd839..69cea5173d04fe071c66795028b193bcae21dbd4 100644 --- a/Examples/OldDrawGrid/CMakeLists.txt +++ b/Examples/OldDrawGrid/CMakeLists.txt @@ -29,7 +29,9 @@ set(SOURCE main.cpp picker.h picker.cpp clipper.h - clipper.cpp) + clipper.cpp + vector.h + vector.cpp) find_package(OpenGL) find_package(GLUT) diff --git a/Examples/OldDrawGrid/clipboard.cpp b/Examples/OldDrawGrid/clipboard.cpp index b44d6a1bdbf686f5cd406cbb7a18b44fb0833d78..6c9da4ee1723d33cfd5dac1d1f9c2c04b5a3cf0f 100644 --- a/Examples/OldDrawGrid/clipboard.cpp +++ b/Examples/OldDrawGrid/clipboard.cpp @@ -73,14 +73,14 @@ bool setTextToPasteboard(std::string str) PasteboardCreate( kPasteboardClipboard, &pasteboard ); err = PasteboardClear( pasteboard ); - require_noerr( err, PasteboardClear_FAILED ); + //require_noerr( err, PasteboardClear_FAILED ); CFDataRef data; data = CFDataCreate(kCFAllocatorDefault, (UInt8*)byteArrayIndex, strlen(byteArrayIndex)+1); err = PasteboardPutItemFlavor( pasteboard, (PasteboardItemID)1, kUTTypeUTF8PlainText, data, 0); - require_noerr( err, PasteboardPutItemFlavor_FAILED ); + //require_noerr( err, PasteboardPutItemFlavor_FAILED ); PasteboardPutItemFlavor_FAILED: PasteboardClear_FAILED: @@ -101,14 +101,14 @@ std::string getTextFromPasteboard() err = badPasteboardSyncErr; err = PasteboardGetItemCount( inPasteboard, &itemCount ); - require_noerr( err, CantGetPasteboardItemCount ); + //require_noerr( err, CantGetPasteboardItemCount ); for( int itemIndex = 1; itemIndex <= itemCount; itemIndex++ ) { PasteboardItemID itemID; CFDataRef flavorData; err = PasteboardGetItemIdentifier( inPasteboard, itemIndex, &itemID ); - require_noerr( err, CantGetPasteboardItemIdentifier ); + //require_noerr( err, CantGetPasteboardItemIdentifier ); err = PasteboardCopyItemFlavorData( inPasteboard, itemID, CFSTR("public.utf8-plain-text"), &flavorData ); if(err==noErr)data = (char*)CFDataGetBytePtr(flavorData); if( data!=NULL && err==noErr ) diff --git a/Examples/OldDrawGrid/main.cpp b/Examples/OldDrawGrid/main.cpp index 65cbb7c5e2adb3b1bd65e18c6fb3b02dd41a102a..89a244485a61d71c69a38216c95a261d57e15d46 100644 --- a/Examples/OldDrawGrid/main.cpp +++ b/Examples/OldDrawGrid/main.cpp @@ -23,6 +23,7 @@ #include "tga.h" #include "screenshot.h" #include "volumetric.h" +#include "vector.h" #include "input.h" #include "picker.h" #include "clipper.h" @@ -206,6 +207,7 @@ std::vector clip_boundary; volumetric * CommonVolumetricView; +Vectors * CommonVectors = NULL; @@ -1294,6 +1296,9 @@ void draw_screen() if( bndupdate ) CommonVolumetricView->camera(campos,interactive); CommonVolumetricView->draw(interactive); } + + if( CommonVectors ) + CommonVectors->Draw(interactive); for(int k = 0; k < streamlines.size(); ++k) @@ -1504,6 +1509,11 @@ void draw_screen() correct_input = true; streamlines.clear(); glutPostRedisplay(); + if( CommonVectors ) + { + delete CommonVectors; + CommonVectors = NULL; + } } } @@ -1515,79 +1525,90 @@ void draw_screen() visualization_type = NONE; for (size_t q = 0; q < stype.size(); ++q) stype[q] = tolower(stype[q]); - if (stype == "node") visualization_type = NODE; - else if (stype == "edge") visualization_type = EDGE; - else if (stype == "face") visualization_type = FACE; - else if (stype == "cell") visualization_type = CELL; - else if (stype == "eset") visualization_type = ESET; - if (visualization_type == NONE) - printf("unknown element type %s\n", typen); + if (stype == "scale") + { + double scale = atof(visualization_prompt + k + 1); + std::cout << "input scale: " << scale << std::endl; + if( CommonVectors ) + CommonVectors->SetScale(scale); + correct_input = true; + glutPostRedisplay(); + } else { - disp_e = InvalidElement(); - bool is_number = true; - for (l = k + 1; l < slen; ++l) - if (!isdigit(visualization_prompt[l])) - is_number = false; - if (is_number) + if (stype == "node") visualization_type = NODE; + else if (stype == "edge") visualization_type = EDGE; + else if (stype == "face") visualization_type = FACE; + else if (stype == "cell") visualization_type = CELL; + else if (stype == "eset") visualization_type = ESET; + if (visualization_type == NONE) + printf("unknown element type %s\n", typen); + else { - comp = atoi(visualization_prompt + k + 1); - visualization_prompt[k] = ':'; + disp_e = InvalidElement(); + bool is_number = true; + for (l = k + 1; l < slen; ++l) + if (!isdigit(visualization_prompt[l])) + is_number = false; + if (is_number) + { + comp = atoi(visualization_prompt + k + 1); + visualization_prompt[k] = ':'; - correct_input = true; + correct_input = true; - if (mesh->isValidElement(visualization_type, comp)) - { - printf("Display data for %s:%d\n", typen, comp); - disp_e = mesh->ElementByLocalID(visualization_type, comp); + if (mesh->isValidElement(visualization_type, comp)) + { + printf("Display data for %s:%d\n", typen, comp); + disp_e = mesh->ElementByLocalID(visualization_type, comp); + } + else + printf("No valid element at %s:%d\n", typen, comp); } - else - printf("No valid element at %s:%d\n", typen, comp); - } - else if (visualization_type == ESET) - { - std::string name = std::string(visualization_prompt + k + 1); - visualization_prompt[k] = ':'; - correct_input = true; - disp_e = mesh->GetSet(name); - if (disp_e.isValid()) - printf("Display data for %s:%d\n", typen, disp_e.LocalID()); - else - printf("Cannot find set with name %s\n", name.c_str()); - } - - if (disp_e.isValid()) - { - if (disp_e.GetElementType() != ESET) + else if (visualization_type == ESET) { - disp_e->Centroid(shift); - for (int r = 0; r < 3; ++r) - shift[r] = -shift[r]; + std::string name = std::string(visualization_prompt + k + 1); + visualization_prompt[k] = ':'; + correct_input = true; + disp_e = mesh->GetSet(name); + if (disp_e.isValid()) + printf("Display data for %s:%d\n", typen, disp_e.LocalID()); + else + printf("Cannot find set with name %s\n", name.c_str()); } - else + + if (disp_e.isValid()) { - shift[0] = shift[1] = shift[2] = 0; - ElementSet s = disp_e.getAsSet(); - int nelem = 0; - for (ElementSet::iterator it = s.Begin(); it != s.End(); ++it) + if (disp_e.GetElementType() != ESET) { - double cnt[3]; - it->Centroid(cnt); - shift[0] += cnt[0]; - shift[1] += cnt[1]; - shift[2] += cnt[2]; - nelem++; + disp_e->Centroid(shift); + for (int r = 0; r < 3; ++r) + shift[r] = -shift[r]; + } + else + { + shift[0] = shift[1] = shift[2] = 0; + ElementSet s = disp_e.getAsSet(); + int nelem = 0; + for (ElementSet::iterator it = s.Begin(); it != s.End(); ++it) + { + double cnt[3]; + it->Centroid(cnt); + shift[0] += cnt[0]; + shift[1] += cnt[1]; + shift[2] += cnt[2]; + nelem++; + } + shift[0] /= (double)nelem; + shift[1] /= (double)nelem; + shift[2] /= (double)nelem; + for (int r = 0; r < 3; ++r) + shift[r] = -shift[r]; } - shift[0] /= (double)nelem; - shift[1] /= (double)nelem; - shift[2] /= (double)nelem; - for (int r = 0; r < 3; ++r) - shift[r] = -shift[r]; } } } - } if( k < slen && l < slen && l+1 < slen ) @@ -1606,16 +1627,19 @@ void draw_screen() comp = ENUMUNDEF; else if (std::string(visualization_prompt + l + 1) == "streamline") comp = ENUMUNDEF-2; + else if (std::string(visualization_prompt + l + 1) == "vector" || + std::string(visualization_prompt + l + 1) == "vec") + comp = ENUMUNDEF-3; else { - std::cout << "unknown name for component, expected number or 'mag'" << std::endl; + std::cout << "unknown name for component, expected number or 'mag' or 'vec' or 'streamline'" << std::endl; comp = ENUMUNDEF-1; } visualization_prompt[k] = ':'; visualization_prompt[l] = ':'; printf("type %s name %s comp %d\n",typen,name,comp); std::string stype(typen), sname(name); - if (mesh->HaveTag(sname) && comp == ENUMUNDEF - 2) + if (mesh->HaveTag(sname) && (comp == ENUMUNDEF - 2 || comp == ENUMUNDEF - 3)) { ElementType vel_def = NONE; for (size_t q = 0; q < stype.size(); ++q) @@ -1628,8 +1652,16 @@ void draw_screen() else if (stype == "face") vel_def = FACE; else if (stype == "cell") vel_def = CELL; - streamlines.clear(); - BuildStreamlines(mesh,mesh->GetTag(sname),vel_def,streamlines); + if( comp == ENUMUNDEF - 2 ) + { + streamlines.clear(); + BuildStreamlines(mesh,mesh->GetTag(sname),vel_def,streamlines); + } + else if( comp == ENUMUNDEF - 3 ) + { + if( CommonVectors ) delete CommonVectors; + CommonVectors = new Vectors(mesh,mesh->GetTag(sname),vel_def); + } } else if( mesh->HaveTag(sname) && comp != ENUMUNDEF-1) { @@ -2017,6 +2049,9 @@ void svg_draw(std::ostream & file) //if( bndupdate ) CommonVolumetricView->camera(campos,interactive); //CommonVolumetricView->draw(interactive); } + + if( CommonVectors ) + CommonVectors->SVGDraw(file,modelview,projection,viewport); for(int k = 0; k < streamlines.size(); ++k) diff --git a/Examples/OldDrawGrid/vector.cpp b/Examples/OldDrawGrid/vector.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0efd413a741cdff805723ddf90feee8266bd869f --- /dev/null +++ b/Examples/OldDrawGrid/vector.cpp @@ -0,0 +1,107 @@ +#include "vector.h" +#include "inc_glut.h" +#include "color_bar.h" +#include "svg_line.h" + +namespace INMOST +{ + + + + static GLUquadric * cylqs = NULL; + static void drawcylinder(coord a, coord b, double width) + { + double matrix[16]; + if (cylqs == NULL) + { + cylqs = gluNewQuadric(); + gluQuadricNormals(cylqs, GLU_SMOOTH); + gluQuadricOrientation(cylqs, GLU_OUTSIDE); + gluQuadricDrawStyle(cylqs, GLU_FILL);//GLU_SILHOUETTE + } + glPushMatrix(); + glTranslated(a[0], a[1], a[2]); + get_matrix(a, b, matrix); + glMultMatrixd(matrix); + gluCylinder(cylqs, width, width, sqrt((b - a) ^ (b - a)), 4, 2); + glPopMatrix(); + } + + Vectors::Vectors(Mesh * m, TagRealArray t, ElementType etype) : m(m), etype(etype) + { + double cmin[3] = { 1.0e20, 1.0e20, 1.0e20}; + double cmax[3] = {-1.0e20,-1.0e20,-1.0e20}; + for(int k = 0; k < m->NodeLastLocalID(); ++k) if( m->isValidNode(k) ) + { + Node n = m->NodeByLocalID(k); + for(int l = 0; l < 3; ++l) + { + if( cmin[l] > n->Coords()[l] ) cmin[l] = n->Coords()[l]; + if( cmax[l] < n->Coords()[l] ) cmax[l] = n->Coords()[l]; + } + } + double diag = 0; + for(int l = 0; l < 3; ++l) + diag += (cmax[l]-cmin[l])*(cmax[l]-cmin[l]); + diag = sqrt(diag); //estimate mesh size + max_length = 0; + for(int k = 0; k < m->LastLocalID(etype); ++k) if( m->isValidElement(etype,k) ) + { + Element e = m->ElementByLocalID(etype,k); + if( e.HaveData(t) ) + { + vec_t add; + e->Centroid(add.cnt.data()); + add.eid = k; + for(int l = 0; l < std::min(t[e].size(),3u); ++l) + add.dir[l] = t[e][l]; + double l = add.dir.length(); + if( l > max_length ) max_length = l; + add.dir /= l; + add.length = l; + vecs.push_back(add); + } + } + scale = diag/100.0; + std::cout << "scale: " << scale << std::endl; + } + + void Vectors::Draw(int reduced) + { + glBegin(GL_LINES); + int pace = 1; + if( reduced ) + pace = std::max(1,std::min(15,(unsigned)vecs.size()/100)); + for (unsigned int i = 0; i < vecs.size(); i+=pace) + { + if (isColorBarEnabled()) + { + color_t c = GetColorBar()->pick_color(m->ElementByLocalID(etype,vecs[i].eid)->RealDF(GetVisualizationTag())); + c.set_color(); + } + glVertex3dv(vecs[i].cnt.data()); + glVertex3dv((vecs[i].cnt+vecs[i].dir*scale*vecs[i].length/max_length).data()); + } + glEnd(); + } + + + void Vectors::SVGDraw(std::ostream & file, double modelview[16], double projection[16], int viewport[4]) + { + for (unsigned int i = 0; i < vecs.size(); i++) + { + coord p2 = (vecs[i].cnt+vecs[i].dir*scale*vecs[i].length/max_length); + double * v0 = vecs[i].cnt.data(); + double * v1 = p2.data(); + + if (isColorBarEnabled()) + { + color_t c = GetColorBar()->pick_color(m->ElementByLocalID(etype,vecs[i].eid)->RealDF(GetVisualizationTag())); + file << "" << std::endl; + } + else file << "" << std::endl; + svg_line(file, v0[0], v0[1], v0[2], v1[0], v1[1], v1[2], modelview, projection, viewport); + file << "" << std::endl; + } + } +} diff --git a/Examples/OldDrawGrid/vector.h b/Examples/OldDrawGrid/vector.h new file mode 100644 index 0000000000000000000000000000000000000000..96a6752d869a206c55e1bfc0f54d40543a1bf2c2 --- /dev/null +++ b/Examples/OldDrawGrid/vector.h @@ -0,0 +1,46 @@ +#ifndef _VECTOR_H +#define _VECTOR_H + +#include "inmost.h" +#include "coord.h" + + +namespace INMOST +{ + + struct vec_t + { + coord cnt, dir; + double length; + int eid; + }; + + class Vectors + { + private: + double max_length; + double scale; + Mesh * m; + ElementType etype; + std::vector vecs; + public: + Vectors(Mesh * m, TagRealArray t, ElementType etype); + Vectors(const Vectors & b) + : max_length(b.max_length), scale(b.scale), m(b.m), etype(b.etype), vecs(b.vecs){} + Vectors & operator =(Vectors const & b) + { + max_length = b.max_length; + scale = b.scale; + m = b.m; + etype = b.etype; + vecs = b.vecs; + return *this; + } + void SetScale(double s) {scale = s;} + ~Vectors() { vecs.clear(); } + void Draw(int reduced); + void SVGDraw(std::ostream & file, double modelview[16], double projection[16], int viewport[4]); + }; + +} +#endif diff --git a/Source/Autodiff/CMakeLists.txt b/Source/Autodiff/CMakeLists.txt index ce56940c80a8690e197c636976040e28a000e103..bae0ee4c44f06bac5660eb6b7e094d6d40a5beae 100644 --- a/Source/Autodiff/CMakeLists.txt +++ b/Source/Autodiff/CMakeLists.txt @@ -1,6 +1,9 @@ set(SOURCE ${SOURCE} ${CMAKE_CURRENT_SOURCE_DIR}/autodiff.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/residual.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/model.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/operator.cpp PARENT_SCOPE ) diff --git a/Source/Autodiff/autodiff.cpp b/Source/Autodiff/autodiff.cpp index dee1ffa188b4ca89501f3aad071b0a8192662eaf..576bdabccfeec2daf6e7923e74579d0fe96e4058 100644 --- a/Source/Autodiff/autodiff.cpp +++ b/Source/Autodiff/autodiff.cpp @@ -105,7 +105,7 @@ namespace INMOST { Mesh * m = b.GetMeshLink(); - ElementType sparse = NONE;// b.GetElementType(); + ElementType sparse = NONE;// b.GetElementType(); for (ElementType q = NODE; q <= MESH; q = NextElementType(q)) if (q & b.GetElementType()) { for(unsigned unk = 0; unk < b.Size(); ++unk) @@ -141,7 +141,7 @@ namespace INMOST void Automatizator::UnregisterEntry(INMOST_DATA_ENUM_TYPE ind) { - assert(act_blocks[ind]); + assert(reg_blocks[ind]); if( reg_blocks[ind]->GetOffsetTag().isValid() ) reg_blocks[ind]->GetOffsetTag().GetMeshLink()->DeleteTag(reg_blocks[ind]->GetOffsetTag()); delete reg_blocks[ind]; @@ -149,6 +149,21 @@ namespace INMOST del_blocks.push_back(ind); act_blocks[ind] = false; } + + void Automatizator::DeactivateEntry(INMOST_DATA_ENUM_TYPE ind) + { + assert(reg_blocks[ind] != NULL); ///This block was not deleted + reg_blocks[ind]->reg_index = ENUMUNDEF; + act_blocks[ind] = false; + } + + void Automatizator::ActivateEntry(INMOST_DATA_ENUM_TYPE ind) + { + assert(reg_blocks[ind] != NULL); ///This block was not deleted + reg_blocks[ind]->reg_index = ind; + act_blocks[ind] = true; + } + void Automatizator::EnumerateEntries() { first_num = last_num = 0; @@ -159,7 +174,7 @@ namespace INMOST AbstractEntry & b = *reg_blocks[it]; TagInteger offset_tag = b.GetOffsetTag(); Mesh * m = offset_tag.GetMeshLink(); - for (ElementType etype = NODE; etype <= MESH; etype = etype << 1) + for (ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype)) if (offset_tag.isDefined(etype) && offset_tag.isSparse(etype)) { for(int kt = 0; kt < m->LastLocalID(etype); ++kt) if( m->isValidElement(etype,kt) ) @@ -397,7 +412,7 @@ namespace INMOST INMOST_DATA_ENUM_TYPE MultiEntry::GetValueComp(INMOST_DATA_ENUM_TYPE unk) const { unsigned pos = 0, k = 0; - while( pos + entries[k]->Size() < unk ) pos += entries[k++]->Size(); + while( pos + entries[k]->Size() <= unk ) pos += entries[k++]->Size(); assert(k < entries.size()); return entries[k]->GetValueComp(unk-pos); } @@ -405,7 +420,7 @@ namespace INMOST TagRealArray MultiEntry::GetValueTag(INMOST_DATA_ENUM_TYPE unk) const { unsigned pos = 0, k = 0; - while( pos + entries[k]->Size() < unk ) pos += entries[k++]->Size(); + while( pos + entries[k]->Size() <= unk ) pos += entries[k++]->Size(); assert(k < entries.size()); return entries[k]->GetValueTag(unk-pos); } @@ -464,114 +479,44 @@ namespace INMOST ret->AddTag(unknown_tags[k],unknown_comp[k]); return ret; } -#endif //USE_MESH - -#if defined(USE_SOLVER) - void Residual::GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const - { - start = residual.GetFirstIndex(); - end = residual.GetLastIndex(); - } - void Residual::SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end) - { - jacobian.SetInterval(beg,end); - residual.SetInterval(beg,end); - } - void Residual::ClearResidual() - { - for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) (*it) = 0.0; - } - void Residual::ClearJacobian() - { - for(Sparse::Matrix::iterator it = jacobian.Begin(); it != jacobian.End(); ++it) it->Clear(); - } - void Residual::ClearHessian() - { - for(Sparse::HessianMatrix::iterator it = hessian.Begin(); it != hessian.End(); ++it) it->Clear(); - } - void Residual::Clear() + + void AbstractEntry::SynchronizeData() { -#if defined(USE_OMP) -#pragma omp for -#endif //USE_OMP - for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) + //synchronize indices { - residual[k] = 0.0; - if( !jacobian.Empty() ) jacobian[k].Clear(); - if( !hessian.Empty() ) hessian[k].Clear(); + Mesh * m = NULL; + std::set exch_tags; + for(unsigned jt = 0; jt < Size(); ++jt) + { + if( m == NULL ) + m = GetValueTag(jt).GetMeshLink(); + else assert(m == GetValueTag(jt).GetMeshLink()); + exch_tags.insert(GetValueTag(jt)); + } + m->ExchangeData(std::vector(exch_tags.begin(),exch_tags.end()), GetElementType(),0); } } - INMOST_DATA_REAL_TYPE Residual::Norm() - { - INMOST_DATA_REAL_TYPE ret = 0; -#if defined(USE_OMP) -#pragma omp parallel for reduction(+:ret) -#endif //USE_OMP - for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) - ret += residual[k]*residual[k]; -#if defined(USE_MPI) - INMOST_DATA_REAL_TYPE tmp = ret; - MPI_Allreduce(&tmp, &ret, 1, INMOST_MPI_DATA_REAL_TYPE, MPI_SUM, jacobian.GetCommunicator()); -#endif - return sqrt(ret); - } - Residual & Residual::operator =(Residual const & other) - { - hessian = other.hessian; - jacobian = other.jacobian; - residual = other.residual; - return *this; - } - void Residual::Rescale(INMOST_DATA_ENUM_TYPE p) + + + void Automatizator::SynchronizeData() { -#if defined(USE_OMP) -#pragma omp parallel for -#endif //USE_OMP - for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) + //synchronize indices { - INMOST_DATA_REAL_TYPE norm = 0.0; - if( p == ENUMUNDEF ) //infinite norm - { - for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q) - if( norm < fabs(jacobian[k].GetValue(q)) ) - norm = fabs(jacobian[k].GetValue(q)); - } - else //p-norm + std::map > exch_tags; + ElementType exch_mask = NONE; + for (unsigned it = 0; it < reg_blocks.size(); ++it) if( act_blocks[it] ) { - for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q) - norm += pow(fabs(jacobian[k].GetValue(q)),p); - norm = pow(norm,1.0/p); + exch_mask |= reg_blocks[it]->GetElementType(); + for(unsigned jt = 0; jt < reg_blocks[it]->Size(); ++jt) + exch_tags[reg_blocks[it]->GetValueTag(jt).GetMeshLink()].insert(reg_blocks[it]->GetValueTag(jt)); } - if( norm ) + for(std::map >::iterator it = exch_tags.begin(); it != exch_tags.end(); ++it) { - norm = 1.0/norm; - residual[k] *= norm; - for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q) - jacobian[k].GetValue(q) *= norm; - if( !hessian.Empty() ) - for(INMOST_DATA_ENUM_TYPE q = 0; q < hessian[k].Size(); ++q) - hessian[k].GetValue(q) *= norm; + it->first->ExchangeData(std::vector(it->second.begin(),it->second.end()), exch_mask,0); } } } - Residual::Residual(std::string name, INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end, INMOST_MPI_Comm _comm) - : jacobian(name,start,end,_comm),residual(name,start,end,_comm), hessian(name,0,0,_comm) - { - } - Residual::Residual(const Residual & other) - : jacobian(other.jacobian), residual(other.residual), hessian(other.hessian) - { - } - - Matrix Residual::operator [](const AbstractMatrix & rows) - { - Matrix ret(rows.Rows(),rows.Cols()); - for(INMOST_DATA_ENUM_TYPE i = 0; i < rows.Rows(); ++i) - for(INMOST_DATA_ENUM_TYPE j = 0; j < rows.Cols(); ++j) - ret(i,j) = multivar_expression_reference(residual[rows(i,j)],&jacobian[rows(i,j)]); - return ret; - } -#endif //USE_SOLVER -}; +#endif //USE_MESH +} //namespace INMOST #endif //USE_AUTODIFF diff --git a/Source/Autodiff/model.cpp b/Source/Autodiff/model.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c370c389e350b2f907e8e67daff1d71a8ed83471 --- /dev/null +++ b/Source/Autodiff/model.cpp @@ -0,0 +1,233 @@ +#include "inmost_model.h" + + +namespace INMOST +{ + + void Model::AddEntry(std::string name, AbstractEntry & entry) + { + assert( !isInitialized() ); // do not add new data after initialization + assert( GetEntry(name) == NULL ); // do not overwrite entry with the same name + Entries.push_back(std::make_pair(name,&entry)); + } + + void Model::AddMesh(std::string name, Mesh & m) + { + assert( !isInitialized() ); // do not add new data after initialization + assert( GetMesh(name) == NULL ); // do not overwrite mesh with the same name + Meshes.push_back(std::make_pair(name,&m)); + } + + void Model::AddSubModel(std::string name, AbstractSubModel & submodel) + { + assert( !isInitialized() ); // do not add new data after initialization + assert( GetSubModel(name) == NULL ); // do not overwrite submodel with the same name + SubModels.push_back(std::make_pair(name,&submodel)); + } + + void Model::AddOperator(std::string name, AbstractOperator & op) + { + assert( !isInitialized() ); // do not add new data after initialization + assert( GetOperator(name) == NULL ); // do not overwrite submodel with the same name + Operators.push_back(std::make_pair(name,&op)); + } + + Mesh * Model::GetMesh(std::string name) + { + for(std::vector< std::pair >::iterator it = Meshes.begin(); + it != Meshes.end(); ++it) + if( it->first == name ) + return it->second; + return NULL; + } + + const Mesh * Model::GetMesh(std::string name) const + { + for(std::vector< std::pair >::const_iterator it = Meshes.begin(); + it != Meshes.end(); ++it) + if( it->first == name ) + return it->second; + return NULL; + } + + std::vector Model::GetMeshesNames() const + { + std::vector ret; + for(std::vector< std::pair >::const_iterator it = Meshes.begin(); + it != Meshes.end(); ++it) + ret.push_back(it->first); + return ret; + } + + AbstractEntry * Model::GetEntry(std::string name) + { + for(std::vector< std::pair >::iterator it = Entries.begin(); + it != Entries.end(); ++it) + if( it->first == name ) + return it->second; + return NULL; + } + + const AbstractEntry * Model::GetEntry(std::string name) const + { + for(std::vector< std::pair >::const_iterator it = Entries.begin(); + it != Entries.end(); ++it) + if( it->first == name ) + return it->second; + return NULL; + } + + std::vector Model::GetEntriesNames() const + { + std::vector ret; + for(std::vector< std::pair >::const_iterator it = Entries.begin(); + it != Entries.end(); ++it) + ret.push_back(it->first); + return ret; + } + + AbstractSubModel * Model::GetSubModel(std::string name) + { + for(std::vector< std::pair >::iterator it = SubModels.begin(); + it != SubModels.end(); ++it) + if( it->first == name ) + return it->second; + return NULL; + } + + const AbstractSubModel * Model::GetSubModel(std::string name) const + { + for(std::vector< std::pair >::const_iterator it = SubModels.begin(); + it != SubModels.end(); ++it) + if( it->first == name ) + return it->second; + return NULL; + } + + std::vector Model::GetSubModelsNames() const + { + std::vector ret; + for(std::vector< std::pair >::const_iterator it = SubModels.begin(); + it != SubModels.end(); ++it) + ret.push_back(it->first); + return ret; + } + + AbstractOperator * Model::GetOperator(std::string name) + { + for(std::vector< std::pair >::iterator it = Operators.begin(); + it != Operators.end(); ++it) + if( it->first == name ) + return it->second; + return NULL; + } + + const AbstractOperator * Model::GetOperator(std::string name) const + { + for(std::vector< std::pair >::const_iterator it = Operators.begin(); + it != Operators.end(); ++it) + if( it->first == name ) + return it->second; + return NULL; + } + + std::vector Model::GetOperatorsNames() const + { + std::vector ret; + for(std::vector< std::pair >::const_iterator it = Operators.begin(); + it != Operators.end(); ++it) + ret.push_back(it->first); + return ret; + } + + bool Model::PrepareEntries() + { + bool success = true; + //first initialize submodels + for(std::vector< std::pair >::iterator it = SubModels.begin(); + it != SubModels.end(); ++it) + success &= it->second->PrepareEntries(*this); + //second initialize operators + for(std::vector< std::pair >::iterator it = Operators.begin(); + it != Operators.end(); ++it) + success &= it->second->PrepareEntries(*this); + return success; + } + + bool Model::Initialize() + { + bool success = true; + //first initialize submodels + for(std::vector< std::pair >::iterator it = SubModels.begin(); + it != SubModels.end(); ++it) + success &= it->second->Initialize(*this); + //second initialize operators + for(std::vector< std::pair >::iterator it = Operators.begin(); + it != Operators.end(); ++it) + success &= it->second->Initialize(*this); + //third register all the entries + for(std::vector< std::pair >::iterator it = Entries.begin(); + it != Entries.end(); ++it) + { + it->second->reg_index = aut.RegisterEntry(*it->second); + it->second->SetOffsetTag(aut.GetEntry(it->second->reg_index).GetOffsetTag()); + } + initialized = success; + return success; + } + + bool Model::FillResidual(Residual & R) const + { + bool success = true; + R.Clear(); + Automatizator::MakeCurrent(&aut); + for(std::vector< std::pair >::const_iterator it = SubModels.begin(); + it != SubModels.end(); ++it) + success &= it->second->FillResidual(R); + Automatizator::RemoveCurrent(); + return success; + } + + bool Model::UpdateSolution(const Sparse::Vector & sol) + { + bool success = true; + for(std::vector< std::pair >::const_iterator it = SubModels.begin(); + it != SubModels.end(); ++it) + success &= it->second->UpdateSolution(sol); + for(std::vector< std::pair >::const_iterator it = Operators.begin(); + it != Operators.end(); ++it) + success &= it->second->UpdateSolution(sol); + return success; + } + + bool Model::UpdateTimeStep() + { + bool success = true; + for(std::vector< std::pair >::const_iterator it = SubModels.begin(); + it != SubModels.end(); ++it) + success &= it->second->UpdateTimeStep(); + for(std::vector< std::pair >::const_iterator it = Operators.begin(); + it != Operators.end(); ++it) + success &= it->second->UpdateTimeStep(); + return success; + } + + bool Model::SetTimeStep(double dt) + { + bool success = true; + for(std::vector< std::pair >::const_iterator it = SubModels.begin(); + it != SubModels.end(); ++it) + success &= it->second->SetTimeStep(dt); + return success; + + } + + bool Model::RestoreTimeStep() + { + bool success = true; + for(std::vector< std::pair >::const_iterator it = SubModels.begin(); + it != SubModels.end(); ++it) + success &= it->second->RestoreTimeStep(); + return success; + } +} diff --git a/Source/Autodiff/operator.cpp b/Source/Autodiff/operator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/Source/Autodiff/residual.cpp b/Source/Autodiff/residual.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4c8225cfce2c53d0fe78ed99b48228066f47a924 --- /dev/null +++ b/Source/Autodiff/residual.cpp @@ -0,0 +1,115 @@ +#include "inmost.h" + +#if defined(USE_AUTODIFF) + +namespace INMOST +{ +#if defined(USE_SOLVER) + void Residual::GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const + { + start = residual.GetFirstIndex(); + end = residual.GetLastIndex(); + } + void Residual::SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end) + { + jacobian.SetInterval(beg,end); + residual.SetInterval(beg,end); + } + void Residual::ClearResidual() + { + for(Sparse::Vector::iterator it = residual.Begin(); it != residual.End(); ++it) (*it) = 0.0; + } + void Residual::ClearJacobian() + { + for(Sparse::Matrix::iterator it = jacobian.Begin(); it != jacobian.End(); ++it) it->Clear(); + } + void Residual::ClearHessian() + { + for(Sparse::HessianMatrix::iterator it = hessian.Begin(); it != hessian.End(); ++it) it->Clear(); + } + void Residual::Clear() + { +#if defined(USE_OMP) +#pragma omp for +#endif //USE_OMP + for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) + { + residual[k] = 0.0; + if( !jacobian.Empty() ) jacobian[k].Clear(); + if( !hessian.Empty() ) hessian[k].Clear(); + } + } + INMOST_DATA_REAL_TYPE Residual::Norm() + { + INMOST_DATA_REAL_TYPE ret = 0; +#if defined(USE_OMP) +#pragma omp parallel for reduction(+:ret) +#endif //USE_OMP + for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) + ret += residual[k]*residual[k]; +#if defined(USE_MPI) + INMOST_DATA_REAL_TYPE tmp = ret; + MPI_Allreduce(&tmp, &ret, 1, INMOST_MPI_DATA_REAL_TYPE, MPI_SUM, jacobian.GetCommunicator()); +#endif + return sqrt(ret); + } + Residual & Residual::operator =(Residual const & other) + { + hessian = other.hessian; + jacobian = other.jacobian; + residual = other.residual; + return *this; + } + void Residual::Rescale(INMOST_DATA_ENUM_TYPE p) + { +#if defined(USE_OMP) +#pragma omp parallel for +#endif //USE_OMP + for(int k = (int)GetFirstIndex(); k < (int)GetLastIndex(); ++k) + { + INMOST_DATA_REAL_TYPE norm = 0.0; + if( p == ENUMUNDEF ) //infinite norm + { + for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q) + if( norm < fabs(jacobian[k].GetValue(q)) ) + norm = fabs(jacobian[k].GetValue(q)); + } + else //p-norm + { + for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q) + norm += pow(fabs(jacobian[k].GetValue(q)),p); + norm = pow(norm,1.0/p); + } + if( norm ) + { + norm = 1.0/norm; + residual[k] *= norm; + for(INMOST_DATA_ENUM_TYPE q = 0; q < jacobian[k].Size(); ++q) + jacobian[k].GetValue(q) *= norm; + if( !hessian.Empty() ) + for(INMOST_DATA_ENUM_TYPE q = 0; q < hessian[k].Size(); ++q) + hessian[k].GetValue(q) *= norm; + } + } + } + Residual::Residual(std::string name, INMOST_DATA_ENUM_TYPE start, INMOST_DATA_ENUM_TYPE end, INMOST_MPI_Comm _comm) + : jacobian(name,start,end,_comm),residual(name,start,end,_comm), hessian(name,0,0,_comm) + { + } + Residual::Residual(const Residual & other) + : jacobian(other.jacobian), residual(other.residual), hessian(other.hessian) + { + } + + Matrix Residual::operator [](const AbstractMatrix & rows) + { + Matrix ret(rows.Rows(),rows.Cols()); + for(INMOST_DATA_ENUM_TYPE i = 0; i < rows.Rows(); ++i) + for(INMOST_DATA_ENUM_TYPE j = 0; j < rows.Cols(); ++j) + new (&ret(i,j)) multivar_expression_reference(residual[rows(i,j)],&jacobian[rows(i,j)]); + return ret; + } +#endif //USE_SOLVER +} //namespace INMOST + +#endif // USE_AUTODIFF diff --git a/Source/Data/tag.cpp b/Source/Data/tag.cpp index aac3e5b10dfa96c711f865199dba61d3b227db32..c913840d5da63f6c7a446ba1cf67a2544a4604cc 100644 --- a/Source/Data/tag.cpp +++ b/Source/Data/tag.cpp @@ -35,6 +35,22 @@ namespace INMOST } return 0; } + + __INLINE static INMOST_DATA_ENUM_TYPE DataTypePackedBytesSize(DataType t) + { + switch(t) + { + case DATA_BULK: return sizeof(INMOST_DATA_BULK_TYPE); + case DATA_INTEGER: return sizeof(INMOST_DATA_INTEGER_TYPE); + case DATA_REAL: return sizeof(INMOST_DATA_REAL_TYPE); + case DATA_REMOTE_REFERENCE: throw -1; //todo, exchange of this data type is not yet supported + case DATA_REFERENCE: throw -1; //todo, exchange of this data type is not yet supported +#if defined(USE_AUTODIFF) + case DATA_VARIABLE: return sizeof(Sparse::Row::entry); +#endif + } + return 0; + } __INLINE static INMOST_DATA_ENUM_TYPE VariableDataSize(DataType t) @@ -52,6 +68,11 @@ namespace INMOST } return 0; } + + INMOST_DATA_ENUM_TYPE Tag::GetPackedBytesSize() const + { + return DataTypePackedBytesSize(GetDataType()); + } Tag::~Tag() { diff --git a/Source/Headers/CMakeLists.txt b/Source/Headers/CMakeLists.txt index d71a85dccf0814e5bef8783d6fe0f914035584cb..97f18a454730d6cc96007e1390ccc9cf38a0dbc6 100644 --- a/Source/Headers/CMakeLists.txt +++ b/Source/Headers/CMakeLists.txt @@ -14,8 +14,11 @@ set(HEADER ${CMAKE_CURRENT_SOURCE_DIR}/inmost_solver.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_partitioner.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_autodiff.h + ${CMAKE_CURRENT_SOURCE_DIR}/inmost_residual.h + ${CMAKE_CURRENT_SOURCE_DIR}/inmost_model.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_expression.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_variable.h + ${CMAKE_CURRENT_SOURCE_DIR}/inmost_block_variable.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_sparse.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_nonlinear.h ${CMAKE_CURRENT_SOURCE_DIR}/inmost_xml.h diff --git a/Source/Headers/container.hpp b/Source/Headers/container.hpp index 273f83d352b92f588b65807d25419a1e38acaaca..17754b9d372ea0fbed3173929a24c0f89febf80e 100644 --- a/Source/Headers/container.hpp +++ b/Source/Headers/container.hpp @@ -2721,6 +2721,13 @@ namespace INMOST m_size = n; } }; + + /// This class is used to replace #pragma omp threadprivate + /// Functionality that is not supported on many older systems. + template + class thread_private + { + }; } #endif diff --git a/Source/Headers/inmost.h b/Source/Headers/inmost.h index e5fc25342ca339f6921cec9c8eccc7617534e1fb..e4c38ee239216d1f0da78491e8e6806fb7abb95f 100644 --- a/Source/Headers/inmost.h +++ b/Source/Headers/inmost.h @@ -5,10 +5,17 @@ #include "inmost_common.h" #include "inmost_mesh.h" #include "inmost_dense.h" +#include "inmost_sparse.h" #include "inmost_solver.h" #include "inmost_partitioner.h" +#include "inmost_expression.h" #include "inmost_variable.h" +#include "inmost_block_variable.h" +#include "inmost_autodiff.h" #include "inmost_nonlinear.h" +#include "inmost_residual.h" +#include "inmost_operator.h" +#include "inmost_model.h" #include "inmost_xml.h" #endif // INMOST_H_INCLUDED diff --git a/Source/Headers/inmost_autodiff.h b/Source/Headers/inmost_autodiff.h index 83d3335477a33a9f444d9900abd0a7c382298c41..8054ac8406e7d9617c16ec4ad990622fa20e5d21 100644 --- a/Source/Headers/inmost_autodiff.h +++ b/Source/Headers/inmost_autodiff.h @@ -1,119 +1,23 @@ #ifndef INMOST_AUTODIFF_H_INCLUDED #define INMOST_AUTODIFF_H_INCLUDED + #include "inmost_common.h" #include "inmost_mesh.h" -#include "inmost_solver.h" -#include "inmost_variable.h" -#if defined(USE_AUTODIFF) -#include +#if defined(USE_AUTODIFF) && defined(USE_MESH) namespace INMOST { class Automatizator; //forward declaration - -#if defined(USE_OMP) -#define OMP_THREAD omp_get_thread_num() -#define MAX_THREADS omp_get_max_threads() -#else //USE_OMP -#define OMP_THREAD 0 -#define MAX_THREADS 1 -#endif //USE_OMP - -#if defined(USE_SOLVER) - /// The Residual class provides a representation for array of residuals of nonlinear equations. - /// By working with the residual class you automatically assemble right hand side and - /// the jacobian of a nonlinear system of equation. - /// Jacobian matrix has a sparse representation. - /// \todo - /// 1. Extend for hessian calculation. - class Residual - { - Sparse::HessianMatrix hessian; ///< Hessian matrix - Sparse::Matrix jacobian; ///< Jacobian matrix. - Sparse::Vector residual; ///< Right hand side vector. - Sparse::LockService locks; ///< Array of locks for openmp shared access. - public: - /// Constructor - /// @param name Name for the matrix and right hand side. Can be used to set options for linear solver. - /// @param start First index of the equation in the local partition. Use Automatizator::GetFirstIndex. - /// @param end Last index of the equation in the local partition. Use Automatizator::GetLastIndex. - /// @param _comm MPI Communicator. - Residual(std::string name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD); - /// Copy constructor. - /// \warning May be expensive if matrices are large. - Residual(const Residual & other); - /// Assignment operator. - /// \warning May be expensive if matrices are large. - Residual & operator =(Residual const & other); - /// Retrive the first index of the equations in the local partition. - INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return residual.GetFirstIndex();} - /// Retrive the last index of the equations in the local partition. - INMOST_DATA_ENUM_TYPE GetLastIndex() const {return residual.GetLastIndex();} - /// Retrive the first and the last indices of the equations in the local partition - /// @param start The first index of the equations will be recorded here. - /// @param end The last index of the equations will be recorded here. - void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const; - /// Assign the new first and last indices of the equations in the local partition. - /// @param start The new first index of the equations. - /// @param end The new last index of the equations. - void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end); - /// Retrive a residual value and a jacobian row corresponding to certain equation. - /// @param row Equation number. - /// @return A structure that can be used in or assigned an automatic differentiation expression. - __INLINE multivar_expression_reference operator [](INMOST_DATA_ENUM_TYPE row) - {return multivar_expression_reference(residual[row],&jacobian[row]);} - /// Retrive a vector of entries in residual, corresponding to a set of equations. - /// @param rows A row-vector of equation numbers. - /// @param A structure that can be used in or assigned an automatic differentiation matrix expression. - Matrix operator [](const AbstractMatrix & rows); - /// Retrive hessian matrix. Use in nonlinear solver. - Sparse::HessianMatrix & GetHessian() {return hessian;} - /// Retrive hessian matrix without right of modificaiton. - const Sparse::HessianMatrix & GetHessian() const {return hessian;} - /// Retrive jacobian matrix. Use in Sparse::Solver::Solve function. - Sparse::Matrix & GetJacobian() {return jacobian;} - /// Retrive jacobian matrix without right of modificaiton. - const Sparse::Matrix & GetJacobian() const {return jacobian;} - /// Retrive right hand side vector. Use in Sparse::Solver::Solve function. - Sparse::Vector & GetResidual() {return residual;} - /// Retrive right hand side vector without right of modification. - const Sparse::Vector & GetResidual() const {return residual;} - /// Zero out right hand side vector. - void ClearResidual(); - /// Remove all entries in jacobian matrix. - void ClearJacobian(); - /// Remove all entries in hessian matrix. - void ClearHessian(); - /// Zero out right hand side vector and remove all entries in jacobian matrix. - void Clear(); - /// Compute the second norm of the right hand side vector over all of the processors. - INMOST_DATA_REAL_TYPE Norm(); - /// Normalize jacobian rows to unit p-norms and scale right hand side accordingly. - /// Use ENUMUNDEF as p to scale to infinite-norm. - void Rescale(INMOST_DATA_ENUM_TYPE p = 2); - /// Initialize openmp locks. - void InitLocks() {locks.SetInterval(GetFirstIndex(),GetLastIndex());} - /// Lock an equation to avoid simultaneous shared access. - /// @param pos Equation number. - __INLINE void Lock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.Lock(pos);} - /// UnLock an equation to allow simultaneous shared access. - /// @param pos Equation number. - __INLINE void UnLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.UnLock(pos);} - /// Try to lock the equation. - /// @param pos Equation number. - /// @return True if equation was locked. - __INLINE bool TestLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) return locks.TestLock(pos); return false;} - }; -#endif //USE_SOLVER -#if defined(USE_MESH) /// This class is used to organize unknowns in abstract way, /// it should be registered with and managed by class Automatizator. /// \todo /// 1. Is there a need for layout on how matrices are returned? /// 2. Is there a need for layout on how unknowns and equations are arrenged? + /// 3. Function for update of variables. + /// 4. Function for synchronization of variables. class AbstractEntry { INMOST_DATA_ENUM_TYPE reg_index; ///< Index of block registry with Automatizator. @@ -121,7 +25,8 @@ namespace INMOST MarkerType mask; ///< Allows to enable or disable the entire block. ElementType etype; ///< Type of elements on which unknowns are defined. public: - AbstractEntry(ElementType etype = NONE, MarkerType mask = 0) : reg_index(ENUMUNDEF), offset_tag(Tag()), mask(mask), etype(etype) {} + AbstractEntry(ElementType etype = NONE, MarkerType mask = 0) + : reg_index(ENUMUNDEF), offset_tag(Tag()), mask(mask), etype(etype) {} /// Retrive mask of the block. MarkerType GetMask() const {return mask;} /// Set mask for the block. @@ -135,7 +40,7 @@ namespace INMOST /// Retrive tag that stores enumeration offset on each element. void SetOffsetTag(TagInteger tag) {offset_tag = tag;} /// Check that the block is valid on given element. - bool isValid(const Storage & e) const {return (e.GetElementType() & etype) && (mask == 0 || e->GetMarker(mask));} + bool isValid(const Storage & e) const {return reg_index != ENUMUNDEF && (e.GetElementType() & etype) && (mask == 0 || e->GetMarker(mask));} /// Return value in vector of unknowns of the block at certain position. /// @param pos Position for which to extract the value, should be no larger the MatrixSize. virtual INMOST_DATA_REAL_TYPE Value(const Storage & e, INMOST_DATA_ENUM_TYPE pos) const = 0; @@ -154,6 +59,8 @@ namespace INMOST /// Return vector filled with indices of unknowns of the block. virtual iMatrix Index(const Storage & e) const = 0; /// Return vector filled with unknowns of the block with their derivatives. + virtual uMatrix Unknown(const Storage & e) const = 0; + /// Return vector filled with unknowns of the block with their derivatives. virtual uMatrix operator [](const Storage & e) const = 0; /// The intended size of the matrix for this entry. virtual INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const = 0; @@ -175,10 +82,14 @@ namespace INMOST virtual AbstractEntry * Copy() const = 0; /// Retrive a registration index. INMOST_DATA_ENUM_TYPE GetRegistrationIndex() const {return reg_index;} + /// Update variables contained in block on ghost elements of the grid. + /// For synchronization of data in all blocks see Automatizator::SynhronizeData. + void SynchronizeData(); /// Destructor. virtual ~AbstractEntry() {} friend class Automatizator; //provide registration index from inside of Automatizator + friend class Model; //provide registration index from inside of Model }; /// This class is used to organize unknowns into blocks, /// blocks enumeration are managed by class Automatizator. @@ -206,6 +117,8 @@ namespace INMOST /// Return vector filled with indices of unknowns of the block. iMatrix Index(const Storage & e) const {iMatrix ret(MatrixSize(e),1); for(unsigned k = 0; k < Size(); ++k) ret(k,0) = Index(e,k); return ret; } /// Return vector filled with unknowns of the block with their derivatives. + uMatrix Unknown(const Storage & e) const {return BlockEntry::operator [](e);} + /// Return vector filled with unknowns of the block with their derivatives. uMatrix operator [](const Storage & e) const {uMatrix ret(MatrixSize(e),1); for(unsigned k = 0; k < Size(); ++k) ret(k,0) = Unknown(e,k); return ret; } /// The intended size of the matrix for this entry. INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const {return Size();} @@ -247,6 +160,8 @@ namespace INMOST /// Return vector filled with indices of unknowns of the block. iMatrix Index(const Storage & e) const { iMatrix ret(1,1); ret(0,0) = Index(e,0); return ret; } /// Return vector filled with unknowns of the block with their derivatives. + uMatrix Unknown(const Storage & e) const {return SingleEntry::operator [](e);} + /// Return vector filled with unknowns of the block with their derivatives. uMatrix operator [](const Storage & e) const { uMatrix ret(1,1); ret(0,0) = Unknown(e,0); return ret; } /// The intended size of the matrix for this entry. INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const {return 1;} @@ -287,6 +202,8 @@ namespace INMOST /// Return vector filled with indices of unknowns of the block. iMatrix Index(const Storage & e) const { iMatrix ret(MatrixSize(e),1); for(int k = 0; k < (int)unknown_tag[e].size(); ++k) ret(k,0) = Index(e,k); return ret; } /// Return vector filled with unknowns of the block with their derivatives. + uMatrix Unknown(const Storage & e) const {return VectorEntry::operator [](e);} + /// Return vector filled with unknowns of the block with their derivatives. uMatrix operator [](const Storage & e) const { uMatrix ret(MatrixSize(e),1); for(int k = 0; k < (int)unknown_tag[e].size(); ++k) ret(0,0) = Unknown(e,k); return ret; } /// The intended size of the matrix for this entry. INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const {return (INMOST_DATA_ENUM_TYPE)unknown_tag[e].size();} @@ -341,6 +258,8 @@ namespace INMOST /// Return vector filled with indices of unknowns of the block. iMatrix Index(const Storage & e) const {iMatrix ret(MatrixSize(e),1); for(INMOST_DATA_ENUM_TYPE k = 0; k < Size(); ++k) ret(k,0) = Index(e,k); return ret; } /// Return vector filled with unknowns of the block with their derivatives. + uMatrix Unknown(const Storage & e) const {return StatusBlockEntry::operator [](e);} + /// Return vector filled with unknowns of the block with their derivatives. uMatrix operator [](const Storage & e) const {uMatrix ret(MatrixSize(e),1); for(INMOST_DATA_ENUM_TYPE k = 0; k < Size(); ++k) ret(k,0) = Unknown(e,k); return ret; } /// The intended size of the matrix for this entry. INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const {return Size();} @@ -392,6 +311,8 @@ namespace INMOST /// Return vector filled with indices of unknowns of the block. iMatrix Index(const Storage & e) const; /// Return vector filled with unknowns of the block with their derivatives. + uMatrix Unknown(const Storage & e) const {return MultiEntry::operator [](e);} + /// Return vector filled with unknowns of the block with their derivatives. uMatrix operator [](const Storage & e) const; /// The intended size of the matrix for this entry. INMOST_DATA_ENUM_TYPE MatrixSize(const Storage & e) const; @@ -474,11 +395,21 @@ namespace INMOST /// \warning /// 1. Have to call Automatizator::EnumerateEntries to recompute indices. void UnregisterEntry(INMOST_DATA_ENUM_TYPE ind); + /// Swith a registered tag to be non-active, in this case it's unknowns are considered to be constant. + /// @param ind Integer returned from Automatizator::RegisterTag. + /// \warning + /// 1. Have to call Automatizator::EnumerateEntries to recompute indices. + void DeactivateEntry(INMOST_DATA_ENUM_TYPE ind); + /// Swith a registered tag to be active, in this case it's unknowns are considered to be variable. + /// @param ind Integer returned from Automatizator::RegisterTag. + /// \warning + /// 1. Have to call Automatizator::EnumerateEntries to recompute indices. + void ActivateEntry(INMOST_DATA_ENUM_TYPE ind); /// Set index for every data entry of dynamic tag. void EnumerateEntries(); /// Check whether the tag is still registered. /// @param True if tag is still registered. - __INLINE bool isRegisteredEntry(INMOST_DATA_ENUM_TYPE ind) const {return act_blocks[ind];} + __INLINE bool isRegisteredEntry(INMOST_DATA_ENUM_TYPE ind) const {return reg_blocks[ind] != NULL;} /// Get index of the unknown associated with the entry on element. INMOST_DATA_ENUM_TYPE GetIndex(const Storage & e, INMOST_DATA_ENUM_TYPE reg_index, INMOST_DATA_ENUM_TYPE pos = 0) const {return GetEntry(reg_index).Index(e,pos);} /// Get value of the unknown associated with the entry on element. @@ -503,9 +434,12 @@ namespace INMOST /// Lists all the indices of registered tags. /// @return An array with indices corresponding to all registered tags. std::vector ListRegisteredEntries() const; + /// Update variables contained in all block of automatizator on ghost elements of the grid. + /// For synchronization of data in individual blocks see AbstractEntry::SynhronizeData. + void SynchronizeData(); }; -#endif //USE_MESH } //namespace INMOST -#endif //USE_AUTODIFF +#endif //USE_AUTODIFF && USE_MESH + #endif //INMOST_AUTODIFF_H_INCLUDED diff --git a/Source/Headers/inmost_block_variable.h b/Source/Headers/inmost_block_variable.h new file mode 100644 index 0000000000000000000000000000000000000000..0803a3bd91967d37ea9e2686dc7909bb12a1ab3e --- /dev/null +++ b/Source/Headers/inmost_block_variable.h @@ -0,0 +1,880 @@ + +#ifndef INMOST_AUTODIFF_ETBVAR_H_INCLUDED +#define INMOST_AUTODIFF_ETBVAR_H_INCLUDED +#include "inmost_common.h" +#include "inmost_expression.h" +#include "inmost_mesh.h" +#include "inmost_autodiff.h" +#include "inmost_solver.h" +#include "inmost_variable.h" +#include //for debug +#include + +#if defined(USE_AUTODIFF) && defined(USE_MESH) + +//TODO: +// 1. Add operations ConcatRows, ConcatCols that help assemble matrix from pieces + + +//This should stop Visual Studio from complaining of very long auto-generated class types +#ifdef _MSC_VER +#pragma warning(disable : 4503) +#endif + + + +namespace INMOST +{ + + class abstract_dynamic_block_variable + { + public: + virtual rMatrix Value (const Storage & e) const = 0; + virtual vMatrix Variable(const Storage & e) const = 0; + virtual abstract_dynamic_block_variable * Copy() const = 0; + virtual ~abstract_dynamic_block_variable() {} + }; + + template + class get_block_variable + { + public: + virtual RetType operator()(const Storage & e) const = 0; + }; + + template<> + class get_block_variable + { + const abstract_dynamic_block_variable & var; + public: + typedef vMatrix type; + get_block_variable(const abstract_dynamic_block_variable & var) : var(var) {} + vMatrix operator()(const Storage & e) const {return var.Variable(e);} + }; + + template<> + class get_block_variable + { + const abstract_dynamic_block_variable & var; + public: + typedef rMatrix type; + get_block_variable(const abstract_dynamic_block_variable & var) : var(var) {} + rMatrix operator()(const Storage & e) const {return var.Value(e);} + }; + + + class stored_block_variable_expression : public abstract_dynamic_block_variable + { + abstract_dynamic_block_variable * var; + public: + stored_block_variable_expression() : var(NULL) {} + stored_block_variable_expression(const abstract_dynamic_block_variable & pvar) : var(pvar.Copy()) {} + stored_block_variable_expression(const stored_block_variable_expression & other) : var(other.var->Copy()) {} + ~stored_block_variable_expression() {delete var; var = NULL;} + stored_block_variable_expression operator =(stored_block_variable_expression const & other) {var = other.var->Copy(); return *this;} + stored_block_variable_expression operator =(const abstract_dynamic_block_variable & pvar) {var = pvar.Copy(); return *this;} + rMatrix Value(const Storage & e) const {return var->Value(e);} + vMatrix Variable(const Storage & e) const {return var->Variable(e);} + + template + get_variable get_variable() {return get_variable(*var);} + abstract_dynamic_block_variable & retrive_expression() {return *var;} + const abstract_dynamic_block_variable & retrive_expression() const {return *var;} + abstract_dynamic_block_variable * Copy() const {return static_cast(new stored_block_variable_expression(*this));} + /// Checks that the stored expresison was defined. + bool isDefined() const {return var != NULL;} + }; + + + class dynamic_block_variable : public abstract_dynamic_block_variable + { + private: + const AbstractEntry * entry; + public: + dynamic_block_variable() :entry(NULL) {} + dynamic_block_variable(Automatizator & aut, INMOST_DATA_ENUM_TYPE reg_index) : entry(reg_index==ENUMUNDEF?NULL:&aut.GetEntry(reg_index)) {} + dynamic_block_variable(const AbstractEntry * re) : entry(re) {} + dynamic_block_variable(const dynamic_block_variable & other) : entry(other.entry) {} + dynamic_block_variable & operator =(const dynamic_block_variable & other) + { + entry = other.entry; + return * this; + } + rMatrix Value(const Storage & e) const {return entry->Value(e);} + //iMatrix Index(const Storage & e) const {return entry->isValid(e) ? entry->Index(e):iMatrix(entry->MatrixSize(e),1,ENUMUNDEF);} + vMatrix Variable(const Storage & e) const {return entry->Unknown(e);} + bool isUnknown(const Storage & e) const {return entry->isValid(e)?true:false;} + abstract_dynamic_block_variable * Copy() const {return static_cast(new dynamic_block_variable(*this));} + }; + + class const_block_variable : public abstract_dynamic_block_variable + { + private: + rMatrix value; + public: + const_block_variable(INMOST_DATA_REAL_TYPE _value) : value(_value) {} + const_block_variable(const const_block_variable & other) : value(other.value) {} + const_block_variable & operator =(const const_block_variable & other) + { + value = other.value; + return * this; + } + rMatrix Value(const Storage & e) const {return value;} + vMatrix Variable(const Storage & e) const { return value; } + abstract_dynamic_block_variable * Copy() const {return static_cast(new const_block_variable(*this));} + }; + + class static_block_variable : public abstract_dynamic_block_variable + { + private: + TagRealArray value_tag; + INMOST_DATA_ENUM_TYPE n,m; + public: + static_block_variable(TagRealArray t) + : value_tag(t), n(t.GetSize()), m(1) {assert(t.GetDataType() == DATA_REAL);} + static_block_variable(TagRealArray t, INMOST_DATA_ENUM_TYPE pn, INMOST_DATA_ENUM_TYPE pm) + : value_tag(t), n(pn), m(pm) {assert(t.GetDataType() == DATA_REAL);} + static_block_variable(const static_block_variable & other) + : value_tag(other.value_tag), n(other.n), m(other.m) {} + static_block_variable & operator =(const static_block_variable & other) + { + value_tag = other.value_tag; + n = other.n; + m = other.m; + return * this; + } + rMatrix Value(const Storage & e) const {return value_tag(e,n,m);} + vMatrix Variable(const Storage & e) const {return value_tag(e,n,m);} + TagRealArray ValueTag() {return value_tag;} + bool isUnknown(const Storage & e) const {(void)e; return false;} + abstract_dynamic_block_variable * Copy() const {return static_cast(new static_block_variable(*this));} + }; + + class stored_block_variable : public abstract_dynamic_block_variable + { + private: + TagVariableArray variable_tag; + INMOST_DATA_ENUM_TYPE n,m; + public: + stored_block_variable(TagVariableArray t) + : variable_tag(t), n(t.GetSize()), m(1) {} + stored_block_variable(TagVariableArray t, INMOST_DATA_ENUM_TYPE pn, INMOST_DATA_ENUM_TYPE pm) + : variable_tag(t), n(pn), m(pm) {assert(t.GetDataType() == DATA_VARIABLE);} + stored_block_variable(const stored_block_variable & other) + : variable_tag(other.variable_tag), n(other.n), m(other.m) {} + stored_block_variable & operator =(const stored_block_variable & other) + { + variable_tag = other.variable_tag; + n = other.n; + m = other.m; + return * this; + } + rMatrix Value(const Storage & e) const { return variable_tag(e,n,m);} + vMatrix Variable(const Storage & e) const { return variable_tag(e,n,m);} + TagVariableArray VariableTag() {return variable_tag;} + bool isUnknown(const Storage & e) const {(void)e; return false;} + abstract_dynamic_block_variable * Copy() const {return static_cast(new stored_block_variable(*this));} + }; + + + /// \todo coefficients could be matrices here, introduce another class? + class stencil_block_variable : public abstract_dynamic_block_variable + { + private: + TagReferenceArray tag_elems; + TagRealArray tag_coefs; + stored_block_variable_expression Arg; + public: + stencil_block_variable(Tag tag_elems, Tag tag_coefs, const abstract_dynamic_block_variable & parg) + : tag_elems(tag_elems), tag_coefs(tag_coefs), Arg(parg) {} + stencil_block_variable(const stencil_block_variable & other) + : tag_elems(other.tag_elems), tag_coefs(other.tag_coefs), Arg(other.Arg) {} + stencil_block_variable & operator =(const stencil_block_variable & other) + { + tag_elems = other.tag_elems; + tag_coefs = other.tag_coefs; + Arg = other.Arg; + return * this; + } + rMatrix Value(const Storage & e) const + { + Storage::real_array coefs = tag_coefs[e]; + Storage::reference_array elems = tag_elems[e]; + assert(coefs.size() == elems.size()); + rMatrix ret = coefs[0]*Arg.Value(elems[0]); + for(INMOST_DATA_ENUM_TYPE k = 1; k < elems.size(); ++k) + ret += coefs[k]*Arg.Value(elems[k]); + return ret; + } + vMatrix Variable(const Storage & e) const + { + Storage::real_array coefs = tag_coefs[e]; + Storage::reference_array elems = tag_elems[e]; + assert(coefs.size() == elems.size()); + vMatrix ret = coefs[0]*Arg.Variable(elems[0]); + for(INMOST_DATA_ENUM_TYPE k = 1; k < elems.size(); ++k) + ret += coefs[k]*Arg.Variable(elems[k]); + return ret; + } + abstract_dynamic_block_variable * Copy() const {return static_cast(new stencil_block_variable(*this));} + }; + + ///Apply table component-wise on argument matrix. + class table_block_variable : public abstract_dynamic_block_variable + { + stored_block_variable_expression Arg; + keyval_table Table; + public: + table_block_variable(const abstract_dynamic_block_variable & parg, const keyval_table & ptable) : Arg(parg), Table(ptable) {} + table_block_variable(const table_block_variable & other) : Arg(other.Arg), Table(other.Table) {} + table_block_variable & operator = (table_block_variable const & other) {Arg = other.Arg; Table = other.Table; return * this;} + rMatrix Value(const Storage & e) const + { + rMatrix ret = Arg.Value(e); + for(int k = 0; k < ret.Rows(); ++k) + for(int l = 0; l < ret.Cols(); ++l) + ret(k,l) = get_table(ret(k,l),Table); + return ret; + } + vMatrix Variable(const Storage & e) const + { + vMatrix ret = Arg.Variable(e); + for(int k = 0; k < ret.Rows(); ++k) + for(int l = 0; l < ret.Cols(); ++l) + ret(k,l) = get_table(ret(k,l),Table); + return ret; + } + abstract_dynamic_block_variable * Copy() const {return static_cast(new table_block_variable(*this));} + }; + + /// This class makes possible to evaluate different expressions on different element types. + /// See etype_branch function. + class etype_branch_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Variable expression to be evaluated when type of provided element matches selected types. + stored_block_variable_expression ArgB; //< Variable expression to be evaluated when type of provided element does not match selected types. + ElementType types_true; //< Selected types of elements. + public: + /// Constructor. Used by etype_branch function. + etype_branch_block_variable(ElementType _types_true, const abstract_dynamic_block_variable & _ArgA, const abstract_dynamic_block_variable & _ArgB) + : types_true(_types_true), ArgA(_ArgA), ArgB(_ArgB) {} + /// Copy constructor. + etype_branch_block_variable(const etype_branch_block_variable & other) + : types_true(other.types_true), ArgA(other.ArgA), ArgB(other.ArgB) {} + /// Assignment operator. + etype_branch_block_variable & operator =(etype_branch_block_variable const & other) + { + types_true = other.types_true; + ArgA = other.ArgA; + ArgB = other.ArgB; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return (e->GetElementType() & types_true) ? ArgA.Value(e) : ArgB.Value(e);} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return (e->GetElementType() & types_true) ? ArgA.Variable(e) : ArgB.Variable(e);} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new etype_branch_block_variable(*this));} + }; + + /// This class makes possible to evaluate different expressions depending on the markers. + /// Works similarly for shared and private markers. + /// See marker_branch function. + class marker_branch_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Variable expression to be evaluated when marker is set on the element. + stored_block_variable_expression ArgB; //< Variable expression to be evaluated when marker is not set on the element. + MarkerType marker; //< Marker. + public: + /// Constructor. Used by marker_branch function. + marker_branch_block_variable(MarkerType _marker, const abstract_dynamic_block_variable & _ArgA, const abstract_dynamic_block_variable & _ArgB) : marker(_marker), ArgA(_ArgA), ArgB(_ArgB) {} + /// Copy constructor. + marker_branch_block_variable(const marker_branch_block_variable & other) : marker(other.marker), ArgA(other.ArgA), ArgB(other.ArgB) {} + /// Assignment operator. + marker_branch_block_variable & operator =(marker_branch_block_variable const & other) + { + marker = other.marker; + ArgA = other.ArgA; + ArgB = other.ArgB; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ( isPrivate(marker) ? e->GetPrivateMarker(marker) : e->GetMarker(marker) ) ? ArgA.Value(e) : ArgB.Value(e);} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ( isPrivate(marker) ? e->GetPrivateMarker(marker) : e->GetMarker(marker) ) ? ArgA.Variable(e) : ArgB.Variable(e);} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new marker_branch_block_variable(*this));} + }; + + /// This class represents addition of two matrices. + class addition_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression on the left. + stored_block_variable_expression ArgB; //< Block variable expression on the right. + public: + /// Constructor. + addition_block_variable(const abstract_dynamic_block_variable & _ArgA, + const abstract_dynamic_block_variable & _ArgB) + : ArgA(_ArgA), ArgB(_ArgB) {} + /// Copy constructor. + addition_block_variable(const addition_block_variable & other) + : ArgA(other.ArgA), ArgB(other.ArgB) {} + /// Assignment operator. + addition_block_variable & operator =(addition_block_variable const & other) + { + ArgA = other.ArgA; + ArgB = other.ArgB; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e) + ArgB.Value(e);} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e) + ArgB.Variable(e);} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new addition_block_variable(*this));} + }; + + /// This class represents subtraction of two matrices. + class subtraction_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression on the left. + stored_block_variable_expression ArgB; //< Block variable expression on the right. + public: + /// Constructor. + subtraction_block_variable(const abstract_dynamic_block_variable & _ArgA, + const abstract_dynamic_block_variable & _ArgB) + : ArgA(_ArgA), ArgB(_ArgB) {} + /// Copy constructor. + subtraction_block_variable(const subtraction_block_variable & other) + : ArgA(other.ArgA), ArgB(other.ArgB) {} + /// Assignment operator. + subtraction_block_variable & operator =(subtraction_block_variable const & other) + { + ArgA = other.ArgA; + ArgB = other.ArgB; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e) - ArgB.Value(e);} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e) - ArgB.Variable(e);} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new subtraction_block_variable(*this));} + }; + + /// This class represents multiplication of two matrices. + class multiplication_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression on the left. + stored_block_variable_expression ArgB; //< Block variable expression on the right. + public: + /// Constructor. + multiplication_block_variable(const abstract_dynamic_block_variable & _ArgA, + const abstract_dynamic_block_variable & _ArgB) + : ArgA(_ArgA), ArgB(_ArgB) {} + /// Copy constructor. + multiplication_block_variable(const multiplication_block_variable & other) + : ArgA(other.ArgA), ArgB(other.ArgB) {} + /// Assignment operator. + multiplication_block_variable & operator =(multiplication_block_variable const & other) + { + ArgA = other.ArgA; + ArgB = other.ArgB; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e)*ArgB.Value(e);} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e)*ArgB.Variable(e);} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new multiplication_block_variable(*this));} + }; + + /// This class represents division of two matrices, this is technically B^{-1}A. + class division_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression on the left. + stored_block_variable_expression ArgB; //< Block variable expression on the right. + public: + /// Constructor. + division_block_variable(const abstract_dynamic_block_variable & _ArgA, + const abstract_dynamic_block_variable & _ArgB) + : ArgA(_ArgA), ArgB(_ArgB) {} + /// Copy constructor. + division_block_variable(const division_block_variable & other) + : ArgA(other.ArgA), ArgB(other.ArgB) {} + /// Assignment operator. + division_block_variable & operator =(division_block_variable const & other) + { + ArgA = other.ArgA; + ArgB = other.ArgB; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return (ArgA.Value(e)/ArgB.Value(e)).first;} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return (ArgA.Variable(e)/ArgB.Variable(e)).first;} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new division_block_variable(*this));} + }; + + /// This class represents division of two matrices, this is technically B^{-1}A. + class concat_cols_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression on the left. + stored_block_variable_expression ArgB; //< Block variable expression on the right. + public: + /// Constructor. + concat_cols_block_variable(const abstract_dynamic_block_variable & _ArgA, + const abstract_dynamic_block_variable & _ArgB) + : ArgA(_ArgA), ArgB(_ArgB) {} + /// Copy constructor. + concat_cols_block_variable(const concat_cols_block_variable & other) + : ArgA(other.ArgA), ArgB(other.ArgB) {} + /// Assignment operator. + concat_cols_block_variable & operator =(concat_cols_block_variable const & other) + { + ArgA = other.ArgA; + ArgB = other.ArgB; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e).ConcatCols(ArgB.Value(e));} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e).ConcatCols(ArgB.Variable(e));} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new concat_cols_block_variable(*this));} + }; + + /// This class represents division of two matrices, this is technically B^{-1}A. + class concat_rows_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression on the left. + stored_block_variable_expression ArgB; //< Block variable expression on the right. + public: + /// Constructor. + concat_rows_block_variable(const abstract_dynamic_block_variable & _ArgA, + const abstract_dynamic_block_variable & _ArgB) + : ArgA(_ArgA), ArgB(_ArgB) {} + /// Copy constructor. + concat_rows_block_variable(const concat_rows_block_variable & other) + : ArgA(other.ArgA), ArgB(other.ArgB) {} + /// Assignment operator. + concat_rows_block_variable & operator =(concat_rows_block_variable const & other) + { + ArgA = other.ArgA; + ArgB = other.ArgB; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e).ConcatRows(ArgB.Value(e));} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e).ConcatRows(ArgB.Variable(e));} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new concat_rows_block_variable(*this));} + }; + + + + /// This class represents transposition of the matrix, this is A^T. + class transpose_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression + public: + /// Constructor. + transpose_block_variable(const abstract_dynamic_block_variable & _ArgA) + : ArgA(_ArgA) {} + /// Copy constructor. + transpose_block_variable(const transpose_block_variable & other) + : ArgA(other.ArgA) {} + /// Assignment operator. + transpose_block_variable & operator =(transpose_block_variable const & other) + { + ArgA = other.ArgA; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e).Transpose();} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e).Transpose();} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new transpose_block_variable(*this));} + }; + + /// This class represents submatrix of the matrix. + class submatrix_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression + INMOST_DATA_ENUM_TYPE row1,row2,col1,col2; + public: + /// Constructor. + submatrix_block_variable(const abstract_dynamic_block_variable & _ArgA, + INMOST_DATA_ENUM_TYPE row1, + INMOST_DATA_ENUM_TYPE row2, + INMOST_DATA_ENUM_TYPE col1, + INMOST_DATA_ENUM_TYPE col2) + : ArgA(_ArgA), row1(row1), row2(row2), col1(col1), col2(col2) {} + /// Copy constructor. + submatrix_block_variable(const submatrix_block_variable & b) + : ArgA(b.ArgA), row1(b.row1), row2(b.row2), col1(b.col1), col2(b.col2) {} + /// Assignment operator. + submatrix_block_variable & operator =(submatrix_block_variable const & b) + { + row1 = b.row1; + row2 = b.row2; + col1 = b.col1; + col2 = b.col2; + ArgA = b.ArgA; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e)(row1,row2,col1,col2);} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e)(row1,row2,col1,col2);} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new submatrix_block_variable(*this));} + }; + + /// This class represents multiplication of the matrix by the constant. + class multiply_const_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression + INMOST_DATA_REAL_TYPE value; + public: + /// Constructor. + multiply_const_block_variable(const abstract_dynamic_block_variable & _ArgA, + INMOST_DATA_REAL_TYPE val) + : ArgA(_ArgA), value(val) {} + /// Copy constructor. + multiply_const_block_variable(const multiply_const_block_variable & b) + : ArgA(b.ArgA), value(b.value) {} + /// Assignment operator. + multiply_const_block_variable & operator =(multiply_const_block_variable const & b) + { + value = b.value; + ArgA = b.ArgA; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e)*value;} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e)*value;} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new multiply_const_block_variable(*this));} + }; + + /// This class represents multiplication of the matrix by the variable. + class condition_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression + stored_block_variable_expression ArgB; //< Block variable expression + stored_variable_expression ArgC; + public: + /// Constructor. + condition_block_variable(const abstract_dynamic_variable & _ArgC, + const abstract_dynamic_block_variable & _ArgA, + const abstract_dynamic_block_variable & _ArgB) + : ArgA(_ArgA), ArgB(_ArgB), ArgC(_ArgC) {} + /// Copy constructor. + condition_block_variable(const condition_block_variable & b) + : ArgA(b.ArgA), ArgB(b.ArgB), ArgC(b.ArgC) {} + /// Assignment operator. + condition_block_variable & operator =(condition_block_variable const & b) + { + ArgB = b.ArgB; + ArgA = b.ArgA; + ArgC = b.ArgC; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgC.Value(e) > 0.0 ? ArgA.Value(e) : ArgB.Value(e);} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgC.Value(e) > 0.0 ? ArgA.Variable(e) : ArgB.Variable(e);} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new condition_block_variable(*this));} + }; + + /// This class represents multiplication of the matrix by the variable. + class multiply_variable_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression + stored_variable_expression ArgB; + public: + /// Constructor. + multiply_variable_block_variable(const abstract_dynamic_block_variable & _ArgA, + const abstract_dynamic_variable & _ArgB) + : ArgA(_ArgA), ArgB(_ArgB) {} + /// Copy constructor. + multiply_variable_block_variable(const multiply_variable_block_variable & b) + : ArgA(b.ArgA), ArgB(b.ArgB) {} + /// Assignment operator. + multiply_variable_block_variable & operator =(multiply_variable_block_variable const & b) + { + ArgB = b.ArgB; + ArgA = b.ArgA; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e)*ArgB.Value(e);} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e)*ArgB.Variable(e);} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new multiply_variable_block_variable(*this));} + }; + + /// This class represents inverse of the matrix, this is A^{-1}. + class inverse_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression + public: + /// Constructor. + inverse_block_variable(const abstract_dynamic_block_variable & _ArgA) + : ArgA(_ArgA) {} + /// Copy constructor. + inverse_block_variable(const inverse_block_variable & other) + : ArgA(other.ArgA) {} + /// Assignment operator. + inverse_block_variable & operator =(inverse_block_variable const & other) + { + ArgA = other.ArgA; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e).Invert(true).first;} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e).Invert(true).first;} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new inverse_block_variable(*this));} + }; + + /// This class represents pseudo-inverse of the matrix, this is A^{+}. + class pseudo_inverse_block_variable : public abstract_dynamic_block_variable + { + private: + stored_block_variable_expression ArgA; //< Block variable expression + double eps; + public: + /// Constructor. + pseudo_inverse_block_variable(const abstract_dynamic_block_variable & _ArgA, double _eps = 1.0e-13) + : ArgA(_ArgA), eps(_eps) {} + /// Copy constructor. + pseudo_inverse_block_variable(const pseudo_inverse_block_variable & other) + : ArgA(other.ArgA), eps(other.eps) {} + /// Assignment operator. + pseudo_inverse_block_variable & operator =(pseudo_inverse_block_variable const & other) + { + ArgA = other.ArgA; + eps = other.eps; + return *this; + } + /// Get value of variable expression on provided element e. + rMatrix Value(const Storage & e) const + {return ArgA.Value(e).PseudoInvert(eps,true).first;} + /// Get value with derivatives of variable expression on provided element e. + /// This function collapses associated expression tree into multivar_expression. + vMatrix Variable(const Storage & e) const + {return ArgA.Variable(e).PseudoInvert(eps,true).first;} + /// Make a copy of this class, used to reproduce and store a tree of variable expressions. + abstract_dynamic_block_variable * Copy() const {return static_cast(new pseudo_inverse_block_variable(*this));} + }; + + + + + typedef abstract_dynamic_block_variable abstract_block_variable; +} + +/// If control evaluates to non-negative number, returns A, otherwise B +__INLINE +INMOST::condition_block_variable +condition(INMOST::abstract_dynamic_variable const & control, + INMOST::abstract_dynamic_block_variable const & if_ge_zero, + INMOST::abstract_dynamic_block_variable const & if_lt_zero) +{ return INMOST::condition_block_variable(control,if_ge_zero,if_lt_zero); } +/// Matrix operation A+B +__INLINE +INMOST::addition_block_variable +operator+(INMOST::abstract_block_variable const & Left, + INMOST::abstract_block_variable const & Right) +{ return INMOST::addition_block_variable(Left, Right); } +/// Matrix operation A-B +__INLINE +INMOST::subtraction_block_variable +operator-(INMOST::abstract_block_variable const & Left, + INMOST::abstract_block_variable const & Right) +{ return INMOST::subtraction_block_variable(Left, Right); } +/// Matrix operation A*B +__INLINE +INMOST::multiplication_block_variable +operator*(INMOST::abstract_block_variable const & Left, + INMOST::abstract_block_variable const & Right) +{ return INMOST::multiplication_block_variable(Left, Right); } +/// Matrix operation B^{-1}*A +__INLINE +INMOST::division_block_variable +operator/(INMOST::abstract_block_variable const & Left, + INMOST::abstract_block_variable const & Right) +{ return INMOST::division_block_variable(Left, Right); } +/// Attach two matrices by rows, i.e. C = [A,B] +__INLINE +INMOST::concat_rows_block_variable +concat_rows(INMOST::abstract_block_variable const & Left, + INMOST::abstract_block_variable const & Right) +{ return INMOST::concat_rows_block_variable(Left, Right); } +/// Attach two matrices by columns, i.e. C = [A^T,B^T]^T +__INLINE +INMOST::concat_cols_block_variable +concat_cols(INMOST::abstract_block_variable const & Left, + INMOST::abstract_block_variable const & Right) +{ return INMOST::concat_cols_block_variable(Left, Right); } +/// Matrix operation A^T +__INLINE +INMOST::transpose_block_variable +transpose(INMOST::abstract_block_variable const & Left) +{ return INMOST::transpose_block_variable(Left); } +/// Matrix operation A^{-1} +__INLINE +INMOST::inverse_block_variable +inv(INMOST::abstract_block_variable const & Left) +{ return INMOST::inverse_block_variable(Left); } +/// Matrix operation A^{+} +__INLINE +INMOST::pseudo_inverse_block_variable +pinv(INMOST::abstract_block_variable const & Left) +{ return INMOST::pseudo_inverse_block_variable(Left); } +/// Submatrix of a matrix +__INLINE +INMOST::submatrix_block_variable +submatrix(INMOST::abstract_block_variable const & Left, + INMOST_DATA_ENUM_TYPE row1, + INMOST_DATA_ENUM_TYPE row2, + INMOST_DATA_ENUM_TYPE col1, + INMOST_DATA_ENUM_TYPE col2) +{ return INMOST::submatrix_block_variable(Left,row1,row2,col1,col2); } +/// Matrix multiplication by a constant A*a +__INLINE +INMOST::multiply_const_block_variable +operator *(INMOST::abstract_block_variable const & Left, + INMOST_DATA_REAL_TYPE Right) +{ return INMOST::multiply_const_block_variable(Left,Right); } +/// Matrix multiplication by a constant a*A +__INLINE +INMOST::multiply_const_block_variable +operator *(INMOST_DATA_REAL_TYPE Left, + INMOST::abstract_block_variable const & Right) +{ return INMOST::multiply_const_block_variable(Right,Left); } +/// Matrix multiplication by inverse of a constant A/a +__INLINE +INMOST::multiply_const_block_variable +operator /(INMOST::abstract_block_variable const & Left, + INMOST_DATA_REAL_TYPE Right) +{ return INMOST::multiply_const_block_variable(Left,1.0/Right); } +/// Matrix multiplication by a scalar expression A*a +__INLINE +INMOST::multiply_variable_block_variable +operator *(INMOST::abstract_block_variable const & Left, + INMOST::abstract_variable const & Right) +{ return INMOST::multiply_variable_block_variable(Left,Right); } +/// Matrix multiplication by a scalar expression a*A +__INLINE +INMOST::multiply_variable_block_variable +operator *(INMOST::abstract_variable const & Left, + INMOST::abstract_block_variable const & Right) +{ return INMOST::multiply_variable_block_variable(Right,Left); } +/// Matrix multiplication by inverse of a scalar expression A/a +__INLINE +INMOST::multiply_variable_block_variable +operator /(INMOST::abstract_block_variable const & Left, + INMOST::abstract_variable const & Right) +{ return INMOST::multiply_variable_block_variable(Left,1.0/INMOST::stored_variable_expression(Right)); } +/// Calculation of matrix convex combination on stencil +__INLINE +INMOST::stencil_block_variable +stencil(INMOST::Tag tag_elems, + INMOST::Tag tag_coefs, + INMOST::abstract_dynamic_block_variable const & Arg) +{ return INMOST::stencil_block_variable(tag_elems,tag_coefs,Arg); } +/// Operation for key-value table +__INLINE +INMOST::table_block_variable +get_table(INMOST::abstract_dynamic_block_variable const & Arg, + INMOST::keyval_table const & Table) +{return INMOST::table_block_variable(Arg,Table);} +/// Branching expression by element type +__INLINE +INMOST::etype_branch_block_variable +etype_branch(INMOST::ElementType true_type, + INMOST::abstract_dynamic_block_variable const & iftrue, + INMOST::abstract_dynamic_block_variable const & iffalse) +{return INMOST::etype_branch_block_variable(true_type,iftrue,iffalse);} +/// Branching expression by marker +__INLINE +INMOST::marker_branch_block_variable +marker_branch(INMOST::MarkerType marker, + INMOST::abstract_dynamic_block_variable const & iftrue, + INMOST::abstract_dynamic_block_variable const & iffalse) +{return INMOST::marker_branch_block_variable(marker,iftrue,iffalse);} + +#endif //defined(USE_AUTODIFF) && defined(USE_MESH) + + + + +#endif //INMOST_AUTODIFF_ETBVAR_H_INCLUDED + diff --git a/Source/Headers/inmost_common.h b/Source/Headers/inmost_common.h index 886d102b383d7b2d3781d5566997126cf581452c..c0656121e2600ab5d7ec95f2bcc0ad9b02405d86 100644 --- a/Source/Headers/inmost_common.h +++ b/Source/Headers/inmost_common.h @@ -54,7 +54,7 @@ // output xml files for debugging of parallel algorithms // search for style.xsl within examples for comfortable // view of generated xml files -//#define USE_PARALLEL_WRITE_TIME +#define USE_PARALLEL_WRITE_TIME // this will revert Mesh::PrepareReceiveInner to always // use MPI point to point functionality disregarding problem type @@ -100,9 +100,15 @@ #include #include #include + #if defined(USE_OMP) -#include -#endif +#define OMP_THREAD omp_get_thread_num() +#define MAX_THREADS omp_get_max_threads() +#else //USE_OMP +#define OMP_THREAD 0 +#define MAX_THREADS 1 +#endif //USE_OMP + #if defined(min) #undef min @@ -255,4 +261,20 @@ namespace INMOST //#include "io.hpp" +namespace INMOST +{ + template + class AbstractMatrix; + + template + class SubMatrix; + + template > + class Matrix; + + template > + class SymmetricMatrix; +} + + #endif //INMOST_COMMON_INCLUDED diff --git a/Source/Headers/inmost_data.h b/Source/Headers/inmost_data.h index adf8f90586764224e558fb55fa0f3aa9a9382550..844b2c923e9c261a9360f8b009172f22af1fe0c5 100644 --- a/Source/Headers/inmost_data.h +++ b/Source/Headers/inmost_data.h @@ -180,8 +180,8 @@ namespace INMOST friend class Storage; }; - - ///This class provides the access to the individual mesh datum and general information about it. + + ///This class provides the access to the individual mesh datum and general information about it. class Tag //implemented in tag.cpp { private: @@ -206,7 +206,13 @@ namespace INMOST __INLINE Tag & operator =(Tag const & other); __INLINE DataType GetDataType() const; __INLINE INMOST_MPI_Type GetBulkDataType() const; + /// Amount of bytes necessery to support one record + /// referred by the tag on one element. Used internally + /// to allocate, manage and copy data of the mesh. __INLINE INMOST_DATA_ENUM_TYPE GetBytesSize() const; + /// Amount of bytes necessery for one record in packed + /// form that is used in GetData. + INMOST_DATA_ENUM_TYPE GetPackedBytesSize() const; __INLINE INMOST_DATA_ENUM_TYPE GetSize() const; __INLINE std::string GetTagName() const; __INLINE bool isDefined(ElementType type) const; diff --git a/Source/Headers/inmost_dense.h b/Source/Headers/inmost_dense.h index da4a049e182f61bc52e760628f7386361e56b8f8..ca629a6f9f66130d3d881f849eddc423b78bb39b 100644 --- a/Source/Headers/inmost_dense.h +++ b/Source/Headers/inmost_dense.h @@ -27,6 +27,15 @@ namespace INMOST template<> struct Promote {typedef variable type;}; template<> struct Promote {typedef hessian_variable type;}; template<> struct Promote {typedef hessian_variable type;}; + //for unknown + //For INMOST_DATA_INTEGER_TYPE + template<> struct Promote {typedef variable type;}; + template<> struct Promote {typedef variable type;}; + template<> struct Promote {typedef variable type;}; + template<> struct Promote {typedef variable type;}; + template<> struct Promote {typedef variable type;}; + template<> struct Promote {typedef hessian_variable type;}; + template<> struct Promote {typedef hessian_variable type;}; //For multivar_expression_reference template<> struct Promote {typedef variable type;}; template<> struct Promote {typedef variable type;}; @@ -61,14 +70,7 @@ namespace INMOST template<> struct Promote {typedef hessian_variable type;}; #endif - template - class SubMatrix; - - template > - class Matrix; - template > - class SymmetricMatrix; /// Abstract class for a matrix, /// used to abstract away all the data storage and access @@ -1514,7 +1516,8 @@ namespace INMOST /// @return Returns a unit matrix. static Matrix Unit(enumerator pn, const Var & c = 1.0) { - Matrix ret(pn,pn,0.0); + Matrix ret(pn,pn); + ret.Zero(); for(enumerator i = 0; i < pn; ++i) ret(i,i) = c; return ret; } @@ -2362,13 +2365,30 @@ namespace INMOST typedef Matrix iMatrix; /// shortcut for matrix of real values. typedef Matrix rMatrix; + /// shortcut for matrix of integer values in existing array. + typedef Matrix > iaMatrix; + /// shortcut for matrix of real values in existing array. + typedef Matrix > raMatrix; + /// return a matrix + static iaMatrix iaMatrixMake(INMOST_DATA_INTEGER_TYPE * p, iaMatrix::enumerator n, iaMatrix::enumerator m) {return iaMatrix(shell(p,n*m),n,m);} + static raMatrix raMatrixMake(INMOST_DATA_REAL_TYPE * p, raMatrix::enumerator n, raMatrix::enumerator m) {return raMatrix(shell(p,n*m),n,m);} #if defined(USE_AUTODIFF) -/// shortcut for matrix of variables with single unit entry of first order derivative. + /// shortcut for matrix of variables with single unit entry of first order derivative. typedef Matrix uMatrix; /// shortcut for matrix of variables with first order derivatives. typedef Matrix vMatrix; //< shortcut for matrix of variables with first and second order derivatives. typedef Matrix hMatrix; + /// shortcut for matrix of unknowns in existing array. + typedef Matrix > uaMatrix; + /// shortcut for matrix of variables in existing array. + typedef Matrix > vaMatrix; + /// shortcut for matrix of variables in existing array. + typedef Matrix > haMatrix; + + static uaMatrix uaMatrixMake(unknown * p, uaMatrix::enumerator n, uaMatrix::enumerator m) {return uaMatrix(shell(p,n*m),n,m);} + static vaMatrix vaMatrixMake(variable * p, vaMatrix::enumerator n, vaMatrix::enumerator m) {return vaMatrix(shell(p,n*m),n,m);} + static haMatrix vaMatrixMake(hessian_variable * p, haMatrix::enumerator n, haMatrix::enumerator m) {return haMatrix(shell(p,n*m),n,m);} #endif } /// Multiplication of matrix by constant from left. @@ -2379,6 +2399,14 @@ template INMOST::Matrix::type> operator *(INMOST_DATA_REAL_TYPE coef, const INMOST::AbstractMatrix & other) {return other*coef;} #if defined(USE_AUTODIFF) +/// Multiplication of matrix by a unknown from left. +/// Takes account for derivatives of variable. +/// @param coef Variable coefficient multiplying matrix. +/// @param other Matrix to be multiplied. +/// @return Matrix, each entry multiplied by a variable. +template +INMOST::Matrix::type> operator *(const INMOST::unknown & coef, const INMOST::AbstractMatrix & other) +{return other*coef;} /// Multiplication of matrix by a variable from left. /// Takes account for derivatives of variable. /// @param coef Variable coefficient multiplying matrix. diff --git a/Source/Headers/inmost_expression.h b/Source/Headers/inmost_expression.h index e7dbeabc42fa2bff668dc15ae888f550ccb8a02e..a46f63c27ce444916579c3b22546d812699f16c2 100644 --- a/Source/Headers/inmost_expression.h +++ b/Source/Headers/inmost_expression.h @@ -120,8 +120,6 @@ namespace INMOST #if defined(PACK_ARRAY) #pragma pack(push,r1,4) #endif - - /// A class that represents a variable with multiple /// first order variations. /// Short type name is variable. @@ -554,12 +552,12 @@ namespace INMOST } friend class hessian_multivar_expression_reference; }; - #if defined(PACK_ARRAY) #pragma pack(pop,r1) #endif static INMOST_DATA_REAL_TYPE stub_multivar_expression_reference_value; //for default constructor in multivar_expression_reference + class multivar_expression_reference : public shell_expression { @@ -1878,47 +1876,6 @@ namespace INMOST } }; - template - class stencil_expression : public shell_expression > - { - dynarray< const_multiplication_expression, 64 > arg; - INMOST_DATA_REAL_TYPE value; - public: - stencil_expression(const dynarray< const_multiplication_expression, 64 > & parg) : arg(parg) - { - value = 0.0; - for(typename dynarray< const_multiplication_expression, 64 >::iterator it = arg.begin(); it != arg.end(); ++it) - value += it->GetValue(); - } - stencil_expression(const stencil_expression & other) : arg(other.arg), value(other.value) {} - __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; } - __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const - { - for(typename dynarray< const_multiplication_expression, 64 >::iterator it = arg.begin(); it != arg.end(); ++it) - it->GetJacobian(mult,r); - } - __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const - { - for(typename dynarray< const_multiplication_expression, 64 >::iterator it = arg.begin(); it != arg.end(); ++it) - it->GetJacobian(mult,r); - } - __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const - { - Sparse::Row tmpJ, curJ; - Sparse::HessianRow tmpH, curH; - for(typename dynarray< const_multiplication_expression, 64 >::iterator it = arg.begin(); it != arg.end(); ++it) - { - curJ.Clear(); - curH.Clear(); - it->GetHessian(multJ,curJ,multH,curH); - Sparse::Row::MergeSortedRows(1.0,curJ,1.0,J,tmpJ); - Sparse::HessianRow::MergeSortedRows(1.0,curH,1.0,H,tmpH); - J.Swap(tmpJ); - H.Swap(tmpH); - } - } - }; - template class function_expression : public shell_expression< function_expression > diff --git a/Source/Headers/inmost_model.h b/Source/Headers/inmost_model.h new file mode 100644 index 0000000000000000000000000000000000000000..0a132573be8d5f5ed7f5d2a42a4346eaa9ada332 --- /dev/null +++ b/Source/Headers/inmost_model.h @@ -0,0 +1,128 @@ +#ifndef INMOST_MODEL_INCLUDED +#define INMOST_MODEL_INCLUDED + +#include "inmost_common.h" +#include "inmost_autodiff.h" +#include "inmost_residual.h" +#include "inmost_operator.h" + +#if defined(USE_AUTODIFF) && defined(USE_MESH) && defined(USE_SOLVER) + +namespace INMOST +{ + + class Model; + class AbstractOperator; + + /// A class to manage a submodel within a model. + /// Each submodel is responsible to define unknowns + /// of the model on the mesh. + class AbstractSubModel + { + public: + /// Let the submodel introduce it's unknowns + virtual bool PrepareEntries(Model & P) = 0; + /// Initialize coupling and dependent unknowns. + virtual bool Initialize(Model & P) = 0; + /// Fill part of the residual related to my unknowns. + virtual bool FillResidual(Residual & R) const = 0; + /// Update solution. + virtual bool UpdateSolution(const Sparse::Vector & sol) = 0; + /// Update time step. + virtual bool UpdateTimeStep() = 0; + /// Provide time step. + virtual bool SetTimeStep(double dt) = 0; + /// Roll back to previous step. + virtual bool RestoreTimeStep() = 0; + }; + + /// A class to organize a model. + class Model + { + Automatizator & aut; ///< Automatizator that is used to manage all unknowns of the model. + //todo: decide later on how it should be stored + std::vector< std::pair > SubModels; ///< A set of submodels of the model. + std::vector< std::pair< std::string, AbstractEntry *> > Entries; ///< A set of entries for blocks of unknowns of the model. + std::vector< std::pair > Meshes; ///< A set of meshes of the model. + std::vector< std::pair > Operators; ///< A set of operators used by submodels + bool initialized; ///< Indicates whether a model was initialized. + public: + Model(Automatizator & aut) : aut(aut), initialized(false) {} + //todo: + //Model(const Model & b) aut(b.aut) {} + //todo: + //Model & operator =(Model const & b); + virtual ~Model() {} + /// Add an entry of block unknowns to a model. + /// The model stores a link to the entry and may modify it contents. + /// The intries should be added from Model::Initialize function, + /// either by model or by any of the submodels. + void AddEntry(std::string name, AbstractEntry & entry); + /// Add a mesh to a model. + /// The model stores a link to the provided mesh, so it should not + /// be deallocated. The meshes are provided by the user from outside + /// before Model::Initialize function was called. + /// The meshes are requested by name by each submodel. + /// Same mesh can be added with different names for submodels. + void AddMesh(std::string name, Mesh & m); + /// Add a submodel to a model. + /// Submodels are added by the user from outside. + /// All submodels are initialized on Model::Initialize function. + void AddSubModel(std::string name, AbstractSubModel & submodel); + /// Add an operator to a model. + void AddOperator(std::string name, AbstractOperator & op); + /// Retrive an automatizator. + Automatizator & GetAutomatizator() {return aut;} + /// Retrive an automatizator. + const Automatizator & GetAutomatizator() const {return aut;} + /// Retrive a mesh by name. + Mesh * GetMesh(std::string); + /// Retrive a mesh by name. + const Mesh * GetMesh(std::string) const; + /// Retrive all names of meshes. + std::vector GetMeshesNames() const; + /// Retrive an entry that describe unknowns of the model by name. + AbstractEntry * GetEntry(std::string name); + /// Retrive an entry that describe unknowns of the model by name. + const AbstractEntry * GetEntry(std::string name) const; + /// Retrive all names of entries. + std::vector GetEntriesNames() const; + /// Retrive a submodel of the model by name. + AbstractSubModel * GetSubModel(std::string name); + /// Retrive a submodel of the model by name. + const AbstractSubModel * GetSubModel(std::string name) const; + /// Retrive all names of submodules. + std::vector GetSubModelsNames() const; + /// Retrive an operator of the model by name. + AbstractOperator * GetOperator(std::string name); + /// Retrive an operator of the model by name. + const AbstractOperator * GetOperator(std::string name) const; + /// Retrive all names of operators. + std::vector GetOperatorsNames() const; + /// Each submodel introduces it's unknowns into the model + /// so that later it can be accessed + bool PrepareEntries(); + /// Initialze all entries and submodels. + bool Initialize(); + /// Compute the residual of the model. + bool FillResidual(Residual & R) const; + /// Update solution. + bool UpdateSolution(const Sparse::Vector & sol); + /// Move to the next time step + bool UpdateTimeStep(); + /// Provide new time step. + bool SetTimeStep(double dt); + /// Roll back to previous time step + bool RestoreTimeStep(); + /// Check was the model initialized. + bool isInitialized() const {return initialized;} + /// Update variables contained in all block of automatizator on ghost elements of the grid. + /// For synchronization of data in individual blocks see AbstractEntry::SynhronizeData. + void SynchronizeData() { aut.SynchronizeData(); } + }; +} + +#endif //USE_AUTODIFF && USE_MESH && USE_SOLVER + +#endif //INMOST_MODEL_INCLUDED + diff --git a/Source/Headers/inmost_nonlinear.h b/Source/Headers/inmost_nonlinear.h index 3ec727e65ca0e9bda5b012fb67e32b0eda7eb70c..e61d5e7b015e844e871604f7d962602f10107b8d 100644 --- a/Source/Headers/inmost_nonlinear.h +++ b/Source/Headers/inmost_nonlinear.h @@ -2,31 +2,35 @@ #define INMOST_NONLINEAR_INCLUDED #include "inmost_common.h" - +#include "inmost_autodiff.h" #if defined(USE_NONLINEAR) namespace INMOST { - - typedef INMOST_DATA_BULK_TYPE RequestedAction; - - static const RequestedAction COMPUTE_FUNCTION = 0x01; - static const RequestedAction COMPUTE_JACOBIAN = 0x02; - static const RequestedAction COMPUTE_HESSIAN = 0x04; - static const RequestedAction FINISHED = 0x08; - - - class NonlinearSolver - { - public: - RequestedAction GetAction() const; - INMOST_DATA_REAL_TYPE GetResidual() const; - INMOST_DATA_ENUM_TYPE GetIterations() const; - std::string GetReason() const; - - }; - + + typedef INMOST_DATA_BULK_TYPE RequestedAction; + + static const RequestedAction COMPUTE_FUNCTION = 0x01; + static const RequestedAction COMPUTE_JACOBIAN = 0x02; + static const RequestedAction COMPUTE_HESSIAN = 0x04; + static const RequestedAction FINISHED = 0x08; + + class NonlinearSolver + { + Automatizator & aut; + public: + NonlinearSolver(Automatizator & aut) : aut(aut) {} + NonlinearSolver(const NonlinearSolver & b) : aut(b.aut) {} + NonlinearSolver & operator =(NonlinearSolver const & b) {aut = b.aut; return *this;} + ~NonlinearSolver() {} + + RequestedAction GetAction() const; + INMOST_DATA_REAL_TYPE GetResidual() const; + INMOST_DATA_ENUM_TYPE GetIterations() const; + std::string GetReason() const; + }; + }; #endif -#endif //INMOST_NONLINEAR_INCLUDED \ No newline at end of file +#endif //INMOST_NONLINEAR_INCLUDED diff --git a/Source/Headers/inmost_operator.h b/Source/Headers/inmost_operator.h new file mode 100644 index 0000000000000000000000000000000000000000..67a5a31cce132d5f9bee072a0994b88dcf4e57d1 --- /dev/null +++ b/Source/Headers/inmost_operator.h @@ -0,0 +1,116 @@ +#ifndef INMOST_OPERATOR_INCLUDED +#define INMOST_OPERATOR_INCLUDED + +#include "inmost_variable.h" + +#if defined(USE_AUTODIFF) && defined(USE_MESH) + +namespace INMOST +{ + class Model; + /// This class is responsible to unite access to various point-wise + /// implementations of discrete operators, such as grad, curl. + /// \todo + /// 1. Different types of operators: + /// time-stepping, + /// local point-wise (curl,grad on element), + /// global integrators (div,curl on domain), + /// interpolators, + /// inter-mesh interpolators. + /// Each has it's own functions. Implementation should be + /// flexible enough to prevent limitation. + /// 2. Ultimately operators should stack together: + /// for staggered incompressible navier-stokes: + /// Time(nU) + Projection(Divergence(ConvectionDiffusion(nU,\mu,Reconstruction(nU)))) - Grad(P) = f + /// Divergence(nU) = 0 + /// + class AbstractOperator + { + public: + /// Destroy all the data of the operator. + virtual ~AbstractOperator() {}; + /// Operator may introduce it's unknowns on this step + virtual bool PrepareEntries(Model & m) = 0; + /// Initialize all the data necessery to evalute the operator. + virtual bool Initialize(Model & m) = 0; + /// Let operator prepare data on the mesh before evaluation. + virtual bool Prepare() = 0; + /// This allows operator to update introduced variables. + /// May be useful for mimetic finite differences operators. + virtual bool UpdateSolution(const Sparse::Vector & V) = 0; + /// This allows operator to update time-dependent variables. + virtual bool UpdateTimeStep() = 0; + /// Check, whether we need to compute operator on this element. + virtual bool isValid(const Storage & e) const = 0; + /// Provide data necessery for initialization of the operator. + virtual bool SetTag(std::string name, Tag t) = 0; + /// Provide data necessery for initialization of the operator. + virtual bool SetMarker(std::string name, MarkerType marker) = 0; + }; + + class AbstractSpaceOperator : virtual public AbstractOperator + { + public: + virtual ~AbstractSpaceOperator() {} + /// Rebuild the discretization of the operator due to update of boundary conditions. + //virtual bool UpdateBoundaryCondition() = 0; + }; + + struct EquationPosition + { + const AbstractEntry & entry; + INMOST_DATA_ENUM_TYPE pos; + EquationPosition(const AbstractEntry & entry, INMOST_DATA_ENUM_TYPE pos) : entry(entry), pos(pos) {} + INMOST_DATA_ENUM_TYPE GetPosition(const Storage & e) const {return entry.Index(e,pos);} + }; + + class AbstractScalarOperator : virtual public AbstractSpaceOperator + { + public: + virtual ~AbstractScalarOperator() {} + /// The main function provided by operator. + /// @param R Residual that is filled in. + /// @param e Element on which operator is evaluated. + /// @param param Expression to be evaluated by operator. + /// @param eqpos Equation position in residual that operator affects. + /// @return Success or failure of operator application on provied expression. + virtual bool Evaluate(Residual & R, const Storage & e, const abstract_variable & param, const EquationPosition & eqpos) const = 0; + }; + + class AbstractVectorOperator : virtual public AbstractSpaceOperator + { + public: + virtual ~AbstractVectorOperator() {} + /// The main function provided by operator. + /// @param expr Operator may have multiple expression types, each is accessable by number. + /// @param param Expression to be evaluated by operator. + /// @return Result of operator application on provied expression in the form of 3x1 matrix. + virtual void Evaluate(Residual & R, const Storage & e, const stored_variable_expression param[3], const EquationPosition eqpos[3]) const = 0; + }; + + class AbstractTensorOperator : virtual public AbstractSpaceOperator + { + public: + virtual ~AbstractTensorOperator() {} + /// The main function provided by operator. + /// @param expr Operator may have multiple expression types, each is accessable by number. + /// @param param Expression to be evaluated by operator. + /// @return Result of operator application on provied expression in the form of 3x3 matrix. + virtual bool Evaluate(Residual & R, const Storage & e, const stored_variable_expression param[9], const EquationPosition eqpos[9]) const = 0; + }; + + class AbstractVoigtOperator : virtual public AbstractSpaceOperator + { + public: + virtual ~AbstractVoigtOperator() {} + /// The main function provided by operator. + /// @param expr Operator may have multiple expression types, each is accessable by number. + /// @param param Expression to be evaluated by operator. + /// @return Result of operator application on provied expression in the form of 6x1 matrix. + virtual bool Evaluate(Residual & R, const Storage & e, const stored_variable_expression param[6], const EquationPosition eqpos[6]) const = 0; + }; +} + +#endif //USE_AUTODIFF && USE_MESH + +#endif //INMOST_OPERATOR_INCLUDED diff --git a/Source/Headers/inmost_residual.h b/Source/Headers/inmost_residual.h new file mode 100644 index 0000000000000000000000000000000000000000..5ef0e90408af77c8180de87c10f28c33485bf494 --- /dev/null +++ b/Source/Headers/inmost_residual.h @@ -0,0 +1,99 @@ +#ifndef INMOST_RESIDUAL_INCLUDED +#define INMOST_RESIDUAL_INCLUDED + +#include "inmost_common.h" +#include "inmost_sparse.h" + +#if defined(USE_AUTODIFF) && defined(USE_SOLVER) + +namespace INMOST +{ + /// The Residual class provides a representation for array of residuals of nonlinear equations. + /// By working with the residual class you automatically assemble right hand side and + /// the jacobian of a nonlinear system of equation. + /// Jacobian matrix has a sparse representation. + /// \todo + /// 1. Extend for hessian calculation. + class Residual + { + Sparse::HessianMatrix hessian; ///< Hessian matrix + Sparse::Matrix jacobian; ///< Jacobian matrix. + Sparse::Vector residual; ///< Right hand side vector. + Sparse::LockService locks; ///< Array of locks for openmp shared access. + public: + /// Constructor + /// @param name Name for the matrix and right hand side. Can be used to set options for linear solver. + /// @param start First index of the equation in the local partition. Use Automatizator::GetFirstIndex. + /// @param end Last index of the equation in the local partition. Use Automatizator::GetLastIndex. + /// @param _comm MPI Communicator. + Residual(std::string name = "", INMOST_DATA_ENUM_TYPE start = 0, INMOST_DATA_ENUM_TYPE end = 0, INMOST_MPI_Comm _comm = INMOST_MPI_COMM_WORLD); + /// Copy constructor. + /// \warning May be expensive if matrices are large. + Residual(const Residual & other); + /// Assignment operator. + /// \warning May be expensive if matrices are large. + Residual & operator =(Residual const & other); + /// Retrive the first index of the equations in the local partition. + INMOST_DATA_ENUM_TYPE GetFirstIndex() const {return residual.GetFirstIndex();} + /// Retrive the last index of the equations in the local partition. + INMOST_DATA_ENUM_TYPE GetLastIndex() const {return residual.GetLastIndex();} + /// Retrive the first and the last indices of the equations in the local partition + /// @param start The first index of the equations will be recorded here. + /// @param end The last index of the equations will be recorded here. + void GetInterval(INMOST_DATA_ENUM_TYPE & start, INMOST_DATA_ENUM_TYPE & end) const; + /// Assign the new first and last indices of the equations in the local partition. + /// @param start The new first index of the equations. + /// @param end The new last index of the equations. + void SetInterval(INMOST_DATA_ENUM_TYPE beg, INMOST_DATA_ENUM_TYPE end); + /// Retrive a residual value and a jacobian row corresponding to certain equation. + /// @param row Equation number. + /// @return A structure that can be used in or assigned an automatic differentiation expression. + __INLINE multivar_expression_reference operator [](INMOST_DATA_ENUM_TYPE row) + {return multivar_expression_reference(residual[row],&jacobian[row]);} + /// Retrive a vector of entries in residual, corresponding to a set of equations. + /// @param rows A row-vector of equation numbers. + /// @param A structure that can be used in or assigned an automatic differentiation matrix expression. + Matrix operator [](const AbstractMatrix & rows); + /// Retrive hessian matrix. Use in nonlinear solver. + Sparse::HessianMatrix & GetHessian() {return hessian;} + /// Retrive hessian matrix without right of modificaiton. + const Sparse::HessianMatrix & GetHessian() const {return hessian;} + /// Retrive jacobian matrix. Use in Sparse::Solver::Solve function. + Sparse::Matrix & GetJacobian() {return jacobian;} + /// Retrive jacobian matrix without right of modificaiton. + const Sparse::Matrix & GetJacobian() const {return jacobian;} + /// Retrive right hand side vector. Use in Sparse::Solver::Solve function. + Sparse::Vector & GetResidual() {return residual;} + /// Retrive right hand side vector without right of modification. + const Sparse::Vector & GetResidual() const {return residual;} + /// Zero out right hand side vector. + void ClearResidual(); + /// Remove all entries in jacobian matrix. + void ClearJacobian(); + /// Remove all entries in hessian matrix. + void ClearHessian(); + /// Zero out right hand side vector and remove all entries in jacobian matrix. + void Clear(); + /// Compute the second norm of the right hand side vector over all of the processors. + INMOST_DATA_REAL_TYPE Norm(); + /// Normalize jacobian rows to unit p-norms and scale right hand side accordingly. + /// Use ENUMUNDEF as p to scale to infinite-norm. + void Rescale(INMOST_DATA_ENUM_TYPE p = 2); + /// Initialize openmp locks. + void InitLocks() {locks.SetInterval(GetFirstIndex(),GetLastIndex());} + /// Lock an equation to avoid simultaneous shared access. + /// @param pos Equation number. + __INLINE void Lock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.Lock(pos);} + /// UnLock an equation to allow simultaneous shared access. + /// @param pos Equation number. + __INLINE void UnLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) locks.UnLock(pos);} + /// Try to lock the equation. + /// @param pos Equation number. + /// @return True if equation was locked. + __INLINE bool TestLock(INMOST_DATA_ENUM_TYPE pos) {if(!locks.Empty()) return locks.TestLock(pos); return false;} + }; +} //namespace INMOST + +#endif //USE_SOLVER && USE_AUTODIFF + +#endif //INMOST_RESIDUAL_INCLUDED diff --git a/Source/Headers/inmost_sparse.h b/Source/Headers/inmost_sparse.h index 3ebbc5a0def289673372e60e272cb5307d584130..bded045eb4ffc0c19e9a2cc652d9bd6d63931840 100644 --- a/Source/Headers/inmost_sparse.h +++ b/Source/Headers/inmost_sparse.h @@ -6,8 +6,6 @@ - - namespace INMOST { namespace Sparse @@ -103,6 +101,8 @@ namespace INMOST INMOST_DATA_REAL_TYPE & operator [](INMOST_DATA_ENUM_TYPE i) {return data[i];} /// Return i-th element of the vector. INMOST_DATA_REAL_TYPE operator [](INMOST_DATA_ENUM_TYPE i) const {return data[i];} + /// Return a block of elements. + INMOST::Matrix operator [](const INMOST::AbstractMatrix & rows) const; /// Return the global size of the vector. INMOST_DATA_ENUM_TYPE Size() const { return static_cast(data.size()); } /// Iterator pointing to the first value of the vector. @@ -643,6 +643,7 @@ namespace INMOST interval< INMOST_DATA_ENUM_TYPE, Row::entry > LinkedList; ///< Storage for linked list. std::vector< INMOST_DATA_ENUM_TYPE > NonlocalPre; ///< List of global indices, that are to the left of owned index interval std::vector< INMOST_DATA_ENUM_TYPE > NonlocalPost; ///< List of global indices, that are to the right of owned index interval + std::map< INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE > Nonlocal; ///< Additional space public: /// This function converts global index into local index. /// @param pos Global index. diff --git a/Source/Headers/inmost_variable.h b/Source/Headers/inmost_variable.h index 331fbe360674eb65d96bc60a756fc10e64410556..32ce35fa08d0175665692426818b87bb77ff706b 100644 --- a/Source/Headers/inmost_variable.h +++ b/Source/Headers/inmost_variable.h @@ -381,12 +381,58 @@ namespace INMOST abstract_dynamic_variable * Copy() const {return static_cast(new stored_variable(*this));} }; + template + class stencil_expression : public shell_expression > + { + public: + typedef const_multiplication_expression argument; + typedef unary_pool_expression< argument, A> pool; + typedef dynarray< argument, 64 > container; + private: + container arg; + INMOST_DATA_REAL_TYPE value; + public: + stencil_expression(const container & parg) : arg(parg) + { + value = 0.0; + for(typename container::iterator it = arg.begin(); it != arg.end(); ++it) + value += it->GetValue(); + } + stencil_expression(const stencil_expression & other) : arg(other.arg), value(other.value) {} + __INLINE INMOST_DATA_REAL_TYPE GetValue() const { return value; } + __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::RowMerger & r) const + { + for(typename container::iterator it = arg.begin(); it != arg.end(); ++it) + it->GetJacobian(mult,r); + } + __INLINE void GetJacobian(INMOST_DATA_REAL_TYPE mult, Sparse::Row & r) const + { + for(typename container::iterator it = arg.begin(); it != arg.end(); ++it) + it->GetJacobian(mult,r); + } + __INLINE void GetHessian(INMOST_DATA_REAL_TYPE multJ, Sparse::Row & J, INMOST_DATA_REAL_TYPE multH, Sparse::HessianRow & H) const + { + Sparse::Row tmpJ, curJ; + Sparse::HessianRow tmpH, curH; + for(typename container::iterator it = arg.begin(); it != arg.end(); ++it) + { + curJ.Clear(); + curH.Clear(); + it->GetHessian(multJ,curJ,multH,curH); + Sparse::Row::MergeSortedRows(1.0,curJ,1.0,J,tmpJ); + Sparse::HessianRow::MergeSortedRows(1.0,curH,1.0,H,tmpH); + J.Swap(tmpJ); + H.Swap(tmpH); + } + } + }; + template class stencil_variable : public shell_dynamic_variable< stencil_expression, stencil_variable > { private: - Tag tag_elems; - Tag tag_coefs; + TagReferenceArray tag_elems; + TagRealArray tag_coefs; A Arg; public: stencil_variable(Tag tag_elems, Tag tag_coefs, const shell_dynamic_variable & parg) : tag_elems(tag_elems), tag_coefs(tag_coefs), Arg(parg) {} @@ -400,13 +446,17 @@ namespace INMOST } stencil_expression operator [](const Storage & e) const { - dynarray< const_multiplication_expression, 64> tmp; - Storage::real_array coefs = e.RealArray(tag_coefs); - Storage::reference_array elems = e.RealArray(tag_elems); + typename stencil_expression::container tmp; + Storage::real_array coefs = tag_coefs[e]; + Storage::reference_array elems = tag_elems[e]; assert(coefs.size() == elems.size()); tmp.resize(elems.size()); for(INMOST_DATA_ENUM_TYPE k = 0; k < elems.size(); ++k) - tmp[k] = const_multiplication_expression(Arg[elems[k]],coefs[k]); + { + typename stencil_expression::argument arg(Arg[elems[k]],coefs[k]); + typename stencil_expression::pool pool(arg); + tmp[k] = pool; + } return stencil_expression(tmp); } void GetVariation(const Storage & e, Sparse::Row & r) const { (*this)[e].GetJacobian(1.0,r); } @@ -423,6 +473,7 @@ namespace INMOST table_variable(const shell_dynamic_variable & parg, const keyval_table & ptable) : Arg(parg), Table(ptable) {} table_variable(const table_variable & other) : Arg(other.Arg), Table(other.Table) {} table_variable & operator = (table_variable const & other) {Arg = other.Arg; Table = other.Table; return * this;} + INMOST_DATA_REAL_TYPE Value(const Storage & e) const {return (*this)[e].GetValue();} multivar_expression Variable(const Storage & e) const { multivar_expression ret = (*this)[e]; diff --git a/Source/Mesh/mesh.cpp b/Source/Mesh/mesh.cpp index 2e5cbc839e31956d7847d18ff970c7c155f7c7c0..150594d0bbe54acbc23be20b3dc74ed1d9212daa 100644 --- a/Source/Mesh/mesh.cpp +++ b/Source/Mesh/mesh.cpp @@ -2149,17 +2149,17 @@ namespace INMOST (void)h; (void)tag; (void)expected; //due to __INLINE these variables considered by compilers as unreferenced } - void Mesh::ClearMarkerSpace(HandleType h) + void Mesh::ClearMarkerSpace(HandleType h) { Storage::bulk * marker_space = static_cast(MGetDenseLink(h,MarkersTag())); for(INMOST_DATA_ENUM_TYPE k = 0; k < MarkerFields; ++k) marker_space[k] = 0; } - void Mesh::GetMarkerSpace(HandleType h, Storage::bulk copy[MarkerFields]) const + void Mesh::GetMarkerSpace(HandleType h, Storage::bulk copy[MarkerFields]) const { const Storage::bulk * marker_space = static_cast(MGetDenseLink(h,MarkersTag())); for(INMOST_DATA_ENUM_TYPE k = 0; k < MarkerFields; ++k) copy[k] = marker_space[k]; } - void Mesh::SetMarkerSpace(HandleType h, Storage::bulk source[MarkerFields]) + void Mesh::SetMarkerSpace(HandleType h, Storage::bulk source[MarkerFields]) { Storage::bulk * marker_space = static_cast(MGetDenseLink(h,MarkersTag())); for(INMOST_DATA_ENUM_TYPE k = 0; k < MarkerFields; ++k) marker_space[k] = source[k]; @@ -2178,65 +2178,72 @@ namespace INMOST case DATA_INTEGER: return static_cast(static_cast(adata)->size()); case DATA_BULK: return static_cast(static_cast(adata)->size()); case DATA_REFERENCE:return static_cast(static_cast(adata)->size()); - case DATA_REMOTE_REFERENCE: - return static_cast(static_cast(adata)->size()); + case DATA_REMOTE_REFERENCE: + return static_cast(static_cast(adata)->size()); #if defined(USE_AUTODIFF) - case DATA_VARIABLE: return static_cast(static_cast(adata)->size()); + case DATA_VARIABLE: return static_cast(static_cast(adata)->size()); #endif } throw BadTag; } return tag.GetSize(); } - INMOST_DATA_ENUM_TYPE Mesh::GetDataCapacity(const INMOST_DATA_BULK_TYPE * adata, INMOST_DATA_ENUM_TYPE size, const Tag & tag) const - { - assert( tag.GetMeshLink() == this ); + INMOST_DATA_ENUM_TYPE Mesh::GetDataCapacity(const INMOST_DATA_BULK_TYPE * adata, INMOST_DATA_ENUM_TYPE size, const Tag & tag) const + { + assert( tag.GetMeshLink() == this ); #if defined(USE_AUTODIFF) - if( tag.GetDataType() == DATA_VARIABLE ) - { - INMOST_DATA_ENUM_TYPE ret = 0; - const Sparse::Row::entry * arr = static_cast(static_cast(adata)); - for(INMOST_DATA_ENUM_TYPE k = 0; k < size; ++k) - ret += variable::RetriveSize(arr+ret); - return ret*sizeof(Sparse::Row::entry); - } - else + if( tag.GetDataType() == DATA_VARIABLE ) + { + INMOST_DATA_ENUM_TYPE ret = 0; + const Sparse::Row::entry * arr = static_cast(static_cast(adata)); + for(INMOST_DATA_ENUM_TYPE k = 0; k < size; ++k) + ret += variable::RetriveSize(arr+ret); + return ret*sizeof(Sparse::Row::entry); + } + else #endif - return size*tag.GetBytesSize(); - assert(false); - return 0; - } - INMOST_DATA_ENUM_TYPE Mesh::GetDataCapacity(HandleType h,const Tag & tag) const + return size*tag.GetBytesSize(); + assert(false); + return 0; + } + INMOST_DATA_ENUM_TYPE Mesh::GetDataCapacity(HandleType h,const Tag & tag) const { assert( tag.GetMeshLink() == this ); - if( tag.GetSize() == ENUMUNDEF ) + INMOST_DATA_ENUM_TYPE s = tag.GetSize(); + if( s == ENUMUNDEF ) { const void * adata = MGetLink(h,tag); assert( adata != NULL ); switch(tag.GetDataType()) { - case DATA_REAL: return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); + case DATA_REAL: return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); case DATA_INTEGER: return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); case DATA_BULK: return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); case DATA_REFERENCE:return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); - case DATA_REMOTE_REFERENCE: - return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); + case DATA_REMOTE_REFERENCE: + return static_cast(static_cast(adata)->size())*tag.GetBytesSize(); #if defined(USE_AUTODIFF) - case DATA_VARIABLE: - { - INMOST_DATA_ENUM_TYPE ret = 0; - const inner_variable_array * arr = static_cast(adata); - for(inner_variable_array::size_type k = 0; k < arr->size(); ++k) - ret += (*arr)[k].RecordSize();//(*arr)[k].GetRow().Size(); - return ret*sizeof(Sparse::Row::entry_s); - } + case DATA_VARIABLE: + { + INMOST_DATA_ENUM_TYPE ret = 0; + const inner_variable_array * arr = static_cast(adata); + for(inner_variable_array::size_type k = 0; k < arr->size(); ++k) + ret += (*arr)[k].RecordSize();//(*arr)[k].GetRow().Size(); + return ret*sizeof(Sparse::Row::entry_s); + } #endif } throw BadTag; } #if defined(USE_AUTODIFF) - if( tag.GetDataType() == DATA_VARIABLE ) - return static_cast(MGetLink(h,tag))->RecordSize()*sizeof(Sparse::Row::entry_s); + if( tag.GetDataType() == DATA_VARIABLE ) + { + INMOST_DATA_ENUM_TYPE ret = 0; + const var * v = static_cast(MGetLink(h,tag)); + for(INMOST_DATA_ENUM_TYPE r = 0; r < s; ++r) + ret += v[r].RecordSize()*sizeof(Sparse::Row::entry_s); + return ret; + } #endif return tag.GetSize()*tag.GetBytesSize(); } @@ -2253,10 +2260,10 @@ namespace INMOST case DATA_INTEGER: static_cast(adata)->resize(new_size); break; case DATA_BULK: static_cast(adata)->resize(new_size); break; case DATA_REFERENCE:static_cast(adata)->resize(new_size); break; - case DATA_REMOTE_REFERENCE: + case DATA_REMOTE_REFERENCE: static_cast(adata)->resize(new_size); break; #if defined(USE_AUTODIFF) - case DATA_VARIABLE: static_cast(adata)->resize(new_size); break; + case DATA_VARIABLE: static_cast(adata)->resize(new_size); break; #endif } return; @@ -2280,30 +2287,30 @@ namespace INMOST case DATA_INTEGER: memcpy(data_out,&(*static_cast(adata))[shift],bytes*size); break; case DATA_BULK: memcpy(data_out,&(*static_cast(adata))[shift],bytes*size); break; case DATA_REFERENCE:memcpy(data_out,&(*static_cast(adata))[shift],bytes*size); break; - case DATA_REMOTE_REFERENCE: + case DATA_REMOTE_REFERENCE: memcpy(data_out,&(*static_cast(adata))[shift],bytes*size); break; #if defined(USE_AUTODIFF) - case DATA_VARIABLE: - { - const inner_variable_array * arr = static_cast(adata); - Sparse::Row::entry_s * data = static_cast(data_out); - int k = 0; - for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) - k += (*arr)[r+shift].Record(data+k); - } - break; + case DATA_VARIABLE: + { + const inner_variable_array * arr = static_cast(adata); + Sparse::Row::entry_s * data = static_cast(data_out); + int k = 0; + for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) + k += (*arr)[r+shift].Record(data+k); + } + break; #endif } } #if defined(USE_AUTODIFF) - else if( tag.GetDataType() == DATA_VARIABLE ) - { - Sparse::Row::entry_s * data = static_cast(data_out); - const var * v = static_cast(MGetLink(h,tag)); - int k = 0; - for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) - k += v[r+shift].Record(data+k); - } + else if( tag.GetDataType() == DATA_VARIABLE ) + { + Sparse::Row::entry_s * data = static_cast(data_out); + const var * v = static_cast(MGetLink(h,tag)); + int k = 0; + for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) + k += v[r+shift].Record(data+k); + } #endif else memcpy(data_out,static_cast(adata)+shift*bytes,size*bytes); return; @@ -2323,30 +2330,30 @@ namespace INMOST case DATA_INTEGER: memcpy(&(*static_cast(adata))[shift],data_in,bytes*size); break; case DATA_BULK: memcpy(&(*static_cast(adata))[shift],data_in,bytes*size); break; case DATA_REFERENCE:memcpy(&(*static_cast(adata))[shift],data_in,bytes*size); break; - case DATA_REMOTE_REFERENCE: - memcpy(&(*static_cast(adata))[shift],data_in,bytes*size); break; + case DATA_REMOTE_REFERENCE: + memcpy(&(*static_cast(adata))[shift],data_in,bytes*size); break; #if defined(USE_AUTODIFF) - case DATA_VARIABLE: - { - inner_variable_array * arr = static_cast(adata); - const Sparse::Row::entry_s * data = static_cast(data_in); - int k = 0; - for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) - k += (*arr)[r+shift].Retrive(data+k); - } - break; + case DATA_VARIABLE: + { + inner_variable_array * arr = static_cast(adata); + const Sparse::Row::entry_s * data = static_cast(data_in); + int k = 0; + for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) + k += (*arr)[r+shift].Retrive(data+k); + } + break; #endif } } #if defined(USE_AUTODIFF) - else if( tag.GetDataType() == DATA_VARIABLE ) - { - const Sparse::Row::entry_s * data = static_cast(data_in); - var * v = static_cast(MGetLink(h,tag)); - int k = 0; - for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) - k += v[r+shift].Retrive(data+k); - } + else if( tag.GetDataType() == DATA_VARIABLE ) + { + const Sparse::Row::entry_s * data = static_cast(data_in); + var * v = static_cast(MGetLink(h,tag)); + int k = 0; + for(INMOST_DATA_ENUM_TYPE r = 0; r < size; ++r) + k += v[r+shift].Retrive(data+k); + } #endif else memcpy(static_cast(adata)+shift*bytes,data_in,size*bytes); } diff --git a/Source/Mesh/parallel.cpp b/Source/Mesh/parallel.cpp index de51a2735796f63e5cd5d329f410d11f2826141c..3761f839eff58096576cdc3516441cca51e84fff 100644 --- a/Source/Mesh/parallel.cpp +++ b/Source/Mesh/parallel.cpp @@ -1341,10 +1341,12 @@ namespace INMOST } if( !mapping.empty() ) std::sort(mapping.begin(),mapping.end(),MappingComparator()); + /* #if defined(USE_PARALLEL_WRITE_TIME) for(std::vector >::iterator it = mapping.begin(); it != mapping.end(); ++it) REPORT_STR("global " << it->first << " local " << it->second << Element::StatusName(ElementByLocalID(PrevElementType(current_mask),it->second)->GetStatus())); #endif + */ time = Timer() - time; REPORT_VAL("mapping size",mapping.size()) REPORT_STR("Compute global to local indexes mapping"); @@ -2243,9 +2245,10 @@ namespace INMOST void Mesh::PackTagData(const Tag & tag, const elements_by_type & elements, ElementType mask, MarkerType select, buffer_type & buffer) { - if( tag.GetDataType() == DATA_REFERENCE || tag.GetDataType() == DATA_REMOTE_REFERENCE ) return; //NOT IMPLEMENTED TODO 14 + if( tag.GetDataType() == DATA_REFERENCE || tag.GetDataType() == DATA_REMOTE_REFERENCE ) return; //NOT IMPLEMENTED TODO 14 ENTER_FUNC(); #if defined(USE_MPI) + REPORT_VAL("Buffer size before pack",buffer.size()); ElementType pack_types[2] = {NONE,NONE}; element_set::const_iterator eit; buffer_type array_data_send; @@ -2276,6 +2279,12 @@ namespace INMOST //array_data_send.resize(had_s+s*tag.GetBytesSize()); if( s ) { + if( tag.GetDataType() == DATA_VARIABLE ) + { + REPORT_VAL("data size: ", s); + REPORT_VAL("data capacity: ", GetDataCapacity(*eit,tag)); + REPORT_VAL("array size: ", had_s); + } array_data_send.resize(had_s+GetDataCapacity(*eit,tag)); GetData(*eit,tag,0,s,&array_data_send[had_s]); //REPORT_VAL("size",s); @@ -2298,6 +2307,13 @@ namespace INMOST //array_data_send.resize(had_s+s*tag.GetBytesSize()); if( s ) { + if( tag.GetDataType() == DATA_VARIABLE ) + { + REPORT_VAL("on element ",Element(this,*eit).GlobalID()); + REPORT_VAL("position: ", had_s); + REPORT_VAL("data capacity: ", GetDataCapacity(*eit,tag)); + REPORT_VAL("size: ", s); + } array_data_send.resize(had_s+GetDataCapacity(*eit,tag)); GetData(*eit,tag,0,s,&array_data_send[had_s]); } @@ -2313,15 +2329,17 @@ namespace INMOST REPORT_VAL("tag sparse on",static_cast(pack_types[1])); REPORT_VAL("size_size",array_size_send[0]); REPORT_VAL("data_size",array_size_send[1]); - int buffer_size = 0,position = static_cast(buffer.size()),temp; + int buffer_size = 0,position = static_cast(buffer.size()),temp,bytes; + bytes = tag.GetPackedBytesSize(); MPI_Pack_size(2 ,INMOST_MPI_DATA_BULK_TYPE,comm,&temp); buffer_size+= temp; MPI_Pack_size(static_cast(array_size_send.size()) ,INMOST_MPI_DATA_ENUM_TYPE,comm,&temp); buffer_size+= temp; - MPI_Pack_size(static_cast(array_data_send.size()/tag.GetBytesSize()),tag.GetBulkDataType() ,comm,&temp); buffer_size+= temp; + MPI_Pack_size(static_cast(array_data_send.size()/bytes),tag.GetBulkDataType() ,comm,&temp); buffer_size+= temp; buffer.resize(position+buffer_size); MPI_Pack(pack_types,2,INMOST_MPI_DATA_BULK_TYPE,&buffer[0],static_cast(buffer.size()),&position,comm); if( !array_size_send.empty() ) MPI_Pack(&array_size_send[0],static_cast(array_size_send.size()),INMOST_MPI_DATA_ENUM_TYPE,&buffer[0],static_cast(buffer.size()),&position,comm); - if( !array_data_send.empty() ) MPI_Pack(&array_data_send[0],static_cast(array_data_send.size()/tag.GetBytesSize()),tag.GetBulkDataType(),&buffer[0],static_cast(buffer.size()),&position,comm); + if( !array_data_send.empty() ) MPI_Pack(&array_data_send[0],static_cast(array_data_send.size()/bytes),tag.GetBulkDataType(),&buffer[0],static_cast(buffer.size()),&position,comm); buffer.resize(position); + REPORT_VAL("Buffer size after pack",buffer.size()); #else (void) tag; (void) elements; @@ -2341,6 +2359,7 @@ namespace INMOST REPORT_VAL("TagName",tag.GetTagName()); REPORT_VAL("select marker",select); #if defined(USE_MPI) + REPORT_VAL("Position before unpack",position); if( !buffer.empty() ) { int pos = 0, k = 0; @@ -2359,7 +2378,18 @@ namespace INMOST array_size_recv.resize(size_recv); array_data_recv.resize(data_recv); if( !array_size_recv.empty() ) MPI_Unpack(&buffer[0],static_cast(buffer.size()),&position,&array_size_recv[0],static_cast(array_size_recv.size()),INMOST_MPI_DATA_ENUM_TYPE,comm); - if( !array_data_recv.empty() ) MPI_Unpack(&buffer[0],static_cast(buffer.size()),&position,&array_data_recv[0],static_cast(array_data_recv.size()/tag.GetBytesSize()),tag.GetBulkDataType(),comm); + if( !array_data_recv.empty() ) + { + int bytes = tag.GetPackedBytesSize(); + REPORT_VAL("occupied by type",bytes); + REPORT_VAL("bytes in entry",sizeof(Sparse::Row::entry)); + REPORT_VAL("stored in type",tag.GetBytesSize()); + REPORT_VAL("all size recv",array_data_recv.size()); + REPORT_VAL("incoming size of data",data_recv); + REPORT_VAL("calculated size of data",array_data_recv.size()/bytes); + REPORT_VAL("calculated size of data",array_data_recv.size()/sizeof(Sparse::Row::entry)); + MPI_Unpack(&buffer[0],static_cast(buffer.size()),&position,&array_data_recv[0],static_cast(array_data_recv.size()/bytes),tag.GetBulkDataType(),comm); + } for(int i = ElementNum(NODE); i <= ElementNum(CELL); i++) if( (recv_mask[0] & ElementTypeFromDim(i)) ) { REPORT_VAL("unpack for type",ElementTypeName(ElementTypeFromDim(i))); @@ -2434,6 +2464,13 @@ namespace INMOST { if( !select || GetMarker(*eit,select) ) { + if( tag.GetDataType() == DATA_VARIABLE ) + { + REPORT_VAL("on element ",Element(this,*eit).GlobalID()); + REPORT_VAL("position ", pos); + REPORT_VAL("capacity ", GetDataCapacity(&array_data_recv[pos],size,tag)); + REPORT_VAL("size ", size); + } op(tag,Element(this,*eit),&array_data_recv[pos],size); pos += GetDataCapacity(&array_data_recv[pos],size,tag); //pos += size*tag.GetBytesSize(); @@ -2447,6 +2484,7 @@ namespace INMOST REPORT_VAL("total unpacked records",total_unpacked); } } + REPORT_VAL("Position after unpack",position); #else (void) tag; (void) elements; @@ -3832,13 +3870,17 @@ namespace INMOST MPI_Attr_get(comm,MPI_TAG_UB,&p_max_tag,&flag); #endif //USE_MPI2 if( flag ) max_tag = *p_max_tag; + REPORT_VAL("max_tag",max_tag); std::vector send_recv_size(send_bufs.size()+recv_bufs.size()); std::vector reqs(send_bufs.size()+recv_bufs.size()); - for(i = 0; i < send_bufs.size(); i++) send_recv_size[i+recv_bufs.size()] = static_cast(send_bufs[i].second.size()); + for(i = 0; i < send_bufs.size(); i++) + send_recv_size[i+recv_bufs.size()] = static_cast(send_bufs[i].second.size()); REPORT_VAL("recv buffers size",recv_bufs.size()); for(i = 0; i < recv_bufs.size(); i++) { mpi_tag = ((parallel_mesh_unique_id+1)*mpisize*mpisize + (mpirank+mpisize+rand_num))%max_tag; + REPORT_VAL("origin",recv_bufs[i].first); + REPORT_VAL("mpi_tag",mpi_tag); //mpi_tag = parallel_mesh_unique_id*mpisize*mpisize+recv_bufs[i].first*mpisize+mpirank; REPORT_MPI(MPI_Irecv(&send_recv_size[i],1,MPI_INT,recv_bufs[i].first,mpi_tag,comm,&reqs[i])); } @@ -3846,6 +3888,9 @@ namespace INMOST for(i = 0; i < send_bufs.size(); i++) { mpi_tag = ((parallel_mesh_unique_id+1)*mpisize*mpisize + (send_bufs[i].first+mpisize+rand_num))%max_tag; + REPORT_VAL("destination",send_bufs[i].first); + REPORT_VAL("mpi_tag",mpi_tag); + REPORT_VAL("size",send_recv_size[i+recv_bufs.size()]); //mpi_tag = parallel_mesh_unique_id*mpisize*mpisize+mpirank*mpisize+send_bufs[i].first; REPORT_MPI(MPI_Isend(&send_recv_size[i+recv_bufs.size()],1,MPI_INT,send_bufs[i].first,mpi_tag,comm,&reqs[i+recv_bufs.size()])); } @@ -3853,7 +3898,13 @@ namespace INMOST { REPORT_MPI(MPI_Waitall(static_cast(recv_bufs.size()),&reqs[0],MPI_STATUSES_IGNORE)); } - for(i = 0; i < recv_bufs.size(); i++) recv_bufs[i].second.resize(send_recv_size[i]); + REPORT_VAL("recieved buffers size",recv_bufs.size()); + for(i = 0; i < recv_bufs.size(); i++) + { + REPORT_VAL("origin",recv_bufs[i].first); + REPORT_VAL("size",send_recv_size[i]); + recv_bufs[i].second.resize(send_recv_size[i]); + } if( !send_bufs.empty() ) { REPORT_MPI(MPI_Waitall(static_cast(send_bufs.size()),&reqs[recv_bufs.size()],MPI_STATUSES_IGNORE)); diff --git a/Source/NonlinearSolver/nonlinear.cpp b/Source/NonlinearSolver/nonlinear.cpp index 1912ef4d736aae9adfbb11787c36d5ec0eb38527..17811180ecbaa8a0801a80cf793f085663f7e13f 100644 --- a/Source/NonlinearSolver/nonlinear.cpp +++ b/Source/NonlinearSolver/nonlinear.cpp @@ -1 +1 @@ -#include "inmost_nonlinear.h" \ No newline at end of file +#include "inmost_nonlinear.h" diff --git a/Source/Solver/SolverInterface.h b/Source/Solver/SolverInterface.h index aa0797670df2c3006739b5f56eb25e5158a499aa..d7ea68390301cccc733850d3217d081e06eb96ea 100644 --- a/Source/Solver/SolverInterface.h +++ b/Source/Solver/SolverInterface.h @@ -1,9 +1,10 @@ #ifndef INMOST_SOLVER_INTERFACE_H #define INMOST_SOLVER_INTERFACE_H +#include "inmost_sparse.h" #include #include -#include "inmost_sparse.h" + #define SILENCE_SET_PARAMETER #if defined(USE_SOLVER) diff --git a/Source/Solver/solver_inner/solver_ddpqiluc2/solver_ddpqiluc2.cpp b/Source/Solver/solver_inner/solver_ddpqiluc2/solver_ddpqiluc2.cpp index 27712bbd58931525acdacfe3fa35b6585b633e04..49e50dbab3b9d93172fe16f07d9e968fb312902b 100644 --- a/Source/Solver/solver_inner/solver_ddpqiluc2/solver_ddpqiluc2.cpp +++ b/Source/Solver/solver_inner/solver_ddpqiluc2/solver_ddpqiluc2.cpp @@ -1,9 +1,9 @@ #define _CRT_SECURE_NO_WARNINGS -#include + #include "inmost_solver.h" #if defined(USE_SOLVER) #include "solver_ddpqiluc2.hpp" - +#include /// \todo /// 1. Correct estimator (as in MPTILUC) diff --git a/Source/Solver/solver_inner/solver_ilu2/solver_ilu2.hpp b/Source/Solver/solver_inner/solver_ilu2/solver_ilu2.hpp index 7168ecb32cdb5f8c1f25b7e51613398964bed0cd..e223c64fced59afb4831372ff36ce13f4d3cd439 100644 --- a/Source/Solver/solver_inner/solver_ilu2/solver_ilu2.hpp +++ b/Source/Solver/solver_inner/solver_ilu2/solver_ilu2.hpp @@ -23,7 +23,7 @@ private: Sparse::Matrix *Alink; Solver::OrderInfo *info; //Sparse::Matrix L,U; - //Sparse::Vector div; + Sparse::Vector div; std::vector luv; std::vector lui; interval ilu, iu; @@ -444,13 +444,13 @@ public: //std::cout << __FUNCTION__ << " done" << std::endl; #endif - /* - info.PrepareVector(div); + /* + info->PrepareVector(div); std::fill(div.Begin(),div.End(),0); for(k = mobeg; k < moend; k++) div[k] = 1.0; - info.Accumulate(div); + info->Accumulate(div); for(k = mobeg; k < moend; k++) div[k] = 1.0/div[k]; - */ + */ init = true; return true; } @@ -505,11 +505,12 @@ public: bool Solve(Sparse::Vector &input, Sparse::Vector &output) { assert(isInitialized()); + INMOST_DATA_ENUM_TYPE mobeg, moend, k; #if defined(USE_OMP) #pragma omp single #endif { - INMOST_DATA_ENUM_TYPE mobeg, moend, r, k, vbeg, vend; //, end; + INMOST_DATA_ENUM_TYPE r, vbeg, vend; //, end; info->GetOverlapRegion(info->GetRank(), mobeg, moend); info->GetVectorRegion(vbeg, vend); for (k = vbeg; k < mobeg; k++) output[k] = 0; //Restrict additive schwartz (maybe do it outside?) @@ -528,9 +529,17 @@ public: output[k - 1] *= luv[iu[k - 1]]; } } - //May assemble partition of unity instead of restriction before accumulation - //assembly should be done instead of initialization - info->Accumulate(output); + //May assemble partition of unity instead of restriction before accumulation + //assembly should be done instead of initialization + info->Accumulate(output); + /* +#if defined(USE_OMP) +#pragma omp single +#endif + { + for (INMOST_DATA_ENUM_TYPE k = mobeg; k < moend; k++) output[k] *= div[k]; + } + */ return true; } diff --git a/Source/Solver/solver_inner/solver_mlmptiluc/solver_mlmtiluc2.cpp b/Source/Solver/solver_inner/solver_mlmptiluc/solver_mlmtiluc2.cpp index dfec800b290efa553143ff7ca74cc7eb6733c49c..0f92068b6e45a2f427f46a234bf05e7c68602fea 100644 --- a/Source/Solver/solver_inner/solver_mlmptiluc/solver_mlmtiluc2.cpp +++ b/Source/Solver/solver_inner/solver_mlmptiluc/solver_mlmtiluc2.cpp @@ -1,10 +1,10 @@ #define _CRT_SECURE_NO_WARNINGS -#include #include "inmost_solver.h" #if defined(USE_SOLVER) #include "solver_mlmtiluc2.hpp" #include #include +#include //#define REPORT_ILU //#undef REPORT_ILU //#define REPORT_ILU_PROGRESS @@ -597,7 +597,7 @@ public: double tlfactor, tlrescale, tlreorder, tlreassamble; ttotal = Timer(); - (*Alink).Save("M.mtx"); + //(*Alink).Save("M.mtx"); //calculate number of nonzeros nzA = 0; @@ -3242,7 +3242,7 @@ public: Li = LineIndecesU[Li]; } LF_Address[k].last = static_cast(LF_Entries.size()); - assert(std::is_sorted(LF_Entries.begin()+LF_Address[k].first,LF_Entries.end())); + //assert(std::is_sorted(LF_Entries.begin()+LF_Address[k].first,LF_Entries.end())); /////////////////////////////////////////////////////////////////////////////////// // clean linked list // /////////////////////////////////////////////////////////////////////////////////// diff --git a/Source/Solver/solver_inner/solver_mptiluc/solver_mtiluc2.cpp b/Source/Solver/solver_inner/solver_mptiluc/solver_mtiluc2.cpp index 24f5e789c112863cf71b3822f44df0d2a1e59f63..4a561b655638babf516be9e9ca344497857641f1 100644 --- a/Source/Solver/solver_inner/solver_mptiluc/solver_mtiluc2.cpp +++ b/Source/Solver/solver_inner/solver_mptiluc/solver_mtiluc2.cpp @@ -1,10 +1,11 @@ #define _CRT_SECURE_NO_WARNINGS -#include + #include "inmost_solver.h" #if defined(USE_SOLVER) #include "solver_mtiluc2.hpp" #include #include +#include //#define REPORT_ILU //#undef REPORT_ILU //#define REPORT_ILU_PROGRESS @@ -3347,8 +3348,16 @@ swap_algorithm: level_interval[k].last = level_interval[k].first + level_size[k]; } */ - condestL = NuL; - condestU = NuU; + condestL = NuL; + condestU = NuU; + //for partition of unity + /* + info->PrepareVector(div); + std::fill(div.Begin(),div.End(),0); + for(k = mobeg; k < moend; k++) div[k] = 1.0; + info->Accumulate(div); + for(k = mobeg; k < moend; k++) div[k] = 1.0/div[k]; + */ return true; } bool MTILUC_preconditioner::Finalize() @@ -3384,11 +3393,12 @@ swap_algorithm: { assert(&input != &output); // + INMOST_DATA_ENUM_TYPE mobeg, moend, k; #if defined(USE_OMP) #pragma omp single #endif { - INMOST_DATA_ENUM_TYPE k, mobeg, moend, vbeg, vend; + INMOST_DATA_ENUM_TYPE vbeg, vend; info->GetOverlapRegion(info->GetRank(), mobeg, moend); info->GetVectorRegion(vbeg, vend); @@ -3441,6 +3451,14 @@ swap_algorithm: for (k = moend; k < vend; ++k) output[k] = 0; } info->Accumulate(output); + /* +#if defined(USE_OMP) +#pragma omp single +#endif + { + for (k = mobeg; k < moend; ++k) output[k] *= div[k]; + } + */ return true; } bool MTILUC_preconditioner::ReplaceMAT(Sparse::Matrix & A){ if (isInitialized()) Finalize(); Alink = &A; return true; } diff --git a/Source/Solver/solver_inner/solver_mptiluc/solver_mtiluc2.hpp b/Source/Solver/solver_inner/solver_mptiluc/solver_mtiluc2.hpp index 9fbc11802f536c185f7daf63c303e1de802f2432..ba298a45245198b5d8a55b1be6008fb9db4186df 100644 --- a/Source/Solver/solver_inner/solver_mptiluc/solver_mtiluc2.hpp +++ b/Source/Solver/solver_inner/solver_mptiluc/solver_mtiluc2.hpp @@ -23,6 +23,7 @@ class MTILUC_preconditioner : public Method } row_col; typedef dynarray levels_t; + //Sparse::Vector div; std::vector LU_Entries, B_Entries; interval LU_Diag; interval U_Address, L_Address, B_Address; diff --git a/Source/Solver/sparse.cpp b/Source/Solver/sparse.cpp index 6eab46e31016a9954e504bc6ae61d3f850d46974..ee26d8c0826a83f27b0d937fffb7bd324c3de0c2 100644 --- a/Source/Solver/sparse.cpp +++ b/Source/Solver/sparse.cpp @@ -1,5 +1,6 @@ #define _CRT_SECURE_NO_WARNINGS #include "inmost_sparse.h" +#include "inmost_dense.h" #include #include @@ -62,42 +63,51 @@ namespace INMOST INMOST_DATA_REAL_TYPE & RowMerger::operator[] (INMOST_DATA_ENUM_TYPE pos) { - pos = MapIndex(pos); - if( LinkedList[pos+1].first != UNDEF ) return LinkedList[pos+1].second; + INMOST_DATA_ENUM_TYPE map = MapIndex(pos); + if( map == ENUMUNDEF ) + return Nonlocal[pos]; + if( LinkedList[map+1].first != UNDEF ) return LinkedList[map+1].second; else { INMOST_DATA_ENUM_TYPE index = LinkedList.get_interval_beg(), next; if( Sorted ) { next = index; - while(next < pos+1) + while(next < map+1) { index = next; next = LinkedList[index].first; } - assert(index < pos+1); - assert(pos+1 < next); + assert(index < map+1); + assert(map+1 < next); ++Nonzeros; - LinkedList[index].first = pos+1; - LinkedList[pos+1].first = next; - return LinkedList[pos+1].second; + LinkedList[index].first = map+1; + LinkedList[map+1].first = next; + return LinkedList[map+1].second; } else { INMOST_DATA_ENUM_TYPE index = LinkedList.get_interval_beg(); ++Nonzeros; - LinkedList[pos+1].first = LinkedList[index].first; - LinkedList[index].first = pos+1; - return LinkedList[pos+1].second; + LinkedList[map+1].first = LinkedList[index].first; + LinkedList[index].first = map+1; + return LinkedList[map+1].second; } } } INMOST_DATA_REAL_TYPE RowMerger::operator[] (INMOST_DATA_ENUM_TYPE pos) const { - pos = MapIndex(pos); - if( LinkedList[pos+1].first != UNDEF ) return LinkedList[pos+1].second; - else throw -1; + INMOST_DATA_ENUM_TYPE map = MapIndex(pos); + if( map == ENUMUNDEF ) + { + std::map::const_iterator it = Nonlocal.find(pos); + if( it != Nonlocal.end() ) + return it->second; + } + else + if( LinkedList[map+1].first != UNDEF ) return LinkedList[map+1].second; + throw -1; } RowMerger::RowMerger(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, bool Sorted) @@ -125,6 +135,7 @@ namespace INMOST LinkedList.begin()->first = EOL; Nonzeros = 0; Sorted = _Sorted; + Nonlocal.clear(); } void RowMerger::Resize(INMOST_DATA_ENUM_TYPE interval_begin, INMOST_DATA_ENUM_TYPE interval_end, const std::vector & Pre, const std::vector & Post, bool _Sorted) @@ -132,6 +143,7 @@ namespace INMOST NonlocalPre = Pre; NonlocalPost = Post; Resize(interval_begin,interval_end,_Sorted); + Nonlocal.clear(); } #if defined(USE_SOLVER) void RowMerger::Resize(const Matrix & A, bool _Sorted) @@ -154,8 +166,9 @@ namespace INMOST NonlocalPre.clear(); NonlocalPre.insert(NonlocalPre.end(),Pre.begin(),Pre.end()); NonlocalPost.clear(); - NonlocalPre.insert(NonlocalPost.end(),Post.begin(),Post.end()); + NonlocalPost.insert(NonlocalPost.end(),Post.begin(),Post.end()); Resize(mbeg,mend,_Sorted); + Nonlocal.clear(); } RowMerger::RowMerger(const Matrix & A, bool Sorted) : Sorted(Sorted), Nonzeros(0) @@ -173,6 +186,7 @@ namespace INMOST NonlocalPre = Pre; NonlocalPost = Post; Resize(IntervalBeg,IntervalEnd,Sorted); + Nonlocal.clear(); } @@ -180,17 +194,25 @@ namespace INMOST { if( pos < IntervalBeg ) { - assert(!NonlocalPre.empty()); //there are indices provided + return ENUMUNDEF; + //assert(!NonlocalPre.empty()); //there are indices provided std::vector< INMOST_DATA_ENUM_TYPE >::const_iterator search = std::lower_bound(NonlocalPre.begin(),NonlocalPre.end(),pos); - assert(*search == pos); //is there such index? - return static_cast(IntervalBeg - NonlocalPre.size() + static_cast(search - NonlocalPre.begin())); + //assert(*search == pos); //is there such index? + if( search != NonlocalPre.end() && *search == pos ) + return static_cast(IntervalBeg - NonlocalPre.size() + static_cast(search - NonlocalPre.begin())); + else + return ENUMUNDEF; } else if( pos >= IntervalEnd ) { - assert(!NonlocalPost.empty()); //there are indices provided + return ENUMUNDEF; + //assert(!NonlocalPost.empty()); //there are indices provided std::vector< INMOST_DATA_ENUM_TYPE >::const_iterator search = std::lower_bound(NonlocalPost.begin(),NonlocalPost.end(),pos); - assert(*search == pos); //is there such index? - return IntervalEnd + static_cast(search-NonlocalPost.begin()); + //assert(*search == pos); //is there such index? + if( search != NonlocalPost.end() && *search == pos ) + return IntervalEnd + static_cast(search-NonlocalPost.begin()); + else + return ENUMUNDEF; } return pos; } @@ -198,9 +220,15 @@ namespace INMOST INMOST_DATA_ENUM_TYPE RowMerger::UnmapIndex(INMOST_DATA_ENUM_TYPE pos) const { if( pos < IntervalBeg ) + { + assert(pos >= IntervalBeg-NonlocalPre.size()); return NonlocalPre[pos+NonlocalPre.size()-IntervalBeg]; + } else if( pos >= IntervalEnd ) + { + assert(pos < IntervalEnd+NonlocalPost.size()); return NonlocalPost[pos-IntervalEnd]; + } return pos; } @@ -216,6 +244,7 @@ namespace INMOST LinkedList[i].second = 0.0; i = j; } + Nonlocal.clear(); Nonzeros = 0; } @@ -225,6 +254,23 @@ namespace INMOST if( Sorted && PreSortRow ) std::sort(r.Begin(),r.End()); PushRow(coef,r); } + + void RowMerger::Multiply(INMOST_DATA_REAL_TYPE coef) + { + INMOST_DATA_ENUM_TYPE i = LinkedList.begin()->first; + while( i != EOL ) + { + LinkedList[i].second *= coef; + i = LinkedList[i].first; + } + if( !Nonlocal.empty() ) + { + std::map< INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE >::iterator it; + for(it = Nonlocal.begin(); it != Nonlocal.end(); ++it) + it->second *= coef; + } + } + void RowMerger::PushRow(INMOST_DATA_REAL_TYPE coef, const Row & r) { @@ -235,12 +281,17 @@ namespace INMOST while( it != r.End() ) { INMOST_DATA_ENUM_TYPE pos = it->first; - pos = MapIndex(pos); - LinkedList[index].first = pos+1; - LinkedList[pos+1].first = EOL; - LinkedList[pos+1].second = it->second*coef; - index = pos+1; - ++Nonzeros; + INMOST_DATA_ENUM_TYPE map = MapIndex(pos); + if( map == ENUMUNDEF ) + Nonlocal[pos] = it->second*coef; + else + { + LinkedList[index].first = map+1; + LinkedList[map+1].first = EOL; + LinkedList[map+1].second = it->second*coef; + index = map+1; + ++Nonzeros; + } jt = it; ++it; assert(!Sorted || it == r.End() || jt->first < it->first); @@ -260,30 +311,35 @@ namespace INMOST while( it != r.End() ) { INMOST_DATA_ENUM_TYPE pos = it->first; - pos = MapIndex(pos); - if( LinkedList[pos+1].first != UNDEF ) - LinkedList[pos+1].second += coef*it->second; - else if( Sorted ) + INMOST_DATA_ENUM_TYPE map = MapIndex(pos); + if( map == ENUMUNDEF ) + Nonlocal[pos] += coef*it->second; + else { - next = index; - while(next < pos+1) + if( LinkedList[map+1].first != UNDEF ) + LinkedList[map+1].second += coef*it->second; + else if( Sorted ) { - index = next; - next = LinkedList[index].first; + next = index; + while(next < map+1) + { + index = next; + next = LinkedList[index].first; + } + assert(index < map+1); + assert(map+1 < next); + LinkedList[index].first = map+1; + LinkedList[map+1].first = next; + LinkedList[map+1].second = coef*it->second; + ++Nonzeros; + } + else + { + LinkedList[map+1].first = LinkedList[index].first; + LinkedList[map+1].second = coef*it->second; + LinkedList[index].first = map+1; + ++Nonzeros; } - assert(index < pos+1); - assert(pos+1 < next); - LinkedList[index].first = pos+1; - LinkedList[pos+1].first = next; - LinkedList[pos+1].second = coef*it->second; - ++Nonzeros; - } - else - { - LinkedList[pos+1].first = LinkedList[index].first; - LinkedList[pos+1].second = coef*it->second; - LinkedList[index].first = pos+1; - ++Nonzeros; } jt = it; ++it; @@ -294,7 +350,7 @@ namespace INMOST void RowMerger::RetriveRow(Row & r) { - r.Resize(Nonzeros); + r.Resize(Nonzeros+Nonlocal.size()); INMOST_DATA_ENUM_TYPE i = LinkedList.begin()->first, k = 0; while( i != EOL ) { @@ -306,6 +362,16 @@ namespace INMOST } i = LinkedList[i].first; } + if( !Nonlocal.empty() ) + { + std::map< INMOST_DATA_ENUM_TYPE, INMOST_DATA_REAL_TYPE >::iterator it; + for(it = Nonlocal.begin(); it != Nonlocal.end(); ++it) if( it->second ) + { + r.GetIndex(k) = it->first; + r.GetValue(k) = it->second; + ++k; + } + } r.Resize(k); } ////////class HessianRow @@ -826,6 +892,15 @@ namespace INMOST output << rhs.rdbuf(); #endif } + + INMOST::Matrix Vector::operator [](const INMOST::AbstractMatrix & rows) const + { + INMOST::Matrix ret(rows.Rows(),rows.Cols()); + for(INMOST_DATA_ENUM_TYPE i = 0; i < rows.Rows(); ++i) + for(INMOST_DATA_ENUM_TYPE j = 0; j < rows.Cols(); ++j) + ret(i,j) = data[rows(i,j)]; + return ret; + } ////////class Matrix void Matrix::MatVec(INMOST_DATA_REAL_TYPE alpha, Vector & x, INMOST_DATA_REAL_TYPE beta, Vector & out) const //y = alpha*A*x + beta * y