[cig-commits] r1374 - in branches/s-wang: . for_deal.II/source/lac source/simulator

[s-wang at dealii.org]

Author: s-wang
Date: 2012-11-16 16:20:33 -0700 (Fri, 16 Nov 2012)
New Revision: 1374

Modified:
   branches/s-wang/for_deal.II/source/lac/constraint_matrix.cc
   branches/s-wang/simple_2d_shell.prm
   branches/s-wang/source/simulator/initial_conditions.cc
Log:
corrected error in ConstraintMatrix::distribute (PETScWrappers::MPI::BlockVector &vec).


Modified: branches/s-wang/for_deal.II/source/lac/constraint_matrix.cc
===================================================================
--- branches/s-wang/for_deal.II/source/lac/constraint_matrix.cc	2012-11-16 13:16:35 UTC (rev 1373)
+++ branches/s-wang/for_deal.II/source/lac/constraint_matrix.cc	2012-11-16 23:20:33 UTC (rev 1374)
@@ -2139,10 +2139,111 @@
 void
 ConstraintMatrix::distribute (PETScWrappers::MPI::BlockVector &vec) const   // modified by shuqiangwang
 {
-  Assert (sorted==true, ExcMatrixIsClosed());
-//  AssertThrow (false, ExcNotImplemented());
-  for(unsigned int i=0; i<vec.n_blocks(); i++)
-	  distribute(vec.block(i));
+//  Assert (sorted==true, ExcMatrixIsClosed());
+//  for(unsigned int i=0; i<vec.n_blocks(); i++)
+//	  distribute(vec.block(i));
+
+	Assert (sorted==true, ExcMatrixIsClosed());
+
+//	IndexSet my_indices (vec.size());
+//	IndexSet local_range_is (vec.size());
+	std::vector<IndexSet> my_indices(vec.n_blocks());
+	std::vector<IndexSet> local_range_is(vec.n_blocks());
+	std::vector<unsigned int> block_sizes(vec.n_blocks());
+
+	for(unsigned int block=0; block<vec.n_blocks(); block++)
+	{
+		my_indices[block].set_size(vec.block(block).size());
+		local_range_is[block].set_size(vec.block(block).size());
+//		block_sizes[block] = vec(block).size();
+	}
+
+	for (unsigned int block=0; block<vec.n_blocks(); ++block)
+	{
+		typedef std::vector<ConstraintLine>::const_iterator constraint_iterator;
+		ConstraintLine index_comparison;
+		index_comparison.line = vec.block(block).local_range().first
+				+vec.get_block_indices().block_start(block);
+		const constraint_iterator begin_my_constraints =
+				Utilities::lower_bound (lines.begin(),lines.end(),index_comparison);
+
+		index_comparison.line = vec.block(block).local_range().second
+				+vec.get_block_indices().block_start(block);
+
+		const constraint_iterator end_my_constraints
+		= Utilities::lower_bound(lines.begin(),lines.end(),index_comparison);
+
+		// Here we search all the indices that we
+		// need to have read-access to - the local
+		// nodes and all the nodes that the
+		// constraints indicate. No caching done
+		// yet. would need some more clever data
+		// structures for doing that.
+		const std::pair<unsigned int, unsigned int>
+		local_range = vec.block(block).local_range();
+		local_range_is[block].add_range(local_range.first, local_range.second);
+		block_sizes[block] = local_range.second - local_range.first;
+
+		my_indices[block].add_range (local_range.first, local_range.second);
+
+		std::set<unsigned int> individual_indices;
+		for (constraint_iterator it = begin_my_constraints;
+				it != end_my_constraints; ++it)
+			for (unsigned int i=0; i<it->entries.size(); ++i)
+				if ((it->entries[i].first < local_range.first)
+						||
+						(it->entries[i].first >= local_range.second))
+					individual_indices.insert (it->entries[i].first);
+
+		my_indices[block].add_indices (individual_indices.begin(),
+				individual_indices.end());
+	}
+
+	// create a vector and import those indices
+//	PETScWrappers::MPI::BlockVector ghost_vec (vec.get_mpi_communicator(),
+//			local_range_is,
+//			my_indices);
+//	ghost_vec = vec;
+//	ghost_vec.update_ghost_values();
+	PETScWrappers::MPI::BlockVector ghost_vec;
+	ghost_vec.reinit(block_sizes,vec.get_mpi_communicator());
+	for (unsigned int block=0; block<vec.n_blocks(); ++block)
+		ghost_vec.block(block).reinit(vec.get_mpi_communicator(),local_range_is[block],my_indices[block]);
+	ghost_vec.collect_sizes();
+
+	ghost_vec = vec;
+	ghost_vec.update_ghost_values();
+
+	for (unsigned int block=0; block<vec.n_blocks(); ++block)
+	{
+		typedef std::vector<ConstraintLine>::const_iterator constraint_iterator;
+		ConstraintLine index_comparison;
+		index_comparison.line = vec.block(block).local_range().first
+				+vec.get_block_indices().block_start(block);
+		const constraint_iterator begin_my_constraints =
+				Utilities::lower_bound (lines.begin(),lines.end(),index_comparison);
+
+		index_comparison.line = vec.block(block).local_range().second
+				+vec.get_block_indices().block_start(block);
+
+		const constraint_iterator end_my_constraints
+		= Utilities::lower_bound(lines.begin(),lines.end(),index_comparison);
+
+		for (constraint_iterator it = begin_my_constraints;
+				it != end_my_constraints; ++it)
+		{
+			// fill entry in line
+			// next_constraint.line by adding the
+			// different contributions
+			double new_value = it->inhomogeneity;
+			for (unsigned int i=0; i<it->entries.size(); ++i)
+				new_value += (ghost_vec(it->entries[i].first) *
+						it->entries[i].second);
+			vec(it->line) = new_value;
+		}
+		vec.block(block).compress(::dealii::VectorOperation::insert);
+	}
+
 }
 
 #endif

Modified: branches/s-wang/simple_2d_shell.prm
===================================================================
--- branches/s-wang/simple_2d_shell.prm	2012-11-16 13:16:35 UTC (rev 1373)
+++ branches/s-wang/simple_2d_shell.prm	2012-11-16 23:20:33 UTC (rev 1374)
@@ -247,11 +247,11 @@
 
   # The number of adaptive refinement steps performed after initial global
   # refinement but while still within the first time step.
-  set Initial adaptive refinement        = 3
+  set Initial adaptive refinement        = 0
 
   # The number of global refinement steps performed on the initial coarse
   # mesh, before the problem is first solved there.
-  set Initial global refinement          = 4
+  set Initial global refinement          = 3
 
   # The fraction of cells with the largest error that should be flagged for
   # refinement.

Modified: branches/s-wang/source/simulator/initial_conditions.cc
===================================================================
--- branches/s-wang/source/simulator/initial_conditions.cc	2012-11-16 13:16:35 UTC (rev 1373)
+++ branches/s-wang/source/simulator/initial_conditions.cc	2012-11-16 23:20:33 UTC (rev 1374)
@@ -202,6 +202,7 @@
         constraints.distribute (system_tmp);
         constraints.print(std::cout);
 //        exit(0);
+
         system_tmp.compress();
         system_tmp.print(std::cout,7,false,false);
         exit(0);